|
 轉(zhuǎn)自 http://my.oschina.net/cvnote/blog/166980
Eigen自帶的稀疏矩陣分解功能包括LDLt、LLt分解（即Cholesky分解，這個(gè)功能是LGPL許可，不是Eigen的MPL許可）、LU分解、QR分解（這個(gè)是3.2版本之后正式Release的）、共軛梯度解矩陣等。另外還提供了到第三方稀疏矩陣庫的C++接口，包括著名的SuiteSparse系列（這個(gè)系列非常強(qiáng)大，有機(jī)會(huì)要好好研究一下）的SparseQR、UmfPack等。(歡迎訪問計(jì)算機(jī)視覺研究筆記和關(guān)注新浪微博@cvnote )
基本稀疏矩陣操作
Eigen中使用 Eigen::Triplet<Scalar>來記錄一個(gè)非零元素的行、列、值，填充一個(gè)稀疏矩陣，只需要將所有表示非零元素的Triplet放在一個(gè) std::vector中即可傳入即可。除了求逆等功能外，Eigen::SparseMatrix 有和 Eigen::Matrix幾乎一樣的各種成員操作函數(shù)，并且可以方便混用。
比如這樣：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
#include <iostream>
#include "eigen/Eigen/Eigen"
int
main
(
)
{
// 矩陣：
// 0 6.1 0 0
// 0 0 7.2 0
// 0 0 8.3 0
// 0 0 0 0
using
namespace
Eigen
;
SparseMatrix
<
double
>
A
(
4
,
4
)
;
std
::
vector
<
Triplet
<
double
>
>
triplets
;
int
r
[
3
]
=
{
0
,
1
,
2
}
;
// 非零元素的行號(hào)
int
c
[
3
]
=
{
1
,
2
,
2
}
;
// 非零元素的列號(hào)
double
val
[
3
]
=
{
6.1
,
7.2
,
8.3
}
;
// 非零元素的值
for
(
int
i
=
0
;
i
<
3
;
++
i
)
triplets
.
emplace_back
(
r
[
i
]
,
c
[
i
]
,
val
[
i
]
)
;
// 填充Triplet
A
.
setFromTriplets
(
triplets
.
begin
(
)
,
triplets
.
end
(
)
)
;
// 初始化系數(shù)矩陣
std
::
cout
<<
"A = "
<<
A
<<
std
::
endl
;
MatrixXd
B
=
A
;
// 可以和普通稠密矩陣方便轉(zhuǎn)換
std
::
cout
<<
"B = \n"
<<
B
<<
std
::
endl
;
std
::
cout
<<
"A * B = \n"
<<
A *
B
<<
std
::
endl
;
// 可以各種運(yùn)算
std
::
cout
<<
"A * A = \n"
<<
A *
A
<<
std
::
endl
;
return
0
;
}
|
用Eigen進(jìn)行矩陣分解
首先復(fù)習(xí)一下Cholesky（LLt）、QR和LU分解，說的不對的地方歡迎數(shù)學(xué)大牛和數(shù)值計(jì)算大牛來指教和補(bǔ)充。一般來講LLt分解可以理解成給矩陣開平方，類比于開平方一般針對正數(shù)而言，LLt分解則限定矩陣需為正定的 Hermitian矩陣（自共軛矩陣，即對稱的實(shí)數(shù)矩陣或?qū)ΨQ元素共軛的復(fù)數(shù)矩陣）。LU分解則稍微放松一點(diǎn)，可用于一般的方陣（順便提一句LU分解是圖靈發(fā)明的）。QR則可用于一般矩陣，結(jié)果也是最穩(wěn)定的。分解算法的效率，三者都是O(n^3)的，具體系數(shù)三者大概是Cholesky:LU:QR=1:2:4。Google可以找到很多相關(guān)資料，比如我看了這個(gè)。
下面試一下用Eigen自帶的Eigen::SparseQR 進(jìn)行我最喜歡的QR分解（其實(shí)我更喜歡SVD）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
#include <iostream>
#include "eigen/Eigen/Eigen"
#include "eigen/Eigen/SparseQR"
int
main
(
)
{
using
namespace
Eigen
;
SparseMatrix
<
double
>
A
(
4
,
4
)
;
std
::
vector
<
Triplet
<
double
>
>
triplets
;
// 初始化非零元素
int
r
[
3
]
=
{
0
,
1
,
2
}
;
int
c
[
3
]
=
{
1
,
2
,
2
}
;
double
val
[
3
]
=
{
6.1
,
7.2
,
8.3
}
;
for
(
int
i
=
0
;
i
<
3
;
++
i
)
triplets
.
emplace_back
(
r
[
i
]
,
c
[
i
]
,
val
[
i
]
)
;
// 初始化稀疏矩陣
A
.
setFromTriplets
(
triplets
.
begin
(
)
,
triplets
.
end
(
)
)
;
std
::
cout
<<
"A = \n"
<<
A
<<
std
::
endl
;
// 一個(gè)QR分解的實(shí)例
SparseQR
<
SparseMatrix
<
double
>
,
AMDOrdering
<
int
>
>
qr
;
// 計(jì)算分解
qr
.
compute
(
A
)
;
// 求一個(gè)A x = b
Vector4d
b
(
1
,
2
,
3
,
4
)
;
Vector4d
x
=
qr
.
solve
(
b
)
;
std
::
cout
<<
"x = \n"
<<
x
;
std
::
cout
<<
"A x = \n"
<<
A *
x
;
return
0
;
}
|
這樣就用QR分解解了一個(gè)系數(shù)矩陣，A和上面的例子是一樣的，注意到這個(gè)Ax=b其實(shí)是沒有確定解的，A(1:3, 2:3)是over determined的，剩下的部分又是非滿秩under determined的，這個(gè)QR分解對于A(1:3, 2:3)給了最小二乘解，其他位返回了0。
另一個(gè)注意的地方就是 SparseQR<SparseMatrix<double>, AMDOrdering<int> >的第二個(gè)模板參數(shù)，是一個(gè)矩陣重排列（ordering）的方法，為什么要重排列呢，wikipedia的LU分解詞條給了一個(gè)例子可以大概解釋一下，某些矩陣沒有重排直接分解可能會(huì)失敗。Eigen提供了三種重排列方法，參見OrderingMethods Module。關(guān)于矩陣重排列的細(xì)節(jié)求數(shù)值計(jì)算牛人指點(diǎn)！我一般就隨便選一個(gè)填進(jìn)去了>_<。
除了解方程，這個(gè)QR實(shí)例也可以用下面代碼返回Q和R矩陣：
|
|
SparseMatrix
<
double
>
Q
,
R
;
qr
.
matrixQ
(
)
.
evalTo
(
Q
)
;
R
=
qr
.
matrixR
(
)
;
|
注意到Q和R的返回方法不一樣，猜測是因?yàn)?nbsp;matrixQ() 成員好像是沒有完整保存Q矩陣（元素太多？）。
用 Eigen::SPQR 調(diào)用SuiteSparseQR進(jìn)行QR分解
SuiteSparseQR效率很高，但是C風(fēng)格接口比較不好用，Eigen提供了 Eigen::SPQR 的接口封裝比如和上面同樣的程序可以這樣寫：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
#include <iostream>
#include "eigen/Eigen/Eigen"
#include "eigen/Eigen/SPQRSupport"
int
main
(
)
{
using
namespace
Eigen
;
SparseMatrix
<
double
>
A
(
4
,
4
)
;
std
::
vector
<
Triplet
<
double
>
>
triplets
;
// 初始化非零元素
int
r
[
3
]
=
{
0
,
1
,
2
}
;
int
c
[
3
]
=
{
1
,
2
,
2
}
;
double
val
[
3
]
=
{
6.1
,
7.2
,
8.3
}
;
for
(
int
i
=
0
;
i
<
3
;
++
i
)
triplets
.
emplace_back
(
r
[
i
]
,
c
[
i
]
,
val
[
i
]
)
;
// 初始化稀疏矩陣
A
.
setFromTriplets
(
triplets
.
begin
(
)
,
triplets
.
end
(
)
)
;
std
::
cout
<<
"A = \n"
<<
A
<<
std
::
endl
;
// 一個(gè)QR分解的實(shí)例
SPQR
<
SparseMatrix
<
double
>
>
qr
;
// 計(jì)算分解
qr
.
compute
(
A
)
;
// 求一個(gè)A x = b
Vector4d
b
(
1
,
2
,
3
,
4
)
;
Vector4d
x
=
qr
.
solve
(
b
)
;
std
::
cout
<<
"x = \n"
<<
x
;
std
::
cout
<<
"A x = \n"
<<
A *
x
;
return
0
;
}
|
如果你有用過SuiteSparseQR的話，會(huì)覺得這個(gè)接口真心好用多了。編譯這個(gè)程序除了spqr庫還需要鏈接blas庫、lapack庫、cholmod庫（SuiteSparse的另一組件），有一點(diǎn)麻煩。比如我在ubuntu，使用 apt-get install libsuitesparse-* 安裝了suitesparse頭文件到/usr/include/suitesparse 目錄，使用如下命令編譯。
|
|
g
++
spqr
.
cpp
-
std
=
c
++
11
-
I
/
usr
/
include
/
suitesparse
-
lcholmod
-
lspqr
-
llapack
-
lblas
|
注意lapack要在blas前面，spqr要在lapack前面。用了c++11是因?yàn)樯厦娲a偷懶用了emplace_back ，和矩陣庫沒關(guān)系。
一些效率的經(jīng)驗(yàn)
SuiteSparseQR畢竟實(shí)現(xiàn)更好一些，我的一些經(jīng)驗(yàn)是比自帶Eigen::SparseQR快50%左右吧。
相關(guān)文章
C++矩陣運(yùn)算庫推薦
歡迎訪問計(jì)算機(jī)視覺研究筆記和關(guān)注新浪微博@cvnote
|
