Cholesky: Pivot LDLT on the Schur complement's diagonal, and block it libeigen/eigen!3423 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/CHANGELOG.md b/CHANGELOG.md index 932fe6f..addeddf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md
@@ -22,8 +22,10 @@ - `JacobiSVD` with `FullPivHouseholderQRPreconditioner` accepted `ComputeThinV` and returned a non-orthogonal `V` for wide inputs; it now fails to compile, or asserts when requested at runtime [!3212] - `HouseholderQR` (and the decompositions and `SplineFitting` built on it) and `JacobiSVD` did not compile with `AutoDiffScalar` [#3192, !3265] - Constructing or assigning an `Affine` or `Isometry` `Transform` from a (Dim+1)x(Dim+1) matrix copied its last row verbatim; the last row is now set to `[0 ... 0 1]`, as for Dim x Dim and Dim x (Dim+1) inputs [#2715, !3412] +- `LDLT` pivoted on the diagonal of its input rather than of the remaining Schur complement, so zero pivots of a semidefinite matrix could come before non-zero ones and report `NumericalIssue`, as for `[1 1 0; 1 1 0; 0 0 1]`; a NaN or infinite pivot reported `Success`. `LDLT` now also factors dynamic-size matrices larger than 96 in blocks, with a panel width that grows as sqrt(3 n): 1.7x faster at n = 128, 3.1x at n = 1024 and 6.2x at n = 3000 [!3423] Breaking changes: +- `LDLT` picks a different pivot order for most inputs, which changes `transpositionsP()` and the order of `vectorD()`. A NaN or infinite pivot reports `NumericalIssue`. A preallocated `compute()` of a dynamic-size matrix larger than 96 can allocate the packing buffers of its matrix products, as those of `LLT`, `PartialPivLU` and `BunchKaufman` already do [!3423] - The default `epsfcn` for central mode differences was changed to `cbrt(NumTraits<Scalar>::epsilon())` (used to be `sqrt(NumTraits<Scalar>::epsilon())`). The new default should reduce numerical error in most cases. - Removed `Eigen/QtAlignedMalloc`, which replaced Qt 4's `qMalloc`, `qFree` and `qRealloc` with aligned versions on platforms whose `malloc` is not already sufficiently aligned. Remove the include. [!3413]
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 5594e61..12a377d 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h
@@ -33,6 +33,10 @@ // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; + +// Panel width of the blocked factorization (defined below); declared here for the size constructor. +template <typename MatrixType> +inline Index ldlt_blocksize(Index size); } // namespace internal /** \ingroup Cholesky_Module @@ -49,17 +53,21 @@ * matrix \f$ A \f$ such that \f$ A = P^TLDL^*P \f$, where P is a permutation matrix, L * is lower triangular with a unit diagonal and D is a diagonal matrix. * - * The decomposition uses pivoting to ensure stability, so that D will have - * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root - * on D also stabilizes the computation. + * Each step pivots on the diagonal entry of largest magnitude of the remaining Schur complement + * (complete diagonal pivoting, as in LAPACK's xPSTRF), so that for a semidefinite matrix of rank r + * the n - r entries of D that vanish in exact arithmetic come last, at the level of rounding errors + * and of either sign. Avoiding the square root on D also stabilizes the computation. * - * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky - * decomposition to determine whether a system of equations has a solution. + * LDLT does not determine a numerical rank, though. Also, do not use a Cholesky decomposition to + * determine whether a system of equations has a solution. * * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. * - * D is purely diagonal, so this class cannot factor an indefinite matrix. For a self-adjoint matrix - * that is indefinite, use BunchKaufman, which produces a block-diagonal D. + * D is purely diagonal, so this class is not backward stable on indefinite matrices: it can report + * \c Success for a factorization with a large backward error. If isPositive() and isNegative() both + * return false, D has pivots of both signs: unless that comes from the rounding-level pivots of a + * singular semidefinite matrix, the matrix is indefinite, and BunchKaufman, which produces a + * block-diagonal D, should be used instead. * * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT, class BunchKaufman */ @@ -100,6 +108,10 @@ * * Like the default constructor but with preallocation of the internal data * according to the specified problem \a size. + * + * Above 96 rows, compute() factors a dynamic-size matrix in blocks, and the matrix products of the blocked + * factorization can allocate their packing buffers, as those of LLT and PartialPivLU do. Fixed-size and bounded + * matrices are factored unblocked and allocate nothing. * \sa LDLT() */ explicit LDLT(Index size) @@ -107,6 +119,8 @@ m_l1_norm(0), m_transpositions(size), m_temporary(size), + // The shape compute() gives the panel workspace, so that it does not reallocate it for this size. + m_workspace(size, internal::ldlt_blocksize<MatrixType>(size)), m_sign(internal::ZeroSign), m_isInitialized(false), m_info(InvalidInput) {} @@ -311,7 +325,9 @@ /** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was successful, - * \c NumericalIssue if the factorization failed because of a zero pivot. + * \c NumericalIssue if the factorization failed: a pivot is not finite, or is zero above a + * non-zero column. In exact arithmetic only an indefinite matrix has the latter, but rounding + * errors can give it to a singular semidefinite matrix too. */ ComputationInfo info() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); @@ -339,6 +355,8 @@ RealScalar m_l1_norm; TranspositionType m_transpositions; TmpMatrixType m_temporary; + // L21 D for the trailing update of the blocked factorization; it stays empty below the blocking threshold. + Matrix<Scalar, Dynamic, Dynamic> m_workspace; internal::SignMatrix m_sign; bool m_isInitialized; ComputationInfo m_info; @@ -346,21 +364,45 @@ namespace internal { +/** \internal Panel width of the blocked factorization of a \a size x \a size matrix, or 0 where a single + * unblocked panel covers the whole matrix: up to three panel widths, below which the level-3 update does not pay + * off, and for fixed-size or bounded matrices, which must not allocate the panel workspace. + * + * The panel's level-2 updates run at a rate R, and the trailing update at P / (1 + kappa / nb), kappa being the + * fixed cost of a pass in steps of its depth. Then T(nb) = n^3 / (3 P) (1 + kappa / nb) + n^2 nb / (2 R) is least at + * nb = sqrt(c n), c = 2 kappa R / (3 P). Measured c ranges from about 0.7 to 4 with AVX2, AVX-512, NEON and SME + * products, and grows once the trailing matrix outgrows the last-level cache. c = 3, with nb rounded to a multiple + * of 16 and clamped to [32, 128], the range the measured optimum spans, loses about 2% to the best width on average. + * Can be overridden (mainly for testing the panel logic at small sizes) by defining EIGEN_LDLT_BLOCKSIZE. */ +template <typename MatrixType> +inline Index ldlt_blocksize(Index size) { + if (int(MatrixType::MaxRowsAtCompileTime) != Dynamic || int(MatrixType::MaxColsAtCompileTime) != Dynamic) return 0; +#ifdef EIGEN_LDLT_BLOCKSIZE + const Index nb = Index(EIGEN_LDLT_BLOCKSIZE); +#else + constexpr Index kMinBlockSize = 32; + // Three of the narrowest panels are not blocked either, so small matrices skip the square root. + if (size <= 3 * kMinBlockSize) return 0; + const Index unclamped = 16 * Index(numext::round(numext::sqrt(3.0 * double(size)) / 16)); + const Index nb = numext::mini(numext::maxi(unclamped, kMinBlockSize), Index(128)); +#endif + return nb > 0 && size > 3 * nb ? nb : 0; +} + template <int UpLo> struct ldlt_inplace; template <> struct ldlt_inplace<Lower> { + // Settles the cases that leave nothing to factor, and returns whether it did: an empty or 1x1 matrix, and a + // diagonal that is entirely zero, which offers no pivot and leaves the matrix either zero or indefinite. Otherwise + // it copies the diagonal into temp, where factor_panel() expects it. template <typename MatrixType, typename TranspositionType, typename Workspace> - static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { - using std::abs; + static bool factor_trivial(const MatrixType& mat, TranspositionType& transpositions, Workspace& temp, + SignMatrix& sign, bool& ret) { using Scalar = typename MatrixType::Scalar; using RealScalar = typename MatrixType::RealScalar; - using IndexType = typename TranspositionType::StorageIndex; - EIGEN_INTERNAL_ASSERT_SQUARE(mat); const Index size = mat.rows(); - bool found_zero_pivot = false; - bool ret = true; if (size <= 1) { transpositions.setIdentity(); @@ -372,62 +414,81 @@ sign = NegativeSemiDef; else sign = ZeroSign; + ret = size == 0 || (numext::isfinite)(numext::real(mat.coeff(0, 0))); return true; } - for (Index k = 0; k < size; ++k) { + temp = mat.diagonal().real().template cast<Scalar>(); + if ((temp.array() == Scalar(0)).all()) { + transpositions.setIdentity(); + sign = ZeroSign; + ret = true; + for (Index j = 0; ret && j < size; ++j) ret = (mat.col(j).tail(size - j - 1).array() == Scalar(0)).all(); + return true; + } + return false; + } + + // Factors columns [k0, k0 + kb) left-looking, column k pivoting on max |diag(S_k)|, S_k = A22 - L21 D L21^* + // (complete diagonal pivoting, as in LAPACK's xPSTF2). Invariant: mat(k0:n, k0:n) holds S_k0 off its diagonal and + // temp(k:n) holds diag(S_k). The caller interchanges the rows of the columns before k0. + template <typename MatrixType, typename TranspositionType, typename Workspace> + static bool factor_panel(MatrixType& mat, Index k0, Index kb, TranspositionType& transpositions, Workspace& temp, + SignMatrix& sign) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using IndexType = typename TranspositionType::StorageIndex; + const Index size = mat.rows(); + bool ret = true; + + for (Index k = k0; k < k0 + kb; ++k) { // Find largest diagonal element Index index_of_biggest_in_corner; - mat.diagonal().tail(size - k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); + temp.tail(size - k).real().cwiseAbs().maxCoeff(&index_of_biggest_in_corner); index_of_biggest_in_corner += k; transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); if (k != index_of_biggest_in_corner) { - selfadjoint_swap_lower(mat, k, index_of_biggest_in_corner); + auto trailing = mat.bottomRightCorner(size - k0, size - k0); + selfadjoint_swap_lower(trailing, k - k0, index_of_biggest_in_corner - k0); + numext::swap(temp.coeffRef(k), temp.coeffRef(index_of_biggest_in_corner)); } // partition the matrix: // A00 | - | - // lu = A10 | A11 | - // A20 | A21 | A22 + // where A10 and A20 span the panel columns [k0, k); those before k0 are already applied to A21. The + // product's operand D A10^* takes temp(k0:k), whose diagonal entries of S have all been used as pivots. Index rs = size - k - 1; - Block<MatrixType, Dynamic, 1> A21(mat, k + 1, k, rs, 1); - Block<MatrixType, 1, Dynamic> A10(mat, k, 0, 1, k); - Block<MatrixType, Dynamic, Dynamic> A20(mat, k + 1, 0, rs, k); + Index done = k - k0; + auto A21 = mat.block(k + 1, k, rs, fix<1>); - if (k > 0) { - temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); - mat.coeffRef(k, k) -= (A10 * temp.head(k)).value(); - if (rs > 0) A21.noalias() -= A20 * temp.head(k); + if (done > 0 && rs > 0) { + auto A10 = mat.block(k, k0, fix<1>, done); + auto A20 = mat.block(k + 1, k0, rs, done); + temp.segment(k0, done) = mat.diagonal().real().segment(k0, done).asDiagonal() * A10.adjoint(); + A21.noalias() -= A20 * temp.segment(k0, done); } // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot // was smaller than the cutoff value. However, since LDLT is not rank-revealing // we should only make sure that we do not introduce INF or NaN values. // Remark that LAPACK also uses 0 as the cutoff value. - RealScalar realAkk = numext::real(mat.coeffRef(k, k)); - bool pivot_is_valid = (abs(realAkk) > RealScalar(0)); + RealScalar realAkk = numext::real(temp.coeff(k)); + mat.coeffRef(k, k) = Scalar(realAkk); + bool pivot_is_valid = (numext::abs(realAkk) > RealScalar(0)); + if (!(numext::isfinite)(realAkk)) ret = false; - if (k == 0 && !pivot_is_valid) { - // The entire diagonal is zero, there is nothing more to do - // except filling the transpositions, and checking whether the matrix is zero. - sign = ZeroSign; - for (Index j = 0; j < size; ++j) { - transpositions.coeffRef(j) = IndexType(j); - ret = ret && (mat.col(j).tail(size - j - 1).array() == Scalar(0)).all(); - } - return ret; - } - - if ((rs > 0) && pivot_is_valid) + // A zero pivot is the largest remaining, so the rest of D is zero too, and its column must be zero as well. + if ((rs > 0) && pivot_is_valid) { A21 /= realAkk; - else if (rs > 0) + // diag(S_{k+1}) = diag(S_k) - conj(L21) .* (D_k L21), where the next pivot is searched for. Squaring L21 + // first would underflow for |L21| < sqrt(min) on a graded matrix, where D_k |L21|^2 is representable. + temp.tail(rs).real() -= A21.conjugate().cwiseProduct(realAkk * A21).real(); + } else if (rs > 0) { ret = ret && (A21.array() == Scalar(0)).all(); - - if (found_zero_pivot && pivot_is_valid) - ret = false; // factorization failed - else if (!pivot_is_valid) - found_zero_pivot = true; + } if (sign == PositiveSemiDef) { if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite; @@ -444,6 +505,59 @@ return ret; } + template <typename MatrixType, typename TranspositionType, typename Workspace> + static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { + EIGEN_INTERNAL_ASSERT_SQUARE(mat); + bool ret = true; + if (factor_trivial(mat, transpositions, temp, sign, ret)) return ret; + return factor_panel(mat, 0, mat.rows(), transpositions, temp, sign); + } + + // Blocked (level-3) variant, with the panel/trailing-update structure of LAPACK's xPSTRF: each panel is factored + // by factor_panel(), then a single product applies it to the trailing Schur complement, A22 -= L21 D L21^*. + template <typename MatrixType, typename TranspositionType, typename Workspace, typename BlockWorkspace> + static bool blocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, BlockWorkspace& workspace, + SignMatrix& sign) { + EIGEN_INTERNAL_ASSERT_SQUARE(mat); + const Index size = mat.rows(); + const Index nb = ldlt_blocksize<MatrixType>(size); + if (nb == 0) return unblocked(mat, transpositions, temp, sign); + bool ret = true; + if (factor_trivial(mat, transpositions, temp, sign, ret)) return ret; + + workspace.resize(size, nb); + for (Index k0 = 0; k0 < size; k0 += nb) { + const Index kb = numext::mini(nb, size - k0); + const Index k1 = k0 + kb; + if (!factor_panel(mat, k0, kb, transpositions, temp, sign)) ret = false; + + // The panel's interchanges, applied to the rows of the columns before it. Column-major storage takes them a + // column at a time, where a row swap would stride through memory. + EIGEN_IF_CONSTEXPR (MatrixType::IsRowMajor) { + for (Index k = k0; k < k1; ++k) { + const Index p = transpositions.coeff(k); + if (p != k) mat.row(k).head(k0).swap(mat.row(p).head(k0)); + } + } else { + for (Index c = 0; c < k0; ++c) { + for (Index k = k0; k < k1; ++k) { + const Index p = transpositions.coeff(k); + if (p != k) numext::swap(mat.coeffRef(k, c), mat.coeffRef(p, c)); + } + } + } + + const Index rs = size - k1; + if (rs > 0) { + auto L21 = mat.block(k1, k0, rs, kb); + auto W = workspace.topLeftCorner(rs, kb); + W.noalias() = L21 * mat.diagonal().real().segment(k0, kb).asDiagonal(); + mat.block(k1, k1, rs, rs).template triangularView<Lower>() -= L21 * W.adjoint(); + } + } + return ret; + } + // Reference for the algorithm: Davis and Hager, "Multiple Rank // Modifications of a Sparse Cholesky Factorization" (Algorithm 1) // Trivial rearrangements of their computations (Timothy E. Holy) @@ -512,6 +626,13 @@ return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); } + template <typename MatrixType, typename TranspositionType, typename Workspace, typename BlockWorkspace> + static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, + BlockWorkspace& workspace, SignMatrix& sign) { + Transpose<MatrixType> matt(mat); + return ldlt_inplace<Lower>::blocked(matt, transpositions, temp, workspace, sign); + } + template <typename MatrixType, typename TranspositionType, typename Workspace, typename WType> static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma = 1) { @@ -556,8 +677,9 @@ m_temporary.resize(size); m_sign = internal::ZeroSign; - m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success - : NumericalIssue; + m_info = internal::ldlt_inplace<UpLo>::blocked(m_matrix, m_transpositions, m_temporary, m_workspace, m_sign) + ? Success + : NumericalIssue; m_isInitialized = true; return *this;
diff --git a/benchmarks/Cholesky/bench_cholesky.cpp b/benchmarks/Cholesky/bench_cholesky.cpp index a901aba..e010f6d 100644 --- a/benchmarks/Cholesky/bench_cholesky.cpp +++ b/benchmarks/Cholesky/bench_cholesky.cpp
@@ -15,20 +15,30 @@ static void BM_LDLT(benchmark::State& state) { int n = state.range(0); - typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType; - typedef Matrix<Scalar, Dynamic, Dynamic> SquareMatrixType; + using MatrixType = Matrix<Scalar, Dynamic, Dynamic>; MatrixType a = MatrixType::Random(n, n); - SquareMatrixType covMat = a * a.adjoint(); + MatrixType covMat = a * a.adjoint(); + // In place: each compute() copies covMat into the bound storage once and allocates no factorization storage. + MatrixType storage = covMat; + LDLT<Ref<MatrixType> > cholnosqrt(storage); int r = internal::random<int>(0, n - 1); int c = internal::random<int>(0, n - 1); Scalar acc = 0; for (auto _ : state) { - LDLT<SquareMatrixType> cholnosqrt(covMat); + cholnosqrt.compute(covMat); acc += cholnosqrt.matrixL().coeff(r, c); benchmark::DoNotOptimize(acc); } + double cost = 0; + for (int j = 0; j < n; ++j) { + int rem = std::max(n - j - 1, 0); + cost += 2 * (rem * j + rem + j); + } + state.counters["GFLOPS"] = + benchmark::Counter(cost, benchmark::Counter::kIsIterationInvariantRate, benchmark::Counter::kIs1000); } -BENCHMARK(BM_LDLT)->RangeMultiplier(2)->Range(4, 1500); +// 96 and 97 bracket the unblocked/blocked crossover. +BENCHMARK(BM_LDLT)->RangeMultiplier(2)->Range(4, 1500)->Arg(96)->Arg(97); static void BM_LLT(benchmark::State& state) { int n = state.range(0);
diff --git a/doc/TopicLinearAlgebraDecompositions.dox b/doc/TopicLinearAlgebraDecompositions.dox index 2280c67..78c2802 100644 --- a/doc/TopicLinearAlgebraDecompositions.dox +++ b/doc/TopicLinearAlgebraDecompositions.dox
@@ -144,7 +144,7 @@ <td>-</td> <td>Yes</td> <td>Excellent</td> - <td>-</td> + <td>Blocking</td> </tr> <tr> @@ -335,10 +335,10 @@ \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="note1">\b 1: </a>There exist two variants of the LDLT algorithm. %Eigen's LDLT produces a pure diagonal D matrix and therefore is not stable on indefinite matrices; it pivots on the diagonal entry of largest magnitude of the remaining Schur complement, which keeps the factorization stable when pivots are tiny or zero and puts the negligible pivots of a semidefinite matrix last. 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; \c PreconditionSquareMatrix applies it to square matrices too. 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="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: a symmetric factorization searches only the diagonal of the Schur complement for its pivot, which blocking leaves cheap to maintain. LAPACK thus provides a blocked complete-pivoting Cholesky (\c pstrf) for positive semidefinite matrices, where the pivoting also determines the rank, and %Eigen's LDLT pivots the same way.</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>
diff --git a/test/cholesky.cpp b/test/cholesky.cpp index d313dae..5b7a723 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp
@@ -486,21 +486,27 @@ VERIFY_IS_NOT_APPROX(mat, ldlt.reconstructedMatrix()); VERIFY(ldlt.info() == NumericalIssue); } + + // Pivoting on the input's diagonal took a zero of the Schur complement as a pivot of each of the next three + // indefinite matrices, above a non-zero column, and failed. Pivoting on the Schur complement's diagonal passes over + // that zero: they factor, and D carries their inertia (Sylvester's law of inertia). #if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE_SSE2) { mat.resize(3, 3); mat << -1, -3, 3, -3, -8.9999999999999999999, 1, 3, 1, 0; ldlt.compute(mat); - VERIFY(ldlt.info() == NumericalIssue); - VERIFY_IS_NOT_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY(ldlt.info() == Success); + VERIFY_IS_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY_IS_EQUAL((ldlt.vectorD().array() < 0).count(), Index(2)); } #endif { mat.resize(3, 3); mat << 1, 2, 3, 2, 4, 1, 3, 1, 0; ldlt.compute(mat); - VERIFY(ldlt.info() == NumericalIssue); - VERIFY_IS_NOT_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY(ldlt.info() == Success); + VERIFY_IS_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY_IS_EQUAL((ldlt.vectorD().array() < 0).count(), Index(1)); } { @@ -509,8 +515,9 @@ 0, 8.49333, 0, 2.00333, 0, 0, 0, 0, -0.1, 0, 0.1, 0, 0, 1, 0, 0, 0, 2.00333, 0, 4.24667, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0; ldlt.compute(mat); - VERIFY(ldlt.info() == NumericalIssue); - VERIFY_IS_NOT_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY(ldlt.info() == Success); + VERIFY_IS_APPROX(mat, ldlt.reconstructedMatrix()); + VERIFY_IS_EQUAL((ldlt.vectorD().array() < 0).count(), Index(2)); } // bug 1479 @@ -523,6 +530,228 @@ } } +// Semidefinite matrices of rank r whose arithmetic is exact, so that their zero pivots stay exactly zero. LDLT pivots +// on the Schur complement's diagonal, so the n - r zeros of D come last. Pivoting on the input's diagonal took the zero +// of [1 1 0; 1 1 0; 0 0 1] second, left it inside D, and reported NumericalIssue. +template <typename MatrixType, int UpLo> +void check_ldlt_zero_pivots_last(const MatrixType& A, Index r) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + using RealVectorType = Matrix<RealScalar, Dynamic, 1>; + const Index n = A.rows(); + const MatrixType triangle = A.template triangularView<UpLo>(); + LDLT<MatrixType, UpLo> ldlt(triangle); + VERIFY(ldlt.info() == Success); + VERIFY(ldlt.isPositive()); + const RealVectorType d = ldlt.vectorD().real(); + VERIFY_IS_CWISE_EQUAL(d.head(r), RealVectorType::Ones(r)); + VERIFY_IS_CWISE_EQUAL(d.tail(n - r), RealVectorType::Zero(n - r)); + VERIFY_IS_CWISE_EQUAL(ldlt.reconstructedMatrix(), A); +} + +template <typename MatrixType> +void cholesky_ldlt_zero_pivots_last(const MatrixType& A, Index r) { + check_ldlt_zero_pivots_last<MatrixType, Lower>(A, r); + check_ldlt_zero_pivots_last<MatrixType, Upper>(A, r); +} + +// The panel width at the smallest size that LDLT factors in panels, or 0 if it factors none, as with +// EIGEN_LDLT_BLOCKSIZE=0. The width grows with the size, so the blocked tests below take their sizes from this one. +template <typename MatrixType> +Index ldlt_smallest_blocksize() { + for (Index n = 1; n <= 4096; ++n) + if (internal::ldlt_blocksize<MatrixType>(n) > 0) return internal::ldlt_blocksize<MatrixType>(n); + return 0; +} + +// P (I_m kron [1 1; 1 1]) P^T has rank m, at a size that takes the blocked factorization. +template <typename MatrixType> +void cholesky_ldlt_zero_pivots_last_blocked() { + const Index nb = ldlt_smallest_blocksize<MatrixType>(); + if (nb == 0) return; + const Index m = 2 * nb + 1, n = 2 * m; + VERIFY(internal::ldlt_blocksize<MatrixType>(n) > 0); + MatrixType A = MatrixType::Zero(n, n); + for (Index i = 0; i < m; ++i) A.template block<2, 2>(2 * i, 2 * i).setOnes(); + PermutationMatrix<Dynamic> perm(n); + perm.setIdentity(); + for (Index i = n - 1; i > 0; --i) + std::swap(perm.indices().coeffRef(i), perm.indices().coeffRef(internal::random<Index>(0, i))); + cholesky_ldlt_zero_pivots_last(MatrixType(perm * A * perm.transpose()), m); +} + +// For a positive semidefinite A of rank r, the trailing n - r pivots of complete diagonal pivoting are at rounding +// level and the leading r are not (N. J. Higham, "Analysis of the Cholesky decomposition of a semi-definite matrix", +// 1990, for LAPACK's xPSTRF). +template <typename MatrixType, int UpLo> +void check_ldlt_semidefinite(const MatrixType& A, Index r) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + const Index n = A.rows(); + const RealScalar eps = NumTraits<RealScalar>::epsilon(); + const MatrixType triangle = A.template triangularView<UpLo>(); + // info() is not checked: rounding can leave an exactly zero pivot above a column of rounding errors, which + // NumericalIssue reports as it would for an indefinite matrix. + LDLT<MatrixType, UpLo> ldlt(triangle); + + // Complete pivoting keeps |L| <= 1 on a semidefinite matrix, so each trailing pivot is a sum of about n rounding + // errors of size eps d_0. Over 20 seeds and n <= 300, in single and double precision, they stayed below + // 1.2 n eps d_0, and the reconstruction error below 0.25 n eps |A|. With r = 1 or n/2 random columns in B, the + // leading pivots stay far above the tail bound. + const Matrix<RealScalar, Dynamic, 1> d = ldlt.vectorD().real().cwiseAbs(); + const RealScalar tail_bound = RealScalar(8 * n) * eps * d(0); + if (r < n) VERIFY(d.tail(n - r).maxCoeff() <= tail_bound); + VERIFY(d.head(r).minCoeff() > tail_bound); + VERIFY((ldlt.reconstructedMatrix() - A).norm() <= RealScalar(2 * n) * eps * A.norm()); +} + +// The two rows of B of largest norm are parallel, so that pivoting by the input's diagonal takes a rounding-level pivot +// second. +template <typename MatrixType> +void cholesky_ldlt_semidefinite(Index n, Index r) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + MatrixType B = MatrixType::Random(n, r); + if (n >= 2) { + B.row(0) *= RealScalar(4); + B.row(1) = RealScalar(0.5) * B.row(0); + } + const MatrixType A = B * B.adjoint(); + check_ldlt_semidefinite<MatrixType, Lower>(A, r); + check_ldlt_semidefinite<MatrixType, Upper>(A, r); +} + +// A = G M G with M = B B^* + I and G a shuffled diagonal over 2 * decades orders of magnitude. Complete diagonal +// pivoting follows the grading, so the factorization is componentwise accurate, +// |A - P^T L D L^* P|_ij <= c n eps sqrt(a_ii a_jj). Forming |L_ik|^2 before D_k |L_ik|^2 in the diagonal downdate +// underflowed for |L_ik| < sqrt(min) and lost that accuracy. +template <typename MatrixType, int UpLo> +void check_ldlt_graded(const MatrixType& A) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + using RealVectorType = Matrix<RealScalar, Dynamic, 1>; + const Index n = A.rows(); + const RealScalar eps = NumTraits<RealScalar>::epsilon(); + const MatrixType triangle = A.template triangularView<UpLo>(); + LDLT<MatrixType, UpLo> ldlt(triangle); + VERIFY(ldlt.info() == Success); + const RealVectorType root = A.diagonal().real().cwiseSqrt(); + const RealScalar error = + ((ldlt.reconstructedMatrix() - A).cwiseAbs().array() / (root * root.transpose()).array()).maxCoeff(); + // Measured below 0.1 n eps, against above 1e10 n eps with the underflowing downdate. + VERIFY(error <= RealScalar(2 * n) * eps); +} + +template <typename MatrixType> +void cholesky_ldlt_graded(Index n, int decades) { + using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real; + Matrix<RealScalar, Dynamic, 1> g(n); + for (Index i = 0; i < n; ++i) + g(i) = std::pow(RealScalar(10), RealScalar(decades) * (RealScalar(2 * i) / RealScalar(n - 1) - RealScalar(1))); + for (Index i = n - 1; i > 0; --i) std::swap(g(i), g(internal::random<Index>(0, i))); + const MatrixType B = MatrixType::Random(n, n); + const MatrixType A = g.asDiagonal() * (B * B.adjoint() + MatrixType::Identity(n, n)) * g.asDiagonal(); + check_ldlt_graded<MatrixType, Lower>(A); + check_ldlt_graded<MatrixType, Upper>(A); +} + +// Definite and KKT matrices of size n, which give pivots of both signs. +template <typename MatrixType> +void check_ldlt_blocked(Index n) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename NumTraits<Scalar>::Real; + using RealMatrixType = Matrix<RealScalar, Dynamic, Dynamic>; + const RealScalar eps = NumTraits<RealScalar>::epsilon(); + const MatrixType B = MatrixType::Random(n, n); + const MatrixType pd = B * B.adjoint() + MatrixType::Identity(n, n); + for (RealScalar sign : {RealScalar(1), RealScalar(-1)}) { + const MatrixType A = sign * pd; + LDLT<MatrixType> ldlt(A); + VERIFY(ldlt.info() == Success); + VERIFY(ldlt.isPositive() == (sign > 0) && ldlt.isNegative() == (sign < 0)); + VERIFY((ldlt.reconstructedMatrix() - A).norm() <= RealScalar(2 * n) * eps * A.norm()); + } + + // [H J^*; J 0] with H positive definite and J of full row rank k has k negative eigenvalues. The matrix is + // indefinite, so the backward error is bounded by |L| |D| |L^*| rather than by |A|: |dA| <= c n eps |L| |D| |L^*| + // (Higham, ASNA 2nd ed., section 11.1). Measured below 0.01 n eps |A|, with |L| |D| |L^*| <= 3 |A|. + const Index k = n / 4; + MatrixType kkt = MatrixType::Zero(n, n); + kkt.topLeftCorner(n - k, n - k) = pd.topLeftCorner(n - k, n - k); + kkt.bottomLeftCorner(k, n - k) = MatrixType::Random(k, n - k); + kkt.topRightCorner(n - k, k) = kkt.bottomLeftCorner(k, n - k).adjoint(); + LDLT<MatrixType, Upper> ldlt(kkt); + VERIFY(ldlt.info() == Success); + const RealMatrixType absL = MatrixType(ldlt.matrixL()).cwiseAbs(); + const RealScalar growth = (absL * ldlt.vectorD().real().cwiseAbs().asDiagonal() * absL.transpose()).norm(); + VERIFY((ldlt.reconstructedMatrix() - kkt).norm() <= RealScalar(2 * n) * eps * growth); + VERIFY_IS_EQUAL((ldlt.vectorD().real().array() < 0).count(), k); +} + +// Above three panel widths, LDLT factors panels and applies each to the trailing Schur complement in one level-3 +// product. At the smallest panel width, the sizes leave a trailing panel of one column, none, and a partial one; a +// fourth size is the first with a wider panel. Fixed-size and bounded matrices stay unblocked, so that they allocate +// no panel workspace. +template <typename MatrixType> +void cholesky_ldlt_blocked() { + using Scalar = typename MatrixType::Scalar; + const Index nb = ldlt_smallest_blocksize<MatrixType>(); + if (nb == 0) return; + VERIFY_IS_EQUAL(internal::ldlt_blocksize<MatrixType>(3 * nb), Index(0)); + VERIFY_IS_EQUAL((internal::ldlt_blocksize<Matrix<Scalar, Dynamic, Dynamic, 0, 512, 512> >(4 * nb + 1)), Index(0)); + for (Index n : {4 * nb + 1, 5 * nb, 5 * nb + nb / 2}) { + VERIFY_IS_EQUAL(internal::ldlt_blocksize<MatrixType>(n), nb); + check_ldlt_blocked<MatrixType>(n); + } +#ifndef EIGEN_LDLT_BLOCKSIZE + Index n = 5 * nb + nb / 2; + while (internal::ldlt_blocksize<MatrixType>(n) == nb) ++n; + check_ldlt_blocked<MatrixType>(n); +#endif +} + +// A NaN or an infinity anywhere in the matrix reaches a pivot, so the factorization reports it instead of Success. +template <typename MatrixType> +void check_ldlt_nonfinite(Index n, const typename MatrixType::Scalar& value) { + const Index positions[][2] = {{0, 0}, {n / 2, n / 2}, {n - 1, n - 1}, {n - 1, 0}, {n - 1, n - 2}}; + for (const auto& position : positions) { + const Index i = position[0], j = position[1]; + if (j < 0 || (i == j && NumTraits<typename MatrixType::Scalar>::IsComplex && numext::real(value) == 0)) continue; + MatrixType A = MatrixType::Identity(n, n); + A(i, j) = value; + A(j, i) = numext::conj(value); + VERIFY((LDLT<MatrixType, Lower>(A).info() == NumericalIssue)); + VERIFY((LDLT<MatrixType, Upper>(A).info() == NumericalIssue)); + } +} + +template <typename MatrixType> +void cholesky_ldlt_nonfinite(Index n) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename NumTraits<Scalar>::Real; + check_ldlt_nonfinite<MatrixType>(n, Scalar(std::numeric_limits<RealScalar>::quiet_NaN())); + check_ldlt_nonfinite<MatrixType>(n, Scalar(std::numeric_limits<RealScalar>::infinity())); +} + +// A NaN in the imaginary part of an off-diagonal entry; the diagonal of a Hermitian matrix has none to read. +template <typename MatrixType> +void cholesky_ldlt_imaginary_nan(Index n) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename NumTraits<Scalar>::Real; + check_ldlt_nonfinite<MatrixType>(n, Scalar(RealScalar(0), std::numeric_limits<RealScalar>::quiet_NaN())); +} + +// The blocked factorization takes its panel workspace from the decomposition, so a preallocated one allocates nothing +// at this size, where the trailing product still packs its operands on the stack. +template <typename Scalar> +void cholesky_ldlt_blocked_preallocated_no_malloc() { + using MatrixType = Matrix<Scalar, Dynamic, Dynamic>; + const Index size = 4 * numext::maxi(ldlt_smallest_blocksize<MatrixType>(), Index(32)) + 3; + const MatrixType A = MatrixType::Random(size, size); + const MatrixType spd = A * A.adjoint() + MatrixType::Identity(size, size); + LDLT<MatrixType> ldlt(size); + internal::set_is_malloc_allowed(false); + ldlt.compute(spd); + internal::set_is_malloc_allowed(true); + VERIFY_IS_EQUAL(ldlt.info(), Success); +} + // Accumulates I + sum_k w_k w_k^* out of rank-1 updates alone, which is the path rankUpdate() takes while the // decomposition holds no factorization, and checks the reported status against the factorization it produced. template <typename MatrixType, int UpLo> @@ -998,6 +1227,40 @@ CALL_SUBTEST_2(cholesky_faillure_cases<void>()); + CALL_SUBTEST_3(cholesky_ldlt_zero_pivots_last((Matrix3d() << 1, 1, 0, 1, 1, 0, 0, 0, 1).finished(), 2)); + CALL_SUBTEST_2(cholesky_ldlt_zero_pivots_last(MatrixXd((Matrix3d() << 1, 1, 0, 1, 1, 0, 0, 0, 1).finished()), 2)); + CALL_SUBTEST_6(cholesky_ldlt_zero_pivots_last( + (Matrix3cd() << 1, std::complex<double>(0, 1), 0, std::complex<double>(0, -1), 1, 0, 0, 0, 1).finished(), 2)); + CALL_SUBTEST_2(cholesky_ldlt_zero_pivots_last_blocked<MatrixXd>()); + CALL_SUBTEST_6(cholesky_ldlt_zero_pivots_last_blocked<MatrixXcd>()); + + CALL_SUBTEST_2(cholesky_ldlt_blocked<MatrixXd>()); + CALL_SUBTEST_2((cholesky_ldlt_blocked<Matrix<double, Dynamic, Dynamic, RowMajor> >())); + CALL_SUBTEST_4(cholesky_ldlt_blocked<MatrixXf>()); + CALL_SUBTEST_6(cholesky_ldlt_blocked<MatrixXcd>()); + CALL_SUBTEST_1(cholesky_ldlt_blocked_preallocated_no_malloc<double>()); + CALL_SUBTEST_1(cholesky_ldlt_blocked_preallocated_no_malloc<std::complex<double> >()); + + { + // Two unblocked sizes and two blocked ones, unless EIGEN_LDLT_BLOCKSIZE=0 disables blocking. + const Index nb = numext::maxi(ldlt_smallest_blocksize<MatrixXd>(), Index(32)); + for (Index n : {Index(3), Index(17), 4 * nb + 1, 5 * nb + nb / 2}) { + CALL_SUBTEST_2(cholesky_ldlt_semidefinite<MatrixXd>(n, 1)); + CALL_SUBTEST_2(cholesky_ldlt_semidefinite<MatrixXd>(n, n / 2)); + CALL_SUBTEST_4(cholesky_ldlt_semidefinite<MatrixXf>(n, n / 2)); + CALL_SUBTEST_6(cholesky_ldlt_semidefinite<MatrixXcd>(n, n / 2)); + CALL_SUBTEST_2(cholesky_ldlt_graded<MatrixXd>(n, 80)); + CALL_SUBTEST_4(cholesky_ldlt_graded<MatrixXf>(n, 12)); + CALL_SUBTEST_6(cholesky_ldlt_graded<MatrixXcd>(n, 80)); + CALL_SUBTEST_2(cholesky_ldlt_nonfinite<MatrixXd>(n)); + CALL_SUBTEST_6(cholesky_ldlt_nonfinite<MatrixXcd>(n)); + CALL_SUBTEST_6(cholesky_ldlt_imaginary_nan<MatrixXcd>(n)); + TEST_SET_BUT_UNUSED_VARIABLE(n); + } + CALL_SUBTEST_2(cholesky_ldlt_nonfinite<MatrixXd>(1)); + CALL_SUBTEST_5(cholesky_ldlt_nonfinite<Matrix4d>(4)); + } + // Blocking and vectorization boundary tests (deterministic, outside g_repeat). CALL_SUBTEST_2(cholesky_blocking_boundary<double>()); CALL_SUBTEST_8(cholesky_blocking_boundary<float>());