【高性能】Eigen VS Matlab

Eigen是经典的C++开源模板矩阵库，很多大型库都包含对Eigen的支持，对矩阵计算有良好的优化。Matlab是科学计算领域的佼佼者，对矩阵计算有得天独厚的优势。将Matlab算法改编为C++程序，使用Eigen库会方便很多。

以下内容转载自Eigen官网，我进行了翻译：

// Basic usage

// Eigen          // Matlab           // 注释

x.size()          // length(x)        // 尺寸

C.rows()          // size(C, 1)       // 列数

C.cols()          // size(C, 2)       // 行数

x(i)              // x(i+1)           // Matlab is 1-based

C(i,j)            // C(i+1,j+1)       //最头疼的地方，matlab和C++序号差1

 

                 

A << 1, 2, 3,     //A= [1,2,3;4,5,6;7,8,9]

    4, 5, 6,    

    7, 8, 9;    

B << A, A, A;     // B = [A,A,A];

A.fill(10);       // A = 10 * Ones(3,3);

A.setRandom();    //A = rand(3,3)

A.setIdentity();  // A = eye(3,3);

 

// 矩阵子区域

// 模版形式是很快的。注意matlab从1开始

// Eigen                           // Matlab

x.head(n)                          // x(1:n)

x.head<n>()                        // x(1:n)

x.tail(n)                           // N = rows(x); x(N - n: N)

x.tail<n>()                         // N = rows(x); x(N - n: N)

x.segment(i, n)                      // x(i+1 : i+n)

x.segment<n>(i)                    // x(i+1 : i+n)

P.block(i, j, rows, cols)              //P(i+1 : i+rows, j+1 : j+cols)

P.block<rows, cols>(i, j)             // P(i+1 : i+rows, j+1 : j+cols)

P.topLeftCorner(rows, cols)          // P(1:rows, 1:cols)

P.topRightCorner(rows, cols)         //[m ,n]=size(P); P(1:rows, n-cols+1:n)

P.bottomLeftCorner(rows, cols)       //[m ,n]=size(P); P(m-rows+1:m, 1:cols)

P.bottomRightCorner(rows, cols)     //[m ,n]=size(P); P(m-rows+1:m, n-cols+1:n)

P.topLeftCorner<rows,cols>()       // P(1:rows, 1:cols)

P.topRightCorner<rows,cols>()      // [m ,n]=size(P); P(1:rows, n-cols+1:n)

P.bottomLeftCorner<rows,cols>()    // [m ,n]=size(P); P(m-rows+1:m, 1:cols)

P.bottomRightCorner<rows,cols>()   // [m ,n]=size(P); P(m-rows+1:m, n-cols+1:n)

 

// swap是高度优化的，很快

// Eigen                           // Matlab

R.row(i) = P.col(j);               // R(i, :) = P(:, i)

R.col(j1).swap(mat1.col(j2));       //R(:, [j1, j2]) = R(:, [j2, j1])

 

// 矩阵变换

// Eigen                           // Matlab

R.adjoint()                        // R'

R.transpose()                      // R.' or conj(R')

R.diagonal()                       // diag(R)

x.asDiagonal()                     // diag(x)

 

// matlab 没有 *= 这种操作，其他都很像

// Matrix-vector.  Matrix-matrix.   Matrix-scalar.

y  =M*x;          R  = P*Q;       R  = P*s;

a  =b*M;          R  = P - Q;     R  = s*P;

a *= M;            R = P + Q;      R  = P/s;

                   R *= Q;          R = s*P;

                   R += Q;          R *= s;

                   R -= Q;          R /= s;

 

 // 把矩阵变成数列

// Eigen                        // Matlab

R = P.cwiseProduct(Q);             // R = P .* Q

R = P.array() * s.array();                     // R = P .* s

R = P.cwiseQuotient(Q);            // R = P ./ Q

R = P.array() / Q.array();                    // R = P ./ Q

R = P.array() + s.array();                     // R = P + s

R = P.array() - s.array();                     // R = P - s

R.array() += s;                    // R = R + s

R.array() -= s;                    // R = R - s

R.array() < Q.array();               // R < Q

R.array() <= Q.array();              // R <= Q

R.cwiseInverse();           //1 ./ P

R.array().inverse();          //1 ./ P

R.array().sin()             //sin(P)

R.array().cos()             //cos(P)

R.array().pow(s)          //P .^ s

R.array().square()           //P .^ 2

R.array().cube()            //P .^ 3

R.cwiseSqrt()             // sqrt(P)

R.array().sqrt()             //sqrt(P)

R.array().exp()           //exp(P)

R.array().log()             //log(P)

R.cwiseMax(P)            // max(R, P)

R.array().max(P.array())      //max(R, P)

R.cwiseMin(P)            // min(R, P)

R.array().min(P.array())      //min(R, P)

R.cwiseAbs()             // abs(P)

R.array().abs()             //abs(P)

R.cwiseAbs2()            // abs(P.^2)

R.array().abs2()            //abs(P.^2)

(R.array() < s).select(P,Q);  // (R< s ? P : Q)

 

// 矩阵数据单元的计算.

int r, c;

// Eigen                  // Matlab

R.minCoeff()              // min(R(:))

R.maxCoeff()              // max(R(:))

s = R.minCoeff(&r, &c)    //[aa, bb] = min(R); [cc, dd] = min(aa);

                        // r = bb(dd); c = dd;s = cc

s = R.maxCoeff(&r, &c)           //[aa, bb] = max(R); [cc, dd] = max(aa);

                        // row = bb(dd); col =dd; s = cc

R.sum()                  // sum(R(:))

R.colwise.sum()           // sum(R)

R.rowwise.sum()          // sum(R, 2) or sum(R')'

R.prod()                 // prod(R(:))

R.colwise.prod()          // prod(R)

R.rowwise.prod()          // prod(R, 2) or prod(R')'

R.trace()                 // trace(R)

R.all()                   // all(R(:))

R.colwise().all()            //all(R)

R.rowwise().all()           //all(R, 2)

R.any()                  // any(R(:))

R.colwise().any()           //any(R)

R.rowwise().any()         // any(R, 2)

 

// 乘法，归一化

// Eigen                  // Matlab

x.norm()                  // norm(x).    .

x.squaredNorm()           // dot(x, x)  

x.dot(y)                  // dot(x, y)

x.cross(y)                // cross(x, y)  需要加#include <Eigen/Geometry>

 

// 把现成的数据复制到Eigenvalues矩阵

float array[3];

Map<Vector3f>(array, 3).fill(10);

int data[4] = 1, 2, 3, 4;

Matrix2i mat2x2(data);

MatrixXi mat2x2 =Map<Matrix2i>(data);

MatrixXi mat2x2 = Map<MatrixXi>(data,2, 2);

 

// 求解Ax = b. x存放结果. Matlab: x = A \ b.

bool solved;

solved = A.ldlt().solve(b, &x));  //    #include <Eigen/Cholesky>

solved = A.llt() .solve(b, &x));  //      #include <Eigen/Cholesky>

solved = A.lu() .solve(b, &x));  // 稳定快速. #include <Eigen/LU>

solved = A.qr() .solve(b, &x));  //    #include <Eigen/QR>

solved = A.svd() .solve(b, &x));  // 稳定，慢 #include <Eigen/SVD>

//存放的函数

// .ldlt() -> .matrixL() and .matrixD()

// .llt() -> .matrixL()

// .lu() -> .matrixL() and .matrixU()

// .qr() -> .matrixQ() and .matrixR()

// .svd() -> .matrixU(),.singularValues(), 和.matrixV()

 

// 其他

// Eigen                          // Matlab

A.eigenvalues();                  // eig(A);

EigenSolver<Matrix3d> eig(A);     // [vec, val] = eig(A)

eig.eigenvalues();                // diag(val)