图像处理中往往会遇到大型稀疏矩阵的特征值分解问题,如果采用一般的密集矩阵求解方法,对内存要求比较大。所以在稀疏矩阵的基础上进行特征值分解,更加合理。opencv2.49版本中还没有对稀疏矩阵的特征值分解函数,所以自己对此类问题进行了网上搜索。调研了ietl、lapack、Armadillo、Eigen这类库。文章先对各种库进行总结,后面在详细的介绍Armadillo库在ubuntu14.04系统上的安装方法,最后给出一个利用Armadillo进行稀疏矩阵特征值分解的例子。
1. 各种库的总结
ietl
网站地址:
http://www.comp-phys.org/software/ietl/
该库可以解决稀疏矩阵的特征值分解问题,但是其目前其不经已经没有更新,而且在求解多个特征值时,只有第一个值是正确的,其余的均为0. 这个是其自身算法稳定性问题,所以自己舍弃了这一个工具包
lapack
网站地址:
http://www.netlib.org/lapack/#_lapack_version_3_5_0
该库解决了密集矩阵的分解问题,但是稀疏矩阵的分解并不是它的优势。
Armadillo
网站地址:
http://arma.sourceforge.net/docs.html
该库目前还在更新,并且函数仿照matlab使用习惯,使用方便,容易上手,并且其对稀疏矩阵的运算支持的很全面。它既支持密集矩阵的求解也支持稀疏矩阵的求解。密集矩阵求解主要基于lapack和blas等库,稀疏矩阵求解主要基于superlu和arpack等库。
Eigen
网站地址:
http://eigen.tuxfamily.org/dox/index.html
该库目前也在更新中,其对密集矩阵的支持比较全面,但是对稀疏矩阵的支持比较欠缺,其没有针对稀疏矩阵特征值分解的函数,有关于LU、和QR分解的功能。
2.Armadillo 库的安装
Armadillo库依赖LAPACK、 BLAS、ARPACK 、SuperLU、ATLAS、OpenBLAS等库。如果只需要对密集矩阵进行处理,安装LAPACK、 BLAS库即可,如果需要对稀疏矩阵进行求解则需要包括ARPACK 、SuperLU。OpenBLAS是BLAS的进化版本,它支持多线程,所以为了提高速度可以采用OpenBLAS代替BLAS。
1) LAPACK包和BLAS包的安装
因为LAPACK包安装时可以同时生成BLAS包,所以这两个一起安装。这里使用的是3.5版本。
将make.inc.example文件改名为make.inc,对该文件进行修改,修改如下,加上lapackelib这一项,这个库后面可能会用到。修改后会生成,执行make 命令生成libblas.a liblapack.a liblapacke.a等库文件,然后将生成的所有库都复制到/usr/local/lib中。
2)ATLAS包的安装
3)ARPACK包的安装
按照README内容对ARmake.inc进行修改
这是对用到的源文件路径进行设置。
编译器部分设置如下,主要是“FFLAGS”和“MAKE”选项的设置。
修改完成后可以直接make lib进行编译
编译完成后,将libarpack.a复制到 /usr/local/lib目录中
4)openblas包安装
按照 readme.md文件进行安装。
直接make编译,make install安装。安装文件默认在./OpenBLAS-0.2.14/lib目录中.
./OpenBLAS-0.2.14/lib/lib 中所有库文件复制到/usr/local/lib目录中,./OpenBLAS-0.2.14/lib/include中的所有文件,复制到/usr/local/include/openblas目录下
5)superlu包安装
根据readme文档进行安装,这里采用4.3版本,因为Armadillo说明文档中提示只支持这一版本。
因为暂时用不到matlab相关的接口,所以暂时到此结束。
6) libf2c.a库的安装
如果系统中没有libf2c.a这一个库函数,需要自己编译,并将其复制到/usr/local/lib目录中。在下面的网址中下载libf2c.zip,然后按照readme文件进行安装。
http://www.netlib.org/f2c/
3. 使用Armadillo例子
1) 基本例子
这个是Armadillo自带的一个例子。利用下面的命令进行编译,我的头文件放在目录/home/xuehen/armadillo-6.100.0/include中,这里只是对密集矩阵进行求解工作,所以只需要liblapack.a libblas.a 库文件即可。-lgfortran是加上gortran的库。
g++ example1.cpp -o example1 -O2 -I/home/xuehen/armadillo-6.100.0/include -DARMA_DONT_USE_WRAPPER -llapack -lblas -lgfortran
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
// Armadillo documentation is available at:
// http://arma.sourceforge.net/docs.html
int
main(int argc, char** argv)
{
cout << "Armadillo version: " << arma_version::as_string() << endl;
mat A(2,3); // directly specify the matrix size (elements are uninitialised)
cout << "A.n_rows: " << A.n_rows << endl; // .n_rows and .n_cols are read only
cout << "A.n_cols: " << A.n_cols << endl;
A(1,2) = 456.0; // directly access an element (indexing starts at 0)
A.print("A:");
A = 5.0; // scalars are treated as a 1x1 matrix
A.print("A:");
A.set_size(4,5); // change the size (data is not preserved)
A.fill(5.0); // set all elements to a particular value
A.print("A:");
// endr indicates "end of row"
A << 0.165300 << 0.454037 << 0.995795 << 0.124098 << 0.047084 << endr
<< 0.688782 << 0.036549 << 0.552848 << 0.937664 << 0.866401 << endr
<< 0.348740 << 0.479388 << 0.506228 << 0.145673 << 0.491547 << endr
<< 0.148678 << 0.682258 << 0.571154 << 0.874724 << 0.444632 << endr
<< 0.245726 << 0.595218 << 0.409327 << 0.367827 << 0.385736 << endr;
A.print("A:");
// determinant
cout << "det(A): " << det(A) << endl;
// inverse
cout << "inv(A): " << endl << inv(A) << endl;
// save matrix as a text file
A.save("A.txt", raw_ascii);
// load from file
mat B;
B.load("A.txt");
// submatrices
cout << "B( span(0,2), span(3,4) ):" << endl << B( span(0,2), span(3,4) ) << endl;
cout << "B.row(0): " << endl << B.row(0) << endl;
cout << "B.col(1): " << endl << B.col(1) << endl;
// transpose
cout << "B.t(): " << endl << B.t() << endl;
// maximum from each column (traverse along rows)
cout << "max(B): " << endl << max(B) << endl;
// maximum from each row (traverse along columns)
cout << "max(B,1): " << endl << max(B,1) << endl;
// maximum value in B
cout << "max(max(B)) = " << max(max(B)) << endl;
// sum of each column (traverse along rows)
cout << "sum(B): " << endl << sum(B) << endl;
// sum of each row (traverse along columns)
cout << "sum(B,1) =" << endl << sum(B,1) << endl;
// sum of all elements
cout << "accu(B): " << accu(B) << endl;
// trace = sum along diagonal
cout << "trace(B): " << trace(B) << endl;
// generate the identity matrix
mat C = eye<mat>(4,4);
// random matrix with values uniformly distributed in the [0,1] interval
mat D = randu<mat>(4,4);
D.print("D:");
// row vectors are treated like a matrix with one row
rowvec r;
r << 0.59119 << 0.77321 << 0.60275 << 0.35887 << 0.51683;
r.print("r:");
// column vectors are treated like a matrix with one column
vec q;
q << 0.14333 << 0.59478 << 0.14481 << 0.58558 << 0.60809;
q.print("q:");
// convert matrix to vector; data in matrices is stored column-by-column
vec v = vectorise(A);
v.print("v:");
// dot or inner product
cout << "as_scalar(r*q): " << as_scalar(r*q) << endl;
// outer product
cout << "q*r: " << endl << q*r << endl;
// multiply-and-accumulate operation (no temporary matrices are created)
cout << "accu(A % B) = " << accu(A % B) << endl;
// example of a compound operation
B += 2.0 * A.t();
B.print("B:");
// imat specifies an integer matrix
imat AA;
imat BB;
AA << 1 << 2 << 3 << endr << 4 << 5 << 6 << endr << 7 << 8 << 9;
BB << 3 << 2 << 1 << endr << 6 << 5 << 4 << endr << 9 << 8 << 7;
// comparison of matrices (element-wise); output of a relational operator is a umat
umat ZZ = (AA >= BB);
ZZ.print("ZZ:");
// cubes ("3D matrices")
cube Q( B.n_rows, B.n_cols, 2 );
Q.slice(0) = B;
Q.slice(1) = 2.0 * B;
Q.print("Q:");
// 2D field of matrices; 3D fields are also supported
field<mat> F(4,3);
for(uword col=0; col < F.n_cols; ++col)
for(uword row=0; row < F.n_rows; ++row)
{
F(row,col) = randu<mat>(2,3); // each element in field<mat> is a matrix
}
F.print("F:");
return 0;
}
2). 求解对称稀疏实矩阵的例子
g++ example.cpp -o example -O2 -I/home/xuehen/armadillo-6.100.0/include /usr/local/lib/libarpack.a /usr/local/lib/liblapack.a /usr/local/lib/libopenblas.a /usr/local/lib/libf2c.a -lgfortran
g++ example.cpp -o example -O2 -I/home/xuehen/armadillo-6.100.0/include -L/usr/local/lib -larpack -llapack -lopenblas -lf2c -lgfortran
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
// Armadillo documentation is available at:
// http://arma.sourceforge.net/docs.html
int
main(int argc, char** argv)
{
// generate sparse matrix
sp_fmat A = sprandu<sp_fmat>(1000, 1000, 0.1);
sp_fmat B = A.t()*A;
fvec eigval;
fmat eigvec;
eigs_sym(eigval, eigvec, B, 1); // find 5 eigenvalues/eigenvectors
}
4.总结
由于尝试对Armadillo打包时出现链接错误,所以这里才用了不安装的模式调用Armadillo。