libicr: Linear algebra and decompositions

Catalogue of decompositions offered by Eigen

	Generic information, not Eigen-specific					Eigen-specific
Decomposition	Requirements on the matrix	Speed	Algorithm reliability and accuracy	Rank-revealing	Allows to compute (besides linear solving)	Linear solver provided by Eigen	Maturity of Eigen's implementation	Optimizations
PartialPivLU	Invertible	Fast	Depends on condition number	-	-	Yes	Excellent	Blocking
FullPivLU	-	Slow	Proven	Yes	-	Yes	Excellent	-
HouseholderQR	-	Fast	Depends on condition number	-	Orthogonalization	Yes	Excellent	Blocking
ColPivHouseholderQR	-	Fast	Good	Yes	Orthogonalization	Yes	Excellent	Soon: blocking
FullPivHouseholderQR	-	Slow	Proven	Yes	Orthogonalization	Yes	Average	-
LLT	Positive definite	Very fast	Depends on condition number	-	-	Yes	Excellent	Blocking Soon: meta unroller
LDLT	Positive or negative semidefinite¹	Very fast	Good	-	-	Yes	Excellent	Soon: blocking
Singular values and eigenvalues decompositions
JacobiSVD (two-sided)	-	Slow (but fast for small matrices)	Excellent-Proven³	Yes	Singular values/vectors, least squares	Yes (and does least squares)	Excellent	R-SVD
SelfAdjointEigenSolver	Self-adjoint	Fast-average²	Good	Yes	Eigenvalues/vectors	-	Good	Soon: specializations for 2x2 and 3x3
ComplexEigenSolver	Square	Slow-very slow²	Depends on condition number	Yes	Eigenvalues/vectors	-	Average	-
EigenSolver	Square and real	Average-slow²	Depends on condition number	Yes	Eigenvalues/vectors	-	Average	-
GeneralizedSelfAdjointEigenSolver	Square	Fast-average²	Depends on condition number	-	Generalized eigenvalues/vectors	-	Good	-
Helper decompositions
RealSchur	Square and real	Average-slow²	Depends on condition number	Yes	-	-	Average	-
ComplexSchur	Square	Slow-very slow²	Depends on condition number	Yes	-	-	Average	-
Tridiagonalization	Self-adjoint	Fast	Good	-	-	-	Good	Soon: blocking
HessenbergDecomposition	Square	Average	Good	-	-	-	Good	Soon: blocking

Notes:

1: There exist two variants of the LDLT algorithm. Eigen's one produces a pure diagonal D matrix, and therefore it cannot handle indefinite matrices, unlike Lapack's one which produces a block diagonal D matrix.
2: Eigenvalues, SVD and Schur decompositions rely on iterative algorithms. Their convergence speed depends on how well the eigenvalues are separated.
3: Our JacobiSVD is two-sided, making for proven and optimal precision for square matrices. For non-square matrices, we have to use a QR preconditioner first. The default choice, ColPivHouseholderQR, is already very reliable, but if you want it to be proven, use FullPivHouseholderQR instead.

Terminology

Selfadjoint

For a real matrix, selfadjoint is a synonym for symmetric. For a complex matrix, selfadjoint is a synonym for hermitian. More generally, a matrix $ A $

is selfadjoint if and only if it is equal to its adjoint $ A^* $

. The adjoint is also called the conjugate transpose.

Positive/negative definite

A selfadjoint matrix $ A $

is positive definite if $ v^* A v > 0 $

for any non zero vector $ v $

. In the same vein, it is negative definite if $ v^* A v < 0 $

for any non zero vector $ v $

Positive/negative semidefinite

A selfadjoint matrix $ A $ is positive semi-definite if $v^* A v \ge 0$ for any non zero vector $ v $ . In the same vein, it is negative semi-definite if $v^* A v \le 0$ for any non zero vector $ v $

Blocking

Means the algorithm can work per block, whence guaranteeing a good scaling of the performance for large matrices.

Meta-unroller

Means the algorithm is automatically and explicitly unrolled for very small fixed size matrices.