phynexis-deps/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
master
eigen/doc/TopicLinearAlgebraDecompositions.dox
T
Rasmus Munk LarsenandRasmus Munk Larsen e857597ca0 Docs: Correct three claims the manual audit left behind 
libeigen/eigen!2857

Closes #307 and #732

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
2026-08-17 20:52:31 -07:00

418 lines
18 KiB
Plaintext
Raw Permalink Blame History

namespace Eigen {
/** \eigenManualPage TopicLinearAlgebraDecompositions Catalogue of dense decompositions
This page presents a catalogue of the dense matrix decompositions offered by Eigen.
For an introduction on linear solvers and decompositions, check this \link TutorialLinearAlgebra page \endlink.
To get an overview of the true relative speed of the different decompositions, check this \link DenseDecompositionBenchmark benchmark \endlink.
\section TopicLinAlgBigTable Catalogue of decompositions offered by Eigen
<table class="manual-vl">
<tr>
<th class="meta"></th>
<th class="meta" colspan="5">Generic information, not Eigen-specific</th>
<th class="meta" colspan="3">Eigen-specific</th>
</tr>
<tr>
<th>Decomposition</th>
<th>Requirements on the matrix</th>
<th>Speed</th>
<th>Algorithm reliability and accuracy</th>
<th>Rank-revealing</th>
<th>Allows to compute (besides linear solving)</th>
<th>Linear solver provided by Eigen</th>
<th>Maturity of Eigen's implementation</th>
<th>Optimizations</th>
</tr>
<tr>
<td>PartialPivLU</td>
<td>Invertible</td>
<td>Fast</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>-</td>
<td>-</td>
<td>Yes</td>
<td>Excellent</td>
<td>Blocking, Implicit MT</td>
</tr>
<tr class="alt">
<td>FullPivLU</td>
<td>-</td>
<td>Slow (no blocking)</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Rank, kernel, image</td>
<td>Yes</td>
<td>Excellent</td>
<td>-</td>
</tr>
<tr>
<td>HouseholderQR</td>
<td>-</td>
<td>Fast</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>-</td>
<td>Orthogonalization, least squares for overdetermined systems</td>
<td>Yes (and does least squares)</td>
<td>Excellent</td>
<td>Blocking</td>
</tr>
<tr class="alt">
<td>ColPivHouseholderQR</td>
<td>-</td>
<td>Fast</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Orthogonalization, least squares for overdetermined systems</td>
<td>Yes (and does least squares)</td>
<td>Excellent</td>
<td><em>-</em></td>
</tr>
<tr>
<td>RandColPivHouseholderQR</td>
<td>-</td>
<td>Fast (large matrices)</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Orthogonalization, least squares for overdetermined systems</td>
<td>Yes (and does least squares)</td>
<td>Good</td>
<td>Blocking, randomized pivot selection</td>
</tr>
<tr class="alt">
<td>FullPivHouseholderQR</td>
<td>-</td>
<td>Slow (no blocking)</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Orthogonalization, least squares for overdetermined systems</td>
<td>Yes (and does least squares)</td>
<td>Average</td>
<td>-</td>
</tr>
<tr>
<td>CompleteOrthogonalDecomposition</td>
<td>-</td>
<td>Fast</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Orthogonalization, minimum-norm least squares, pseudo-inverse</td>
<td>Yes (and does least squares)</td>
<td>Excellent</td>
<td><em>-</em></td>
</tr>
<tr class="alt">
<td>RandCompleteOrthogonalDecomposition</td>
<td>-</td>
<td>Fast (large matrices)</td>
<td>Good<sup><a href="#note4">4</a></sup></td>
<td>Yes</td>
<td>Orthogonalization, minimum-norm least squares, pseudo-inverse</td>
<td>Yes (and does least squares)</td>
<td>Good</td>
<td>Blocking, randomized pivot selection</td>
</tr>
<tr>
<td>LLT</td>
<td>Positive definite</td>
<td>Very fast</td>
<td>Good<sup><a href="#note5">5</a></sup></td>
<td>-</td>
<td>-</td>
<td>Yes</td>
<td>Excellent</td>
<td>Blocking</td>
</tr>
<tr class="alt">
<td>LDLT</td>
<td>Positive or negative semidefinite<sup><a href="#note1">1</a></sup></td>
<td>Very fast</td>
<td>Good</td>
<td>-</td>
<td>-</td>
<td>Yes</td>
<td>Excellent</td>
<td>-</td>
</tr>
<tr>
<td>BunchKaufman</td>
<td>Self-adjoint (possibly indefinite)</td>
<td>Very fast</td>
<td>Good</td>
<td>-</td>
<td>Inertia</td>
<td>Yes</td>
<td>Good</td>
<td>Blocking</td>
</tr>
<tr><th class="inter" colspan="9">\n Singular values and eigenvalues decompositions</th></tr>
<tr>
<td>BDCSVD (divide \& conquer)</td>
<td>-</td>
<td>One of the fastest SVD algorithms</td>
<td>Excellent</td>
<td>Yes</td>
<td>Singular values/vectors, least squares</td>
<td>Yes (and does least squares)</td>
<td>Excellent</td>
<td>Blocked bidiagonalization</td>
</tr>
<tr class="alt">
<td>JacobiSVD (two-sided)</td>
<td>-</td>
<td>Slow (but fast for small matrices)</td>
<td>High relative accuracy<sup><a href="#note3">3</a></sup></td>
<td>Yes</td>
<td>Singular values/vectors, least squares</td>
<td>Yes (and does least squares)</td>
<td>Excellent</td>
<td>R-SVD</td>
</tr>
<tr>
<td>SelfAdjointEigenSolver</td>
<td>Self-adjoint</td>
<td>Fast-average<sup><a href="#note2">2</a></sup></td>
<td>Good</td>
<td>-</td>
<td>Eigenvalues/vectors</td>
<td>-</td>
<td>Excellent</td>
<td><em>Closed forms for 2x2 and 3x3</em></td>
</tr>
<tr class="alt">
<td>TridiagonalEigenSolver</td>
<td>Symmetric tridiagonal</td>
<td>Fast</td>
<td>Good</td>
<td>-</td>
<td>Eigenvalues/vectors, subsets of the spectrum</td>
<td>-</td>
<td>Good</td>
<td>Vectorization, Explicit MT</td>
</tr>
<tr>
<td>ComplexEigenSolver</td>
<td>Square</td>
<td>Slow-very slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Eigenvalues/vectors</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr class="alt">
<td>EigenSolver</td>
<td>Square and real</td>
<td>Average-slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Eigenvalues/vectors</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr>
<td>GeneralizedSelfAdjointEigenSolver</td>
<td>Square</td>
<td>Fast-average<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Generalized eigenvalues/vectors</td>
<td>-</td>
<td>Good</td>
<td>-</td>
</tr>
<tr class="alt">
<td>GeneralizedEigenSolver</td>
<td>Square and real</td>
<td>Slow-very slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Generalized eigenvalues/vectors</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr><th class="inter" colspan="9">\n Helper decompositions</th></tr>
<tr>
<td>RealSchur</td>
<td>Square and real</td>
<td>Average-slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr class="alt">
<td>ComplexSchur</td>
<td>Square</td>
<td>Slow-very slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr>
<td>RealQZ</td>
<td>Square and real (pair of matrices)</td>
<td>Average-slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Generalized eigenvalues of a pencil</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr class="alt">
<td>ComplexQZ</td>
<td>Square (pair of matrices)</td>
<td>Slow-very slow<sup><a href="#note2">2</a></sup></td>
<td>Depends on condition number</td>
<td>-</td>
<td>Generalized eigenvalues of a pencil</td>
<td>-</td>
<td>Average</td>
<td>-</td>
</tr>
<tr>
<td>Tridiagonalization</td>
<td>Self-adjoint</td>
<td>Fast</td>
<td>Good</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>Good</td>
<td>-</td>
</tr>
<tr class="alt">
<td>HessenbergDecomposition</td>
<td>Square</td>
<td>Average</td>
<td>Good</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>Good</td>
<td>-</td>
</tr>
</table>
\b Notes:
<ul>
<li><a name="note1">\b 1: </a>There exist two variants of the LDLT algorithm. %Eigen's LDLT produces a pure diagonal D matrix and therefore cannot handle indefinite matrices; it uses symmetric pivoting, which keeps the factorization stable when pivots are tiny or zero. The other variant, which produces the block-diagonal D needed for indefinite self-adjoint matrices, is provided by the BunchKaufman class.</li>
<li><a name="note2">\b 2: </a>Eigenvalues, SVD and Schur decompositions rely on iterative algorithms. Their convergence speed depends on how well the eigenvalues are separated.</li>
<li><a name="note3">\b 3: </a>Our JacobiSVD is two-sided, which buys high relative accuracy on square matrices: where a bidiagonalizing algorithm only computes each singular value accurately relative to the largest one, Jacobi computes it accurately relative to itself, so that tiny singular values are still meaningful. That property is conditional rather than absolute -- the relative error is governed by the condition number of the matrix obtained by scaling its columns to unit norm, so a matrix that stays ill-conditioned under any such scaling loses it (J. Demmel and K. Veselic, "Jacobi's method is more accurate than QR", SIAM J. Matrix Anal. Appl. 13 (1992), pp. 1204-1245). For non-square matrices, we have to use a QR preconditioner first. Householder QR is backward stable with or without pivoting, so this choice is about rank revelation rather than about the accuracy of the factorization itself: the default \c ColPivHouseholderQRPreconditioner reveals rank reliably in practice but without a guarantee, \c FullPivHouseholderQRPreconditioner inspects the whole trailing submatrix at a substantially higher cost, and \c HouseholderQRPreconditioner and \c NoQRPreconditioner do not pivot at all. See the guidance below on complete pivoting.</li>
<li><a name="note4">\b 4: </a>These factorizations are all backward stable, but pivoting buys something different in each. For LU the backward error bound is proportional to the growth factor: partial pivoting bounds it only by \f$ 2^{n-1} \f$, whereas complete pivoting admits a far smaller bound that grows only slightly faster than polynomially in \c n. That partial-pivoting worst case is not observed on real matrices, however, and partial pivoting is more than sufficient in practice. Householder QR is backward stable irrespective of pivoting, with a bound that involves no growth factor at all; there, pivoting serves rank revelation rather than stability, and what separates column from complete pivoting is how reliably the trailing diagonal of R exposes a rank deficiency. Since complete pivoting also rules out the cache-friendly blocked algorithms, LAPACK and %Eigen both treat partial and column pivoting as the defaults for LU and QR; see the guidance below. This is a judgement about those two factorizations rather than about complete pivoting in general: LAPACK does provide a complete-pivoting Cholesky (\c pstrf) for positive semidefinite matrices, where the pivoting also determines the rank.</li>
<li><a name="note5">\b 5: </a>Cholesky needs no pivoting for stability: on a positive definite matrix the growth factor is bounded by 1, so the unpivoted factorization is backward stable as it stands. What conditioning governs is instead whether it completes: on a matrix that is positive definite but nearly singular, rounding can drive a computed pivot to zero or below, and %Eigen reports that through info() rather than returning a meaningless factor. Use LDLT, or a rank-revealing decomposition, for such matrices.</li>
</ul>
\section TopicLinAlgPracticalGuidance Practical guidance
The following recommendations apply to the most common use cases:
\li <b>Symmetric positive definite systems:</b> Use \b LLT. It is the fastest solver and has excellent
numerical properties for this class of problems. For semidefinite or nearly singular symmetric systems,
use \b LDLT.
\li <b>Symmetric indefinite systems:</b> Use \b BunchKaufman. It extends the Cholesky family to
indefinite self-adjoint matrices with a backward-stable pivoting strategy at LDLT-like speed,
and also reveals the inertia of the matrix.
\li <b>General invertible systems:</b> Use \b PartialPivLU. It uses cache-friendly blocking and implicit
multi-threading, making it the fastest general-purpose solver. Partial pivoting is sufficient for
virtually all practical problems.
\li <b>Least squares (over- or under-determined systems):</b> Use \b CompleteOrthogonalDecomposition as
the default. Like the SVD, it robustly computes the minimum-norm solution for rank-deficient and
under-determined problems, but at QR-like speed. Use \b BDCSVD when you also need singular values
or vectors, not just the least squares solution.
\li <b>Large least squares problems:</b> Use \b RandCompleteOrthogonalDecomposition. It produces the
same kind of factorization as \b CompleteOrthogonalDecomposition, but chooses pivots from a small random
sketch so that nearly all of the work happens in level-3 BLAS kernels rather than in the
memory-bound column scan of classical column pivoting. The pivot-quality difference is empirically
minor. On small matrices the sketch overhead dominates, so prefer the non-randomized class there.
\li <b>Pseudo-inverse:</b> \b CompleteOrthogonalDecomposition and \b RandCompleteOrthogonalDecomposition
are the only decompositions that provide pseudoInverse() directly. Prefer solve() over multiplying
by an explicitly formed pseudo-inverse.
\li <b>Full-rank least squares (overdetermined systems):</b> When the matrix is known to be full rank,
\b HouseholderQR is the fastest option. For very tall and skinny well-conditioned matrices,
solving via the normal equations with \b LLT can be faster still.
\li <b>FullPivLU and FullPivHouseholderQR</b> use complete pivoting, which rules out the
cache-friendly blocked algorithms and so makes them substantially slower than their partial- and
column-pivoting counterparts, increasingly so as the matrix grows. What complete pivoting buys --
a smaller growth factor for LU, a more reliable rank decision for QR -- rarely repays that cost on
its own, so PartialPivLU and ColPivHouseholderQR are the right default. LAPACK reflects the same
trade-off for LU and QR: it exposes no complete-pivoting driver for either, standardizing instead
on partial-pivoting LU (\c getrf) and column-pivoted QR (\c geqp3, \c geqp3rk). Complete-pivoting
LU survives there only as the auxiliary routine \c getc2, called by the generalized Sylvester
solver, and complete pivoting has no QR counterpart at all.
\li <b>When complete pivoting is worth it:</b> rank determination. Column pivoting is reliable in
practice but not guaranteed (Kahan's matrix is the standard example of a nearly rank-deficient
matrix whose column-pivoted QR exposes no small trailing diagonal entry), whereas complete
pivoting inspects the whole trailing submatrix at each step. If your application depends on
getting the rank right, FullPivLU or FullPivHouseholderQR is a defensible choice, and the
argument is strongest for small matrices, where the lost blocking costs comparatively little.
When rank is the whole point and the cost is acceptable, BDCSVD remains the most reliable answer,
since the singular values let you choose the threshold yourself.
\section TopicLinAlgTerminology Terminology
<dl>
<dt><b>Selfadjoint</b></dt>
<dd>For a real matrix, selfadjoint is a synonym for symmetric. For a complex matrix, selfadjoint is a synonym for \em hermitian.
More generally, a matrix \f$ A \f$ is selfadjoint if and only if it is equal to its adjoint \f$ A^* \f$. The adjoint is also called the \em conjugate \em transpose. </dd>
<dt><b>Positive/negative definite</b></dt>
<dd>A selfadjoint matrix \f$ A \f$ is positive definite if \f$ v^* A v > 0 \f$ for any non zero vector \f$ v \f$.
In the same vein, it is negative definite if \f$ v^* A v < 0 \f$ for any non zero vector \f$ v \f$ </dd>
<dt><b>Positive/negative semidefinite</b></dt>
<dd>A selfadjoint matrix \f$ A \f$ is positive semi-definite if \f$ v^* A v \ge 0 \f$ for any non zero vector \f$ v \f$.
In the same vein, it is negative semi-definite if \f$ v^* A v \le 0 \f$ for any non zero vector \f$ v \f$ </dd>
<dt><b>Blocking</b></dt>
<dd>Means the algorithm can work per block, whence guaranteeing a good scaling of the performance for large matrices.</dd>
<dt><b>Implicit Multi Threading (MT)</b></dt>
<dd>Means the algorithm can take advantage of multicore processors via OpenMP. "Implicit" means the algorithm itself is not parallelized, but that it relies on parallelized matrix-matrix product routines.</dd>
<dt><b>Explicit Multi Threading (MT)</b></dt>
<dd>Means the algorithm is explicitly parallelized to take advantage of multicore processors via OpenMP.</dd>
<dt><b>Meta-unroller</b></dt>
<dd>Means the algorithm is automatically and explicitly unrolled for very small fixed size matrices.</dd>
</dl>
*/
}
Reference in New Issue View Git Blame Copy Permalink