Back to my personal website

Eigen

Examples

#include <Eigen/Core>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

// declarations
Matrix<float, 2, 3> matrix_23; // type, row, col
Vector3d v_3d;      // Matrix<float, 3, 1>
Matrix3d matrix_33; // Matrix<double, 3, 3>
MatrixXd matrix_x;  // Matrix<double, Dynamic, Dynamic>

// input & output
matrix_33 = Matrix3d::Zero();
matrix_23 << 1, 2, 3, 4, 5, 6;
cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl;
cout << "matrix(1, 2): " << matrix_23(0, 1) << endl;
v_3d << 3, 2, 1;

// requires explicit type conversion (float to double)
Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
cout << "[1,2,3;4,5,6]*[3,2,1]=" << result.transpose() << endl;

// matrix operations
matrix_33 = Matrix3d::Random();
cout << "random matrix: \n" << matrix_33 << endl;
cout << "transpose: \n" << matrix_33.transpose() << endl;
cout << "sum: " << matrix_33.sum() << endl;
cout << "trace: " << matrix_33.trace() << endl;
cout << "times 10: \n" << 10 * matrix_33 << endl;
cout << "inverse: \n" << matrix_33.inverse() << endl;
cout << "det: " << matrix_33.determinant() << endl;

// eigen value
SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);
cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;

// solve equations
Matrix<double, N, N> matrix_NN = MatrixXd::Random(N, N);
matrix_NN = matrix_NN * matrix_NN.transpose();  // for semi-definite
Matrix<double, N, 1> v_Nd = MatrixXd::Random(N, 1);
x = matrix_NN.colPivHouseholderQr().solve(v_Nd);  // QR decomposition
x = matrix_NN.ldlt().solve(v_Nd);  // cholesky decomposition

// timing
clock_t tic = clock();
cout << "time is " << 1000 * (clock() - tic) / (double) CLOCKS_PER_SEC << "ms" << endl;

Linear Solver

We discuss some properties of the matrix ATAA^TA, constructed from a m-by-n matrix AA where m is significantly greater than n. This is normally the case in least-square problems.

  • The matrix ATAA^TA is always symmetric by definition.

  • The matrix ATAA^TA is positive semi-definite if it is invertible, since we havexTATAx=(Ax)TAx=∥Ax∥2≥0x^TA^TAx=(Ax)^TAx=∥Ax∥^2≥0.

  • The matrix ATAA^TA is almost always invertible, because m is significantly greater than n, and hence the constructed matrix is high likely to be full-rank.

  • The dense linear solver ldlt() in Eigen is the fastest for positive or negative semi-definite matrices. Remember to use double precision to avoid numerical instability issues.

References: Quora: Is ATA always positive definite?; StackExchange: Is a matrix multiplied with its transpose something special?; StackExchange: Proof for why a matrix multiplied by its transpose is positive semidefinite; Eigen Documentation: Basic linear solving;

Numerical Stability

The following logs reveal one of the potential issues of using A.ldlt().solve(b);

  • The solver was used for Gauss-Newton optimization in an Iterative Closest Point algorithm. It was tested under the termination condition of max_iteration only.

  • The result suggests that we should check if A and b are healthy before sending to A.ldlt().solve(b);, and it might be better to add epsilon as another termination condition.

  • Note that float type supports a value range of 1.2E-38 to 3.4E+38, and 6-7 decimal digits, while double type supports a value range of 2.3E-308 to 1.7E+308, and 15 decimal digits.

  • Updates: This issue was resolved by switching to double precision. (Tested multiple times towards 100 iterations. The produced result is stable.)

PreviousOther LibrariesNextCeres

Last updated

Was this helpful?