|
 轉(zhuǎn)自 http://my.oschina.net/cvnote/blog/166980
Eigen自帶的稀疏矩陣分解功能包括LDLt、LLt分解（即Cholesky分解，這個功能是LGPL許可，不是Eigen的MPL許可）、LU分解、QR分解（這個是3.2版本之后正式Release的）、共軛梯度解矩陣等。另外還提供了到第三方稀疏矩陣庫的C++接口，包括著名的SuiteSparse系列（這個系列非常強大，有機會要好好研究一下）的SparseQR、UmfPack等。(歡迎訪問計算機視覺研究筆記和關(guān)注新浪微博@cvnote )
基本稀疏矩陣操作
Eigen中使用 Eigen::Triplet<Scalar>來記錄一個非零元素的行、列、值，填充一個稀疏矩陣，只需要將所有表示非零元素的Triplet放在一個 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
}
;
// 非零元素的行號
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
]
)
;
// 填充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
;
// 可以各種運算
std
::
cout
<<
"A * A = \n"
<<
A *
A
<<
std
::
endl
;
return
0
;
}
|
用Eigen進行矩陣分解
首先復(fù)習(xí)一下Cholesky（LLt）、QR和LU分解，說的不對的地方歡迎數(shù)學(xué)大牛和數(shù)值計算大牛來指教和補充。一般來講LLt分解可以理解成給矩陣開平方，類比于開平方一般針對正數(shù)而言，LLt分解則限定矩陣需為正定的 Hermitian矩陣（自共軛矩陣，即對稱的實數(shù)矩陣或?qū)ΨQ元素共軛的復(fù)數(shù)矩陣）。LU分解則稍微放松一點，可用于一般的方陣（順便提一句LU分解是圖靈發(fā)明的）。QR則可用于一般矩陣，結(jié)果也是最穩(wěn)定的。分解算法的效率，三者都是O(n^3)的，具體系數(shù)三者大概是Cholesky:LU:QR=1:2:4。Google可以找到很多相關(guān)資料，比如我看了這個。
下面試一下用Eigen自帶的Eigen::SparseQR 進行我最喜歡的QR分解（其實我更喜歡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
;
// 一個QR分解的實例
SparseQR
<
SparseMatrix
<
double
>
,
AMDOrdering
<
int
>
>
qr
;
// 計算分解
qr
.
compute
(
A
)
;
// 求一個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分解解了一個系數(shù)矩陣，A和上面的例子是一樣的，注意到這個Ax=b其實是沒有確定解的，A(1:3, 2:3)是over determined的，剩下的部分又是非滿秩under determined的，這個QR分解對于A(1:3, 2:3)給了最小二乘解，其他位返回了0。
另一個注意的地方就是 SparseQR<SparseMatrix<double>, AMDOrdering<int> >的第二個模板參數(shù)，是一個矩陣重排列（ordering）的方法，為什么要重排列呢，wikipedia的LU分解詞條給了一個例子可以大概解釋一下，某些矩陣沒有重排直接分解可能會失敗。Eigen提供了三種重排列方法，參見OrderingMethods Module。關(guān)于矩陣重排列的細節(jié)求數(shù)值計算牛人指點！我一般就隨便選一個填進去了>_<。
除了解方程，這個QR實例也可以用下面代碼返回Q和R矩陣：
|
|
SparseMatrix
<
double
>
Q
,
R
;
qr
.
matrixQ
(
)
.
evalTo
(
Q
)
;
R
=
qr
.
matrixR
(
)
;
|
注意到Q和R的返回方法不一樣，猜測是因為 matrixQ() 成員好像是沒有完整保存Q矩陣（元素太多？）。
用 Eigen::SPQR 調(diào)用SuiteSparseQR進行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
;
// 一個QR分解的實例
SPQR
<
SparseMatrix
<
double
>
>
qr
;
// 計算分解
qr
.
compute
(
A
)
;
// 求一個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的話，會覺得這個接口真心好用多了。編譯這個程序除了spqr庫還需要鏈接blas庫、lapack庫、cholmod庫（SuiteSparse的另一組件），有一點麻煩。比如我在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是因為上面代碼偷懶用了emplace_back ，和矩陣庫沒關(guān)系。
一些效率的經(jīng)驗
SuiteSparseQR畢竟實現(xiàn)更好一些，我的一些經(jīng)驗是比自帶Eigen::SparseQR快50%左右吧。
相關(guān)文章
C++矩陣運算庫推薦
歡迎訪問計算機視覺研究筆記和關(guān)注新浪微博@cvnote
|
