StructuredMatrices: Add DiagonalPlusLowRank operator with Woodbury solve libeigen/eigen!2693 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com> Co-authored-by: Rasmus Munk Larsen <rlarsen@nvidia.com>
diff --git a/unsupported/Eigen/StructuredMatrices b/unsupported/Eigen/StructuredMatrices index fabe21e..4e9d160 100644 --- a/unsupported/Eigen/StructuredMatrices +++ b/unsupported/Eigen/StructuredMatrices
@@ -34,7 +34,9 @@ * matrices as an implicit operator; products, solves, least squares, * eigendecomposition, SVD, inverse and determinant all factor through the * operands, and diagonal (in particular identity) factors are stored and - * applied in diagonal form. + * applied in diagonal form; + * - \c DiagonalPlusLowRank : a diagonal matrix plus a rank-k correction, with + * O(nk) products and O(nk^2) Woodbury solves, closed under inversion. * * The operator types derive from \c EigenBase and store only compact generators * or factors. The FFT-backed operators (\c Circulant, \c Toeplitz and \c Hankel) @@ -73,6 +75,7 @@ #include "src/StructuredMatrices/LookAheadLevinson.h" #include "src/StructuredMatrices/Hankel.h" #include "src/StructuredMatrices/KroneckerOperator.h" +#include "src/StructuredMatrices/DiagonalPlusLowRank.h" // IWYU pragma: end_exports #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
diff --git a/unsupported/Eigen/src/StructuredMatrices/DiagonalPlusLowRank.h b/unsupported/Eigen/src/StructuredMatrices/DiagonalPlusLowRank.h new file mode 100644 index 0000000..18df0a0 --- /dev/null +++ b/unsupported/Eigen/src/StructuredMatrices/DiagonalPlusLowRank.h
@@ -0,0 +1,459 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +// References: +// [1] M. A. Woodbury, "Inverting Modified Matrices", Memorandum Report 42, +// Statistical Research Group, Princeton University, 1950. The Woodbury +// identity behind solve() and inverse(). +// [2] W. W. Hager, "Updating the Inverse of a Matrix", SIAM Review, 31(2), +// pp. 221-239, 1989. Review of the Woodbury identity and of the matrix +// determinant lemma used by determinant(). +// [3] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart, "LINPACK +// Users' Guide", SIAM, 1979. determinant()'s balanced accumulation follows +// the convention of its xGEDI routines, which return determinants as a +// (fraction, exponent) pair to avoid spurious overflow/underflow. +// [4] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974. +// Scaling by a power of two is exact, the property the balanced +// accumulation and the normalized Woodbury products rely on. +// [5] N. J. Higham, "Accuracy and Stability of Numerical Algorithms", 2nd ed., +// SIAM, 2002, chapter 27. Overflow-avoiding rescaling of intermediate +// quantities, the technique behind the normalized capacitance products. +// [6] E. L. Yip, "A Note on the Stability of Solving a Rank-p Modification of +// a Linear System by the Sherman-Morrison-Woodbury Formula", SIAM Journal +// on Scientific and Statistical Computing, 7(3), pp. 507-513, 1986. The +// accuracy limitation of the Woodbury solve noted on solve(). + +#ifndef EIGEN_STRUCTURED_DIAGONAL_PLUS_LOW_RANK_H +#define EIGEN_STRUCTURED_DIAGONAL_PLUS_LOW_RANK_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +template <typename Scalar_, int Size_ = Dynamic, int Rank_ = Dynamic> +class DiagonalPlusLowRank; + +namespace internal { + +template <typename Scalar_, int Size_, int Rank_> +struct traits<DiagonalPlusLowRank<Scalar_, Size_, Rank_>> { + using Scalar = Scalar_; + using StorageKind = Dense; + using XprKind = MatrixXpr; + using StorageIndex = int; + static constexpr int RowsAtCompileTime = Size_; + static constexpr int ColsAtCompileTime = Size_; + static constexpr int MaxRowsAtCompileTime = Size_; + static constexpr int MaxColsAtCompileTime = Size_; + // Deliberately no NestByRefBit: transpose(), conjugate(), adjoint(), inverse() + // and makeDiagonalPlusLowRank() return owning temporaries, so Product must nest + // the operator by value for a delayed-evaluated product expression to keep its + // left factor alive. The copy is O(nk), on par with a single product evaluation. + static constexpr unsigned int Flags = 0; +}; + +template <typename Scalar_, int Size_, int Rank_> +struct evaluator_traits<DiagonalPlusLowRank<Scalar_, Size_, Rank_>> { + using Kind = IndexBased; + using Shape = StructuredShape; +}; + +// Entrywise multiplication by the exact power of two 2^e [4], for the factor normalization of the +// Woodbury products [5]. Uses structured_balance_impl::apply_exponent (an ldexp per component) rather +// than a multiplication by 2^e or by two half powers 2^(e/2): ldexp is exact and saturates entrywise +// for every exponent, whereas 2^-eu itself overflows for a subnormal factor bound and a saturated +// half power turns zero entries into NaN through 0 * Inf. +template <typename Scalar> +struct dplr_ldexp_op { + int e; + Scalar operator()(const Scalar& x) const { return structured_balance_impl<Scalar>::apply_exponent(x, e); } +}; + +// The smallest e with n <= 2^e: the inner-dimension term of a product's entry +// magnitude bound (a length-n dot product of entries below 2^a and 2^b stays +// below 2^(a + b + e)). +inline int dplr_index_exponent(Index n) { return log2_ceil(static_cast<std::make_unsigned_t<Index>>(n)); } + +// True when a product of operands with entry-magnitude exponent bounds ea and eb over an inner +// dimension of length \a inner provably cannot overflow, with one bit of headroom for a trailing +// addition. Where it fails the Woodbury kernels rescale their factors exactly [4][5]; where it holds +// they keep the plain association. +template <typename RealScalar> +bool dplr_product_fits(int ea, int eb, Index inner) { + return ea + eb + dplr_index_exponent(inner) + 1 < NumTraits<RealScalar>::max_exponent(); +} + +// Compile-time dispatch keeps fixed rank zero from instantiating a 0 x 0 LU; +// C++14 cannot express this branch with if constexpr. +template <int Rank_> +struct dplr_capacitance_impl { + /** Applies the Woodbury correction + * \f[ x\mathrel{-}=D^{-1}U(I_k+V^H D^{-1}U)^{-1}V^Hx. \f] + * Exact power-of-two scaling is used only when either structural product can + * overflow [4][5], preserving the unscaled result bit-for-bit otherwise. */ + template <typename Op, typename Dinv, typename Workspace> + static void subtractSolveCorrection(const Op& op, const Dinv& dinv, Workspace& x) { + using Scalar = typename Op::Scalar; + using RealScalar = typename Op::RealScalar; + if (op.correctionRank() == 0 || x.size() == 0) return; + PartialPivLU<typename Op::CapacitanceType> cap(op.capacitance()); + const int eu = structured_exponent_bound(op.factorU()); + const int ev = structured_exponent_bound(op.factorV()); + const int ed = structured_exponent_bound(dinv); + const int ex = structured_exponent_bound(x); + if (dplr_product_fits<RealScalar>(ev, ex, op.rows()) && dplr_product_fits<RealScalar>(ed, eu, 1)) { + x.noalias() -= dinv.asDiagonal() * (op.factorU() * cap.solve(op.factorV().adjoint() * x)); + return; + } + Matrix<Scalar, Rank_, Workspace::ColsAtCompileTime> y = + cap.solve(op.factorV().unaryExpr(dplr_ldexp_op<Scalar>{-ev}).adjoint() * x); + y = y.unaryExpr(dplr_ldexp_op<Scalar>{eu + ev}); + x.noalias() -= dinv.asDiagonal() * (op.factorU().unaryExpr(dplr_ldexp_op<Scalar>{-eu}) * y); + } + /** Forms \f$U'=-D^{-1}U(I_k+V^HD^{-1}U)^{-1}\f$. Scaling prevents a + * spurious overflow in \f$D^{-1}U\f$ when the capacitance inverse subsequently + * shrinks it [4][5]. */ + template <typename Op, typename Dinv, typename Factor> + static void inverseFactor(const Op& op, const Dinv& dinv, Factor& Up) { + using Scalar = typename Op::Scalar; + using RealScalar = typename Op::RealScalar; + if (op.correctionRank() == 0) return; + PartialPivLU<typename Op::CapacitanceType> cap(op.capacitance()); + const typename Op::CapacitanceType K = cap.inverse(); + const int eu = structured_exponent_bound(op.factorU()); + const int ed = structured_exponent_bound(dinv); + const int ek = structured_exponent_bound(K); + if (dplr_product_fits<RealScalar>(ed, eu, 1) && dplr_product_fits<RealScalar>(ed + eu, ek, op.correctionRank())) { + Up.noalias() = -(dinv.asDiagonal() * op.factorU() * K); + return; + } + Up.noalias() = -(dinv.asDiagonal() * op.factorU().unaryExpr(dplr_ldexp_op<Scalar>{-eu}) * K); + Up = Up.unaryExpr(dplr_ldexp_op<Scalar>{eu}); + } + /** Accumulates the determinant-lemma factor + * \f$\det(I_k+V^HD^{-1}U)\f$ as a mantissa and power of two. Each LU pivot is + * renormalized before multiplication because \f$\det(D)\f$ and the capacitance + * determinant can occupy opposite ends of the exponent range. */ + template <typename Op> + static typename Op::Scalar balancedCapacitanceDeterminant(const Op& op, Index& exponent) { + using Scalar = typename Op::Scalar; + if (op.correctionRank() == 0) return Scalar(1); + const PartialPivLU<typename Op::CapacitanceType> lu(op.capacitance()); + Scalar det(static_cast<typename Op::RealScalar>(lu.permutationP().determinant())); + for (Index i = 0; i < lu.matrixLU().rows(); ++i) + det = structured_balance(det * structured_balance(lu.matrixLU().coeff(i, i), exponent), exponent); + return det; + } +}; + +template <> +struct dplr_capacitance_impl<0> { + template <typename Op, typename Dinv, typename Workspace> + static void subtractSolveCorrection(const Op&, const Dinv&, Workspace&) {} + template <typename Op, typename Dinv, typename Factor> + static void inverseFactor(const Op&, const Dinv&, Factor&) {} + template <typename Op> + static typename Op::Scalar balancedCapacitanceDeterminant(const Op&, Index&) { + return typename Op::Scalar(1); + } +}; + +} // namespace internal + +/** \ingroup StructuredMatrices_Module + * \class DiagonalPlusLowRank + * \brief An \c n x \c n operator \f$ D + U V^H \f$: a diagonal matrix plus a + * rank-k correction, stored as its diagonal and the two \c n x \c k factors. + * + * Products cost O(nk) instead of O(n^2). Linear systems are solved in O(nk^2) + * through the Woodbury identity [1] + * \f$ (D + UV^H)^{-1} = D^{-1} - D^{-1} U (I_k + V^H D^{-1} U)^{-1} V^H D^{-1} \f$, + * factoring only the \c k x \c k \em capacitance matrix (\ref solve). The + * determinant follows from the matrix determinant lemma [2] + * \f$ \det(D)\,\det(I_k + V^H D^{-1} U) \f$ (\ref determinant), and the class is + * closed under \ref inverse, \ref transpose, \ref conjugate and \ref adjoint. + * + * The rank-k correction may have \c k = 0 (a plain diagonal operator). \c k is + * not required to be small, but every advantage over a dense matrix vanishes as + * \c k approaches \c n. + * + * The operator stores its own copies of the diagonal and the factors and derives + * from \c EigenBase. Because \c operator* returns an Eigen product expression, a + * \c DiagonalPlusLowRank also drops into the matrix-free iterative solvers, and + * it can be assigned to a dense matrix when an explicit representation is + * needed. As with any matrix-free operator, the iterative solvers must be + * instantiated with \c IdentityPreconditioner (e.g. + * \c GMRES<DiagonalPlusLowRank<double>,IdentityPreconditioner>): the default + * preconditioners read individual coefficients through \c col() or + * \c InnerIterator, which the structured operators do not expose. + * + * \tparam Scalar_ the scalar type, real or complex. + * \tparam Size_ the dimension at compile time, or \c Dynamic (the default). + * \tparam Rank_ the correction rank at compile time, or \c Dynamic (the default). + * + * \sa makeDiagonalPlusLowRank() + */ +template <typename Scalar_, int Size_, int Rank_> +class DiagonalPlusLowRank : public EigenBase<DiagonalPlusLowRank<Scalar_, Size_, Rank_>> { + public: + using Scalar = Scalar_; + using RealScalar = typename NumTraits<Scalar>::Real; + using StorageIndex = int; + using DiagonalVector = Matrix<Scalar, Size_, 1>; + using FactorType = Matrix<Scalar, Size_, Rank_>; + using CapacitanceType = Matrix<Scalar, Rank_, Rank_>; + + static constexpr int RowsAtCompileTime = Size_; + static constexpr int ColsAtCompileTime = Size_; + static constexpr int MaxRowsAtCompileTime = Size_; + static constexpr int MaxColsAtCompileTime = Size_; + static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Size_, Size_); + static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime; + static constexpr bool IsRowMajor = false; + // Deliberately no IsVectorAtCompileTime: Ref<const DiagonalPlusLowRank>'s default + // StrideType argument reads it, so its absence makes internal::is_ref_compatible + // SFINAE to false and keeps the iterative solvers on their matrix-free path. + + /** Builds the operator \c diag(d) + \c U*V.adjoint(). The factors must have the + * same number of columns (the correction rank \c k, possibly zero) and as many + * rows as \a d has entries. */ + template <typename DDerived, typename UDerived, typename VDerived> + DiagonalPlusLowRank(const MatrixBase<DDerived>& d, const MatrixBase<UDerived>& U, const MatrixBase<VDerived>& V) + : m_d(d), m_U(U), m_V(V) { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(DDerived) + eigen_assert(m_d.size() > 0 && "DiagonalPlusLowRank must be non-empty"); + eigen_assert(m_U.rows() == m_d.size() && m_V.rows() == m_d.size() && m_U.cols() == m_V.cols() && + "factor dimensions do not match"); + } + + EIGEN_DEVICE_FUNC Index rows() const { return m_d.size(); } + EIGEN_DEVICE_FUNC Index cols() const { return m_d.size(); } + /** \returns the correction rank \c k. */ + Index correctionRank() const { return m_U.cols(); } + + /** \returns the diagonal vector \c d. */ + const DiagonalVector& diagonal() const { return m_d; } + /** \returns the left factor \c U. */ + const FactorType& factorU() const { return m_U; } + /** \returns the right factor \c V. */ + const FactorType& factorV() const { return m_V; } + + /** \returns the coefficient at row \a row and column \a col. */ + Scalar coeff(Index row, Index col) const { + Scalar r = m_V.row(col).dot(m_U.row(row)); // dot() conjugates its first argument + if (row == col) r += m_d.coeff(row); + return r; + } + + /** \returns the transpose of \c *this, itself a \c DiagonalPlusLowRank + * operator: \f$ (D + UV^H)^T = D + \bar V \bar U^H \f$. */ + DiagonalPlusLowRank transpose() const { return DiagonalPlusLowRank(m_d, m_V.conjugate(), m_U.conjugate()); } + + /** \returns the complex conjugate of \c *this, itself a \c DiagonalPlusLowRank + * operator. */ + DiagonalPlusLowRank conjugate() const { + return DiagonalPlusLowRank(m_d.conjugate(), m_U.conjugate(), m_V.conjugate()); + } + + /** \returns the adjoint of \c *this, itself a \c DiagonalPlusLowRank operator: + * \f$ (D + UV^H)^H = \bar D + V U^H \f$. */ + DiagonalPlusLowRank adjoint() const { return DiagonalPlusLowRank(m_d.conjugate(), m_V, m_U); } + + /** \returns the capacitance matrix \f$ I_k + V^H D^{-1} U \f$ of the Woodbury + * identity. Every consumer of the triple product -- \ref solve, \ref inverse + * and \ref determinant -- forms it through this one method. + * + * With extreme factor magnitudes the plain association can overflow even + * though the capacitance itself is representable: for \c d = 1e-200, + * \c U = 1e-200, \c V = 1e200 the operator is essentially the identity and + * the capacitance \c 1 + 1e200 is representable, but \f$ V^H D^{-1} \f$ is + * \c 1e400. The factors are therefore rescaled by exact powers of two [4] + * when a conservative exponent bound detects the danger, the product is + * formed as \f$ \hat V^H (D^{-1} \hat U) \f$ -- whose intermediates are + * bounded by roughly \c n * max|1/d| -- and the removed exponent is folded + * back entrywise before the identity is added [5]. When the plain + * association provably cannot overflow it is kept, so results for moderate + * data are bit-identical to the unnormalized evaluation. + * \warning The diagonal must have no zero entries. */ + CapacitanceType capacitance() const { + const Index k = correctionRank(); + CapacitanceType c = CapacitanceType::Identity(k, k); + if (k == 0) return c; + const DiagonalVector dinv = m_d.cwiseInverse(); + const int eu = internal::structured_exponent_bound(m_U); + const int ev = internal::structured_exponent_bound(m_V); + const int ed = internal::structured_exponent_bound(dinv); + if (internal::dplr_product_fits<RealScalar>(ev, ed, 1) && + internal::dplr_product_fits<RealScalar>(ev + ed, eu, rows())) { + c.noalias() += m_V.adjoint() * dinv.asDiagonal() * m_U; + return c; + } + const FactorType W = dinv.asDiagonal() * m_U.unaryExpr(internal::dplr_ldexp_op<Scalar>{-eu}); + CapacitanceType t(k, k); + t.noalias() = m_V.unaryExpr(internal::dplr_ldexp_op<Scalar>{-ev}).adjoint() * W; + t = t.unaryExpr(internal::dplr_ldexp_op<Scalar>{eu + ev}); + c += t; + return c; + } + + /** \returns the solution of \c (*this) * x = b through the Woodbury identity + * [1], factoring only the \c k x \c k capacitance matrix: O(nk^2 + k^3) setup + * and O(nk) per right-hand side. Supports multiple right-hand sides. + * \warning Requires an invertible diagonal \em and an invertible capacitance + * matrix (equivalently, an invertible operator); this is not checked beyond + * the NaN/Inf propagation of the arithmetic itself. + * \warning The Woodbury splitting routes the solution through \f$ D^{-1} b \f$: + * when \c max|1/d| greatly exceeds the norm of the operator's inverse, the + * correction cancels most of those amplified digits and the achievable + * accuracy degrades accordingly, even for a well-conditioned operator [6]. + * The normalized kernels keep such solves finite; they cannot restore the + * cancelled digits. */ + template <typename Rhs> + Matrix<Scalar, Size_, Rhs::ColsAtCompileTime> solve(const MatrixBase<Rhs>& b) const { + EIGEN_STATIC_ASSERT(RowsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic || + int(RowsAtCompileTime) == int(Rhs::RowsAtCompileTime), + YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES) + eigen_assert(b.rows() == rows() && "right-hand side has the wrong number of rows"); + const DiagonalVector dinv = m_d.cwiseInverse(); + // D^{-1} b can overflow while the solution is representable: at d = 1e-200, b = 1e200 the operator + // is essentially the identity and x is 1e200, yet this term is 1e400. The solve is linear in b, so + // b is rescaled by an exact power of two [4] on the same exponent test the other kernels use and + // the exponent is folded back into the result [5]; the plain form is kept where it provably + // cannot overflow. + const int eb = internal::structured_exponent_bound(b.derived()); + const int ed = internal::structured_exponent_bound(dinv); + const bool rescale = !internal::dplr_product_fits<RealScalar>(ed, eb, 1); + + Matrix<Scalar, Size_, Rhs::ColsAtCompileTime> x = + rescale ? (dinv.asDiagonal() * b.derived().unaryExpr(internal::dplr_ldexp_op<Scalar>{-eb})).eval() + : (dinv.asDiagonal() * b).eval(); + internal::dplr_capacitance_impl<Rank_>::subtractSolveCorrection(*this, dinv, x); + if (rescale) x = x.unaryExpr(internal::dplr_ldexp_op<Scalar>{eb}); + return x; + } + + /** \returns the inverse of \c *this, itself a \c DiagonalPlusLowRank operator: + * by the Woodbury identity [1], \f$ (D + UV^H)^{-1} = D^{-1} + U' V'^H \f$ with + * \f$ U' = -D^{-1} U (I_k + V^H D^{-1} U)^{-1} \f$ and \f$ V' = D^{-H} V \f$. + * \warning Same invertibility requirements as \ref solve. */ + DiagonalPlusLowRank inverse() const { + const DiagonalVector dinv = m_d.cwiseInverse(); + FactorType Up(rows(), correctionRank()); + internal::dplr_capacitance_impl<Rank_>::inverseFactor(*this, dinv, Up); + FactorType Vp = dinv.conjugate().asDiagonal() * m_V; + return DiagonalPlusLowRank(dinv, Up, Vp); + } + + /** \returns the determinant through the matrix determinant lemma [2]: + * \f$ \det(D)\,\det(I_k + V^H D^{-1} U) \f$, in O(nk^2 + k^3) operations. + * Both factors are accumulated in the balanced form \c m * 2^e (the split + * fraction/exponent determinant convention of LINPACK's xGEDI [3]) -- the + * diagonal entries and the LU pivots of the capacitance matrix are + * renormalized to unit magnitude one at a time, with the powers of two tracked + * in a shared exponent + * -- so no partial product, in particular neither ordinary determinant on its + * own, can overflow or underflow when the combined determinant is + * representable, whatever the ordering and magnitudes of the entries. + * Genuinely out-of-range determinants still saturate to (signed) zero or + * infinity. + * \warning The diagonal must have no zero entries (use the lemma symmetrically + * or a dense fallback for that case), and, as with \ref solve and \ref inverse, + * its reciprocals must be finite: a subnormal diagonal entry overflows + * \f$ D^{-1} \f$ -- and with it the capacitance entries -- to infinity. */ + Scalar determinant() const { + Scalar det(1); + Index exponent = 0; + for (Index i = 0; i < rows(); ++i) + det = internal::structured_balance(det * internal::structured_balance(m_d.coeff(i), exponent), exponent); + const Scalar capDet = internal::dplr_capacitance_impl<Rank_>::balancedCapacitanceDeterminant(*this, exponent); + det = internal::structured_balance(det * capDet, exponent); + return internal::structured_ldexp_clamped(det, exponent); + } + + /** \internal Writes the dense representation into \a dst. Invoked through + * \c dense = op; */ + template <typename Dest> + void evalTo(Dest& dst) const { + dst.noalias() = m_U * m_V.adjoint(); + dst.diagonal() += m_d; + } + + /** \internal Computes \c dst += (*this), see evalTo(). */ + template <typename Dest> + void addTo(Dest& dst) const { + dst.noalias() += m_U * m_V.adjoint(); + dst.diagonal() += m_d; + } + + /** \internal Computes \c dst -= (*this), see evalTo(). */ + template <typename Dest> + void subTo(Dest& dst) const { + dst.noalias() -= m_U * m_V.adjoint(); + dst.diagonal() -= m_d; + } + + /** \returns the product expression \c (*this) * \a v, evaluated at O(nk) + * operations without forming the matrix. The expression carries the default + * product tag, so assigning it behaves like any dense product: a temporary + * resolves aliasing between the destination and \a v, and \c .noalias() skips + * it. */ + template <typename Rhs> + Product<DiagonalPlusLowRank, Rhs> operator*(const MatrixBase<Rhs>& v) const { + EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic || + int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime), + INVALID_MATRIX_PRODUCT) + eigen_assert(v.rows() == cols() && "invalid product: dimensions do not match"); + return Product<DiagonalPlusLowRank, Rhs>(*this, v.derived()); + } + + /** \internal Computes \c dst += alpha * (*this) * rhs. \c ProductScalar is the + * promoted scalar of the product (complex when a real operator is applied to a + * complex right-hand side); the workspaces and the accumulation run in the + * promoted type. */ + template <typename Dest, typename Rhs, typename ProductScalar> + void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const { + eigen_assert(rhs.rows() == rows() && "invalid product: dimensions do not match"); + const Matrix<ProductScalar, Size_, Rhs::ColsAtCompileTime> r = rhs; + dst += alpha * (m_d.asDiagonal() * r); + if (correctionRank() > 0) { + const Matrix<ProductScalar, Rank_, Rhs::ColsAtCompileTime> t = m_V.adjoint() * r; + dst.noalias() += alpha * (m_U * t); + } + } + + private: + DiagonalVector m_d; + FactorType m_U; + FactorType m_V; +}; + +/** \ingroup StructuredMatrices_Module + * \returns a \ref DiagonalPlusLowRank operator \c diag(d) + \c U*V.adjoint(). + * The compile-time dimension and rank are deduced from the arguments. */ +template <typename DDerived, typename UDerived, typename VDerived> +DiagonalPlusLowRank<typename DDerived::Scalar, DDerived::SizeAtCompileTime, UDerived::ColsAtCompileTime> +makeDiagonalPlusLowRank(const MatrixBase<DDerived>& d, const MatrixBase<UDerived>& U, const MatrixBase<VDerived>& V) { + return DiagonalPlusLowRank<typename DDerived::Scalar, DDerived::SizeAtCompileTime, UDerived::ColsAtCompileTime>(d, U, + V); +} + +namespace internal { + +// Single product specialization covering every product tag; see the note in +// Circulant.h. +template <typename Scalar_, int Size_, int Rank_, typename Rhs, int ProductTag> +struct generic_product_impl<DiagonalPlusLowRank<Scalar_, Size_, Rank_>, Rhs, StructuredShape, DenseShape, ProductTag> + : structured_product_impl<DiagonalPlusLowRank<Scalar_, Size_, Rank_>, Rhs> {}; + +} // namespace internal + +} // namespace Eigen + +#endif // EIGEN_STRUCTURED_DIAGONAL_PLUS_LOW_RANK_H
diff --git a/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h b/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h index 51e7214..66bafba 100644 --- a/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h +++ b/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h
@@ -225,6 +225,11 @@ bool structured_exponent_bound_finite(const Xpr& x, int& e) { using ScalarTraits = NumTraits<typename Xpr::Scalar>; using RealScalar = typename ScalarTraits::Real; + // maxCoeff() asserts on an empty input, and a degenerate operand -- a rank-0 + // factor, a solve with no right-hand sides -- reaches here legitimately. + // An empty operand bounds nothing, so its exponent bound is 0. + e = 0; + if (x.size() == 0) return true; RealScalar m; if (ScalarTraits::IsComplex) // realView() reduces over both components in one pass, vectorized for
diff --git a/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt b/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt index 4c8aa3b..f87f4ea 100644 --- a/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt +++ b/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt
@@ -7,3 +7,4 @@ eigen_add_benchmark(bench_structured_levinson bench_structured_levinson.cpp) eigen_add_benchmark(bench_structured_transpose bench_structured_transpose.cpp) eigen_add_benchmark(bench_structured_kronecker bench_structured_kronecker.cpp) +eigen_add_benchmark(bench_structured_dplr bench_structured_dplr.cpp)
diff --git a/unsupported/benchmarks/StructuredMatrices/bench_structured_dplr.cpp b/unsupported/benchmarks/StructuredMatrices/bench_structured_dplr.cpp new file mode 100644 index 0000000..e22aa06 --- /dev/null +++ b/unsupported/benchmarks/StructuredMatrices/bench_structured_dplr.cpp
@@ -0,0 +1,72 @@ +// Benchmarks for the DiagonalPlusLowRank operator: the implicit O(nk) product +// and the O(nk^2 + k^3) Woodbury solve against their dense equivalents (the +// O(n^2) materialized product and the O(n^3) LU factor-and-solve). The solve +// benchmarks time factorization plus one right-hand side on both sides, the +// use case the Woodbury identity accelerates. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Core> +#include <Eigen/LU> +#include <unsupported/Eigen/StructuredMatrices> + +using namespace Eigen; + +typedef Matrix<double, Dynamic, 1> Vec; +typedef Matrix<double, Dynamic, Dynamic> Mat; + +// Well-conditioned operator: diagonal bounded away from zero and a contractive +// correction, so the capacitance matrix is safely invertible. +static DiagonalPlusLowRank<double> makeOperator(Index n, Index k) { + Vec d = Vec::Random(n).array() + 3.0; + Mat U = 0.5 * Mat::Random(n, k); + Mat V = 0.5 * Mat::Random(n, k); + return DiagonalPlusLowRank<double>(d, U, V); +} + +// --- Matrix-vector product: implicit O(nk) vs. materialized dense O(n^2) --- +static void BM_DplrProduct(benchmark::State& state) { + const Index n = state.range(0), k = state.range(1); + DiagonalPlusLowRank<double> A = makeOperator(n, k); + Vec x = Vec::Random(n), y(n); + for (auto _ : state) { + y.noalias() = A * x; + benchmark::DoNotOptimize(y.data()); + } +} +BENCHMARK(BM_DplrProduct)->ArgsProduct({{256, 1024, 4096}, {4, 32}}); + +static void BM_DenseProduct(benchmark::State& state) { + const Index n = state.range(0), k = state.range(1); + Mat dense = makeOperator(n, k); + Vec x = Vec::Random(n), y(n); + for (auto _ : state) { + y.noalias() = dense * x; + benchmark::DoNotOptimize(y.data()); + } +} +BENCHMARK(BM_DenseProduct)->ArgsProduct({{256, 1024, 4096}, {4, 32}}); + +// --- Solve: Woodbury (k x k capacitance LU) vs. dense n x n LU --- +static void BM_DplrSolve(benchmark::State& state) { + const Index n = state.range(0), k = state.range(1); + DiagonalPlusLowRank<double> A = makeOperator(n, k); + Vec b = Vec::Random(n), x(n); + for (auto _ : state) { + x = A.solve(b); // factors the k x k capacitance matrix each iteration + benchmark::DoNotOptimize(x.data()); + } +} +BENCHMARK(BM_DplrSolve)->ArgsProduct({{256, 1024, 4096}, {4, 32}}); + +static void BM_DenseLuSolve(benchmark::State& state) { + const Index n = state.range(0), k = state.range(1); + Mat dense = makeOperator(n, k); + Vec b = Vec::Random(n), x(n); + for (auto _ : state) { + x = dense.partialPivLu().solve(b); // factors the n x n matrix each iteration + benchmark::DoNotOptimize(x.data()); + } +} +BENCHMARK(BM_DenseLuSolve)->ArgsProduct({{256, 1024, 4096}, {4, 32}});
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 61b83a6..a84ab5c 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt
@@ -128,6 +128,7 @@ ei_add_test(polynomialsolver) ei_add_test(polynomialutils) ei_add_test(splines) +ei_add_test(structured_dplr) ei_add_test(levenberg_marquardt) ei_add_test(kronecker_product) ei_add_test(structured_matrices)
diff --git a/unsupported/test/structured_dplr.cpp b/unsupported/test/structured_dplr.cpp new file mode 100644 index 0000000..2c915f1 --- /dev/null +++ b/unsupported/test/structured_dplr.cpp
@@ -0,0 +1,696 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "main.h" + +#include <unsupported/Eigen/StructuredMatrices> + +using namespace Eigen; + +// Well-conditioned random operator: diagonal bounded away from zero and a +// contractive correction, so the capacitance matrix is safely invertible. +template <typename Scalar> +void random_dplr(Index n, Index k, Matrix<Scalar, Dynamic, 1>& d, Matrix<Scalar, Dynamic, Dynamic>& U, + Matrix<Scalar, Dynamic, Dynamic>& V) { + typedef typename NumTraits<Scalar>::Real RealScalar; + d = Matrix<Scalar, Dynamic, 1>::Random(n); + d.array() += Scalar(RealScalar(3)); + U = Scalar(RealScalar(0.5)) * Matrix<Scalar, Dynamic, Dynamic>::Random(n, k); + V = Scalar(RealScalar(0.5)) * Matrix<Scalar, Dynamic, Dynamic>::Random(n, k); +} + +// Reference dense matrix built independently, entry by entry. +template <typename Scalar> +Matrix<Scalar, Dynamic, Dynamic> reference_dplr(const Matrix<Scalar, Dynamic, 1>& d, + const Matrix<Scalar, Dynamic, Dynamic>& U, + const Matrix<Scalar, Dynamic, Dynamic>& V) { + const Index n = d.size(); + Matrix<Scalar, Dynamic, Dynamic> dense(n, n); + for (Index j = 0; j < n; ++j) + for (Index i = 0; i < n; ++i) { + Scalar s = (i == j) ? d[i] : Scalar(0); + for (Index t = 0; t < U.cols(); ++t) s += U(i, t) * numext::conj(V(j, t)); + dense(i, j) = s; + } + return dense; +} + +template <typename Scalar> +void test_dplr_product(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + VERIFY_IS_EQUAL(A.correctionRank(), k); + + Mat Ad = A; + VERIFY_IS_APPROX(Ad, dense); + for (Index t = 0; t < 5; ++t) { + Index i = internal::random<Index>(0, n - 1), j = internal::random<Index>(0, n - 1); + VERIFY_IS_APPROX(A.coeff(i, j), dense(i, j)); + } + + Vec x = Vec::Random(n); + VERIFY_IS_APPROX((A * x).eval(), (dense * x).eval()); + + Mat X = Mat::Random(n, 3); + VERIFY_IS_APPROX((A * X).eval(), (dense * X).eval()); + + // Accumulation form exercised by the iterative solvers. + Vec y = Vec::Random(n); + Vec y0 = y; + y.noalias() += A * x; + VERIFY_IS_APPROX(y, (y0 + dense * x).eval()); + + // Dimension mismatches are caught when the product expression is formed. + VERIFY_RAISES_ASSERT(A * Vec(Vec::Random(n + 1))); +} + +template <typename Scalar> +void test_dplr_transpose(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + Mat Td = A.transpose(); + VERIFY_IS_APPROX(Td, Mat(dense.transpose())); + Mat Ad = A.adjoint(); + VERIFY_IS_APPROX(Ad, Mat(dense.adjoint())); + Mat Kd = A.conjugate(); + VERIFY_IS_APPROX(Kd, Mat(dense.conjugate())); + + Vec x = Vec::Random(n); + VERIFY_IS_APPROX((A.transpose() * x).eval(), (dense.transpose() * x).eval()); + VERIFY_IS_APPROX((A.adjoint() * x).eval(), (dense.adjoint() * x).eval()); +} + +template <typename Scalar> +void test_dplr_solve(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + Vec b = Vec::Random(n); + Vec x = A.solve(b); + VERIFY_IS_APPROX((dense * x).eval(), b); + VERIFY_IS_APPROX(x, dense.partialPivLu().solve(b).eval()); + + Mat B = Mat::Random(n, 3); + Mat X = A.solve(B); + VERIFY_IS_APPROX(X, dense.partialPivLu().solve(B).eval()); + + Mat empty_rhs(n, 0); + Mat empty_solution = A.solve(empty_rhs); + VERIFY_IS_EQUAL(empty_solution.rows(), n); + VERIFY_IS_EQUAL(empty_solution.cols(), 0); +} + +template <typename Scalar> +void test_dplr_inverse_determinant(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + // The inverse is itself diagonal-plus-low-rank; materialize and check. + DiagonalPlusLowRank<Scalar> Ainv = A.inverse(); + VERIFY_IS_EQUAL(Ainv.correctionRank(), k); + Mat invd = Ainv; + VERIFY_IS_APPROX((invd * dense).eval(), Mat(Mat::Identity(n, n))); + + // Inverting twice returns to the original operator. + Mat back = Ainv.inverse(); + VERIFY_IS_APPROX(back, dense); + + VERIFY_IS_APPROX(A.determinant(), dense.determinant()); +} + +template <typename Scalar> +void test_dplr_matrix_free_gmres(Index n, Index k) { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + Vec b = Vec::Random(n); + GMRES<DiagonalPlusLowRank<Scalar>, IdentityPreconditioner> gmres; + gmres.setTolerance(RealScalar(1e-12)); + gmres.compute(A); + Vec x = gmres.solve(b); + VERIFY(gmres.info() == Success); + const RealScalar tol = RealScalar(5e6) * NumTraits<RealScalar>::epsilon(); // ~1e-9 in double + VERIFY((dense * x - b).norm() <= tol * b.norm()); +} + +// The products carry the default product tag, so an assignment whose right-hand +// side aliases the destination materializes the product into a temporary exactly +// like a dense product would -- including the aliasing that the former +// is_same_dense check could not see: overlapping views and right-hand-side +// expressions referencing the destination. +template <typename Scalar> +void test_dplr_aliased_product(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + // Pin the routing: the structured products must dispatch through the shared + // structured_product_impl. + STATIC_CHECK( + (std::is_base_of<internal::structured_product_impl<DiagonalPlusLowRank<Scalar>, Vec>, + internal::generic_product_impl<DiagonalPlusLowRank<Scalar>, Vec, internal::StructuredShape, + DenseShape, GemvProduct>>::value)); + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + Vec x = Vec::Random(n); + Vec y = x; + y = A * y; + VERIFY_IS_APPROX(y, (dense * x).eval()); + + y = x; + y += A * y; + VERIFY_IS_APPROX(y, (x + dense * x).eval()); + + // Right-hand-side expression referencing the destination. + y = x; + y = A * (y + Vec::Ones(n)); + VERIFY_IS_APPROX(y, (dense * (x + Vec::Ones(n))).eval()); + + // Overlapping views of a shared buffer, in both directions. + Vec w = Vec::Random(n + 1); + Vec w0 = w; + w.tail(n) = A * w.head(n); + VERIFY_IS_APPROX(w.tail(n).eval(), (dense * w0.head(n)).eval()); + w = w0; + w.head(n) = A * w.tail(n); + VERIFY_IS_APPROX(w.head(n).eval(), (dense * w0.tail(n)).eval()); +} + +// Product expressions must nest the structured operand by value: a +// delayed-evaluated expression has to outlive the temporary operator it was +// built from. The static check pins the value nesting; the behavioral check +// would read freed memory if the product held a reference instead. +template <typename Scalar> +void test_dplr_delayed_product(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + STATIC_CHECK(!std::is_reference<typename internal::ref_selector<DiagonalPlusLowRank<Scalar>>::type>::value); + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + Mat dense = reference_dplr<Scalar>(d, U, V); + + Vec x = Vec::Random(n); + auto expr = A.adjoint() * x; // the adjoint temporary dies with the full expression + Vec scribble = Vec::Random(2 * n); // reuses the temporary's freed heap storage + Vec y = expr; + VERIFY_IS_APPROX(y, (dense.adjoint() * x).eval()); + VERIFY_IS_EQUAL(scribble.size(), 2 * n); // keep the scribble alive across the evaluation +} + +// Mixed-scalar products: a real operator applied to a complex right-hand side +// (and a complex operator applied to a real one) promotes to the complex product +// scalar, so alpha, the workspaces and the accumulation must run in the promoted +// type rather than the operator scalar. +template <typename RealScalar> +void test_dplr_mixed_scalar(Index n, Index k) { + typedef std::complex<RealScalar> Complex; + typedef Matrix<RealScalar, Dynamic, 1> RVec; + typedef Matrix<RealScalar, Dynamic, Dynamic> RMat; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + RVec d; + RMat U, V; + random_dplr<RealScalar>(n, k, d, U, V); + DiagonalPlusLowRank<RealScalar> A(d, U, V); + CMat dense = reference_dplr<RealScalar>(d, U, V).template cast<Complex>(); + + CVec x = CVec::Random(n); + CVec y = A * x; + VERIFY_IS_APPROX(y, (dense * x).eval()); + + CVec y0 = CVec::Random(n); + y = y0; + y.noalias() += A * x; + VERIFY_IS_APPROX(y, (y0 + dense * x).eval()); + + CVec dc; + CMat Uc, Vc; + random_dplr<Complex>(n, k, dc, Uc, Vc); + DiagonalPlusLowRank<Complex> Ac(dc, Uc, Vc); + CMat denseC = reference_dplr<Complex>(dc, Uc, Vc); + RVec xr = RVec::Random(n); + CVec z = Ac * xr; + VERIFY_IS_APPROX(z, (denseC * xr).eval()); +} + +// Wide-dynamic-range determinants: the diagonal product is accumulated in the +// balanced form m * 2^e, so partial products neither overflow nor underflow +// when the determinant itself is representable, whatever the ordering of large +// and small diagonal entries. Genuinely non-representable determinants must +// still saturate to infinity / zero. +void test_dplr_determinant_scaled() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + typedef std::complex<double> Complex; + typedef Matrix<Complex, Dynamic, 1> CVec; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + + { + // The reviewer's magnitude regime: det(D) = 1e200 * 1e200 * 1e-100 = 1e300 is + // finite, but the plain running product hits 1e400 partway. k = 0 at runtime. + Vec d(3); + d << 1e200, 1e200, 1e-100; + DiagonalPlusLowRank<double> A(d, Mat::Zero(3, 0), Mat::Zero(3, 0)); + const double det = A.determinant(); + VERIFY((numext::isfinite)(det)); + VERIFY_IS_APPROX(det, 1e300); + } + { + // Underflow-side analogue: the partial product flushes to zero at 1e-400 + // while det(D) = 1e-300 is representable. + Vec d(3); + d << 1e-200, 1e-200, 1e100; + DiagonalPlusLowRank<double> A(d, Mat::Zero(3, 0), Mat::Zero(3, 0)); + const double det = A.determinant(); + VERIFY(det != 0.0); + VERIFY_IS_APPROX(det, 1e-300); + } + { + // The capacitance determinant enters the same balanced accumulation: with + // 450 entries of 10 first, the running diagonal product tops out at 1e450 + // while det(D) = 10^450 * 10^-150 = 1e300; the rank-1 correction + // U = d[0] e_0, V = e_0 gives the capacitance 1 + d[0]/d[0] = 2, so + // det(A) = 2e300 by the determinant lemma. + const Index n = 600; + Vec d(n); + d.head(450).setConstant(10.0); + d.tail(150).setConstant(0.1); + Mat U = Mat::Zero(n, 1), V = Mat::Zero(n, 1); + U(0, 0) = d[0]; + V(0, 0) = 1.0; + DiagonalPlusLowRank<double> A(d, U, V); + const double det = A.determinant(); + // One rounding multiply per balanced factor, and the decimal constants 0.1 + // and their 150-fold product are themselves inexact: a few n eps total. + const double det_tol = 4.0 * double(n) * NumTraits<double>::epsilon(); + VERIFY((numext::isfinite)(det)); + VERIFY(numext::abs(det / 2e300 - 1.0) <= det_tol); + } + { + // Complex scalars balance on max(|re|, |im|): det = i * 1e200 * 1e100 = 1e300 i. + CVec d(3); + d << Complex(0, 1e200), Complex(1e200, 0), Complex(1e-100, 0); + DiagonalPlusLowRank<Complex> A(d, CMat::Zero(3, 0), CMat::Zero(3, 0)); + VERIFY_IS_APPROX(A.determinant(), Complex(0, 1e300)); + } + { + // A genuinely overflowing determinant still reports infinity... + Vec d(2); + d << 1e200, 1e200; + DiagonalPlusLowRank<double> A(d, Mat::Zero(2, 0), Mat::Zero(2, 0)); + VERIFY((numext::isinf)(A.determinant())); + } + { + // ... and a genuinely vanishing one an exact zero. + Vec d(2); + d << 1e-200, 1e-200; + DiagonalPlusLowRank<double> A(d, Mat::Zero(2, 0), Mat::Zero(2, 0)); + VERIFY_IS_EQUAL(A.determinant(), 0.0); + } +} + +// The two factors of the determinant lemma can sit at opposite ends of the exponent range even when +// the determinant is ~1, since D^{-1} scales the capacitance entries by the reciprocal diagonal. Both +// det(D) and det(capacitance) must therefore accumulate in the balanced mantissa * 2^e form, the +// latter from the LU pivots rather than the plain pivot product .determinant() would form. +void test_dplr_determinant_capacitance() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef Matrix<double, Dynamic, Dynamic> Mat; + typedef std::complex<double> Complex; + typedef Matrix<Complex, Dynamic, Dynamic> CMat; + const double eps = NumTraits<double>::epsilon(); + + { + // The reviewer's reproducer: D + I I^H = diag(1 + 1e-200) is essentially the + // identity, but det(D) = 1e-400 underflows and det(capacitance) ~ 1e400 + // overflows when either is formed on its own. + Vec d(2); + d << 1e-200, 1e-200; + DiagonalPlusLowRank<double> A(d, Mat::Identity(2, 2), Mat::Identity(2, 2)); + const double det = A.determinant(); + VERIFY(numext::abs(det - 1.0) <= 16.0 * eps); + } + { + // Mirrored, in exact powers of two: d = 2^60 and U = -(2^60 - 2^8) I with + // V = I represent diag(2^8), so det(A) = 2^176 exactly. Every intermediate + // is a power-of-two computation (the capacitance is exactly 2^-52 I by + // Sterbenz subtraction), yet det(D) = 2^1342-ish overflows and the plain + // capacitance pivot product 2^-1144 underflows. + const Index n = 22; + Vec d = Vec::Constant(n, std::ldexp(1.0, 60)); + Mat U = (std::ldexp(1.0, 8) - std::ldexp(1.0, 60)) * Mat::Identity(n, n); + Mat V = Mat::Identity(n, n); + DiagonalPlusLowRank<double> A(d, U, V); + VERIFY_IS_EQUAL(A.determinant(), std::ldexp(1.0, 176)); + } + { + // Complex analogue of the reproducer: the mantissa stays complex and is + // balanced on max(|re|, |im|). D + I I^H = diag(1 + 1e-200 i), det ~ 1. + Matrix<Complex, Dynamic, 1> d(2); + d << Complex(0, 1e-200), Complex(0, 1e-200); + DiagonalPlusLowRank<Complex> A(d, CMat::Identity(2, 2), CMat::Identity(2, 2)); + VERIFY_IS_APPROX(A.determinant(), Complex(1, 0)); + } + { + // Sign correctness through the capacitance path at extreme scales: + // D + (-I) I^H = diag(1e-200 - 1) ~ -I_3, det ~ -1; the balanced pivot + // product of the capacitance ~ -1e200 I_3 carries the sign. + Vec d = Vec::Constant(3, 1e-200); + DiagonalPlusLowRank<double> A(d, (-Mat::Identity(3, 3)).eval(), Mat::Identity(3, 3)); + VERIFY_IS_APPROX(A.determinant(), -1.0); + } + { + // Sign correctness across mixed extreme diagonal scales (rank 0). + Vec d(4); + d << -1e200, 1e-200, 1e200, 1e-200; + DiagonalPlusLowRank<double> A(d, Mat::Zero(4, 0), Mat::Zero(4, 0)); + VERIFY_IS_APPROX(A.determinant(), -1.0); + } + { + // Genuine overflow saturates to infinity of the correct sign... + Vec d(3); + d << -1e200, 1e200, 1e200; + DiagonalPlusLowRank<double> A(d, Mat::Zero(3, 0), Mat::Zero(3, 0)); + const double det = A.determinant(); + VERIFY((numext::isinf)(det) && det < 0.0); + } + { + // ... and genuine underflow to a zero of the correct sign. + Vec d(3); + d << -1e-200, 1e-200, 1e-200; + DiagonalPlusLowRank<double> A(d, Mat::Zero(3, 0), Mat::Zero(3, 0)); + const double det = A.determinant(); + VERIFY_IS_EQUAL(det, 0.0); + VERIFY(std::signbit(det)); + } +} + +// The capacitance triple product V^H D^{-1} U can overflow in its plain association while the +// capacitance itself is representable: D^{-1} pairs the reciprocal of a tiny diagonal with a huge +// factor before the other, tiny factor pulls the product back down. The Woodbury kernels form these +// products from exactly rescaled factors so no spurious Inf/NaN reaches determinant(), solve() or +// inverse(). +void test_dplr_capacitance_overflow() { + typedef Matrix<double, Dynamic, 1> Vec; + typedef std::complex<double> Complex; + typedef Matrix<Complex, Dynamic, 1> CVec; + const double eps = NumTraits<double>::epsilon(); + + { + // The reviewer's case: D + UV^H = 1e-200 + 1e-200 * 1e200 ~ 1 and the + // capacitance 1 + 1e200 is representable, but V^H D^{-1} = 1e400 is not; + // determinant() used to return Inf and solve() NaN. + Vec d(1), u(1), v(1); + d << 1e-200; + u << 1e-200; + v << 1e200; + DiagonalPlusLowRank<double> A(d, u, v); + VERIFY((numext::isfinite)(A.capacitance()(0, 0))); + VERIFY_IS_APPROX(A.capacitance()(0, 0), 1e200); + VERIFY(numext::abs(A.determinant() - 1.0) <= 8.0 * eps); + // The Woodbury splitting routes this solve through D^{-1} b ~ 1e200 and + // cancels it back down to x ~ 1, so the amplified digits are lost to any + // association (see the solve() documentation); the fix is that the solve + // stays finite instead of turning into NaN through an infinite capacitance. + Vec b(1); + b << 1.0; + VERIFY(A.solve(b).allFinite()); + + // D^{-1} b is 1e400 here, though the operator is essentially the identity and the solution 1e200 + // is representable; without the power-of-two rescaling the first term is Inf and the correction + // turns it into NaN, which is what this pins. It cannot pin finiteness: A ~ 1 makes the correction + // cancel D^{-1} b over ~200 digits, and whether that leaves 0 or an ulp the undone scaling lifts + // back to Inf is a codegen detail (x86-64 gives the first, aarch64 the second). + Vec big_b(1); + big_b << 1e200; + VERIFY(!A.solve(big_b).array().isNaN().any()); + } + { + // The same overflow pattern in exact powers of two, where the Woodbury + // arithmetic is exact end to end: d = 2^-522, U = 2^-1025, V = 2^503 + // represent A = [2^-521] exactly; V^H D^{-1} = 2^1025 still overflows, and + // the capacitance is exactly 2. Determinant and solve must match the dense + // results for the represented matrix exactly. + Vec d(1), u(1), v(1); + d << std::ldexp(1.0, -522); + u << std::ldexp(1.0, -1025); + v << std::ldexp(1.0, 503); + DiagonalPlusLowRank<double> A(d, u, v); + VERIFY_IS_EQUAL(A.capacitance()(0, 0), 2.0); + VERIFY_IS_EQUAL(A.determinant(), std::ldexp(1.0, -521)); + Vec b(1); + b << 1.0; + VERIFY_IS_EQUAL(A.solve(b)(0), std::ldexp(1.0, 521)); + } + { + // Symmetric mirror (U huge, V tiny): the plain capacitance association is + // safe here, but the inverse's D^{-1} U product overflows spuriously while + // the inverse factor -1e400 / (1 + 1e200) ~ -1e200 is representable. + Vec d(1), u(1), v(1); + d << 1e-200; + u << 1e200; + v << 1e-200; + DiagonalPlusLowRank<double> A(d, u, v); + VERIFY(numext::abs(A.determinant() - 1.0) <= 8.0 * eps); + DiagonalPlusLowRank<double> Ainv = A.inverse(); + VERIFY(Ainv.factorU().allFinite()); + VERIFY_IS_APPROX(Ainv.factorU()(0, 0), -1e200); + Vec b(1); + b << 1.0; + VERIFY(A.solve(b).allFinite()); + } + { + // Exact-power-of-two mirror: d = 2^-522, U = 2^502, V = 2^-1025 represent + // A = [1.5 * 2^-522]; D^{-1} U = 2^1024 overflows, yet the inverse factor + // -2^1024 / 1.5 is representable and the solve is exact to rounding. + Vec d(1), u(1), v(1); + d << std::ldexp(1.0, -522); + u << std::ldexp(1.0, 502); + v << std::ldexp(1.0, -1025); + DiagonalPlusLowRank<double> A(d, u, v); + VERIFY_IS_EQUAL(A.capacitance()(0, 0), 1.5); + VERIFY_IS_EQUAL(A.determinant(), 1.5 * std::ldexp(1.0, -522)); + DiagonalPlusLowRank<double> Ainv = A.inverse(); + VERIFY_IS_APPROX(Ainv.factorU()(0, 0), -std::ldexp(1.0 / 1.5, 1024)); + Vec b(1); + b << 3.0; + VERIFY_IS_EQUAL(A.solve(b)(0), std::ldexp(1.0, 523)); + } + { + // Complex analogue of the reviewer's case: the operator is essentially + // [i], the capacitance 1 + 1e200 is representable, and the balancing and + // rescaling act on component magnitudes. + CVec d(1), u(1), v(1); + d << Complex(0, 1e-200); + u << Complex(1e-200, 0); + v << Complex(0, -1e200); + DiagonalPlusLowRank<Complex> A(d, u, v); + VERIFY((numext::isfinite)(numext::real(A.capacitance()(0, 0)))); + VERIFY_IS_APPROX(A.capacitance()(0, 0), Complex(1e200, 0)); + VERIFY_IS_APPROX(A.determinant(), Complex(0, 1)); + } +} + +// Moderate data must keep the plain product association: the normalized +// Woodbury kernels only engage when a conservative exponent bound says a +// product could overflow, so for well-scaled operators the capacitance and the +// solve are bit-identical to the unnormalized evaluation. +template <typename Scalar> +void test_dplr_capacitance_fast_path(Index n, Index k) { + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + + Vec d; + Mat U, V; + random_dplr<Scalar>(n, k, d, U, V); + DiagonalPlusLowRank<Scalar> A(d, U, V); + + Mat capRef = Mat::Identity(k, k); + capRef.noalias() += V.adjoint() * d.cwiseInverse().asDiagonal() * U; + VERIFY_IS_CWISE_EQUAL(A.capacitance(), capRef); + + Vec b = Vec::Random(n); + PartialPivLU<Mat> capLU(capRef); + Vec xRef = d.cwiseInverse().asDiagonal() * b; + xRef.noalias() -= d.cwiseInverse().asDiagonal() * (U * capLU.solve(V.adjoint() * xRef)); + VERIFY_IS_CWISE_EQUAL(A.solve(b), xRef); +} + +// A fixed correction rank of zero is a compile-level regression: the operator +// must instantiate -- the capacitance solve is compile-time dispatched away, so +// PartialPivLU<Matrix<Scalar, 0, 0>> is never formed -- and behave as the plain +// diagonal it is. +template <typename Scalar, int N> +void test_dplr_fixed_rank0() { + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, N, 1> VecN; + typedef Matrix<Scalar, N, 0> FacN0; + typedef Matrix<Scalar, N, N> MatN; + + VecN d = VecN::Random(); + d.array() += Scalar(RealScalar(3)); + FacN0 U, V; + DiagonalPlusLowRank<Scalar, N, 0> A(d, U, V); + VERIFY_IS_EQUAL(A.correctionRank(), Index(0)); + + MatN dense = MatN::Zero(); + dense.diagonal() = d; + MatN Ad = A; + VERIFY_IS_APPROX(Ad, dense); + + VecN x = VecN::Random(); + VecN y = A * x; + VERIFY_IS_APPROX(y, (d.asDiagonal() * x).eval()); + + VecN b = VecN::Random(); + VecN xs = A.solve(b); + VERIFY_IS_APPROX(xs, (b.array() / d.array()).matrix().eval()); + + VERIFY_IS_APPROX(A.determinant(), d.prod()); + + DiagonalPlusLowRank<Scalar, N, 0> Ainv = A.inverse(); + MatN invd = Ainv; + VERIFY_IS_APPROX((invd * dense).eval(), MatN(MatN::Identity())); + + // Dynamic size with the rank still fixed at zero takes the same dispatch. + Matrix<Scalar, Dynamic, 1> dd = d; + Matrix<Scalar, Dynamic, 0> Ud(N, 0), Vd(N, 0); + DiagonalPlusLowRank<Scalar, Dynamic, 0> B(dd, Ud, Vd); + Matrix<Scalar, Dynamic, 1> yb = B * Matrix<Scalar, Dynamic, 1>(x); + VERIFY_IS_APPROX(yb, (d.asDiagonal() * x).eval()); + VERIFY_IS_APPROX(B.solve(Matrix<Scalar, Dynamic, 1>(b)).eval(), (b.array() / d.array()).matrix().eval()); + VERIFY_IS_APPROX(B.determinant(), d.prod()); +} + +template <typename Scalar, int N, int K> +void test_dplr_fixed() { + typedef Matrix<Scalar, N, 1> VecN; + typedef Matrix<Scalar, N, K> FacNK; + typedef Matrix<Scalar, Dynamic, 1> Vec; + typedef Matrix<Scalar, Dynamic, Dynamic> Mat; + typedef typename NumTraits<Scalar>::Real RealScalar; + + VecN d = VecN::Random(); + d.array() += Scalar(RealScalar(3)); + FacNK U = FacNK::Random(), V = FacNK::Random(); + DiagonalPlusLowRank<Scalar, N, K> A(d, U, V); + STATIC_CHECK((DiagonalPlusLowRank<Scalar, N, K>::RowsAtCompileTime == N)); + STATIC_CHECK((internal::remove_all_t<decltype(makeDiagonalPlusLowRank(d, U, V))>::RowsAtCompileTime == N)); + + Mat dense = reference_dplr<Scalar>(Vec(d), Mat(U), Mat(V)); + Matrix<Scalar, N, N> Adense = A; + VERIFY_IS_APPROX(Mat(Adense), dense); + + VecN x = VecN::Random(); + VecN y = A * x; + VERIFY_IS_APPROX(y, (dense * x).eval()); + + VecN b = VecN::Random(); + VecN xs = A.solve(b); + VERIFY_IS_APPROX((dense * xs).eval(), b); + + // A dynamic right-hand side against the fixed-size operator: the compile-time + // product-dimension check must accept the Dynamic/fixed mix. + Vec xd = Vec(x); + Vec yd = A * xd; + VERIFY_IS_APPROX(yd, (dense * xd).eval()); +} + +EIGEN_DECLARE_TEST(structured_dplr) { + for (int i = 0; i < g_repeat; ++i) { + // Products, dense assignment, coefficient access; k = 0 is a plain diagonal + // and k = n a full-rank correction. + CALL_SUBTEST_1((test_dplr_product<double>(1, 0))); + CALL_SUBTEST_1((test_dplr_product<double>(10, 0))); + CALL_SUBTEST_1((test_dplr_product<double>(10, 1))); + CALL_SUBTEST_1((test_dplr_product<double>(20, 3))); + CALL_SUBTEST_1((test_dplr_product<double>(8, 8))); + CALL_SUBTEST_1((test_dplr_product<float>(12, 2))); + CALL_SUBTEST_1((test_dplr_product<std::complex<double>>(9, 2))); + CALL_SUBTEST_1((test_dplr_product<std::complex<float>>(7, 3))); + CALL_SUBTEST_1((test_dplr_transpose<double>(10, 2))); + CALL_SUBTEST_1((test_dplr_transpose<std::complex<double>>(8, 3))); + + // Woodbury solves against the dense LU. + CALL_SUBTEST_2((test_dplr_solve<double>(1, 0))); + CALL_SUBTEST_2((test_dplr_solve<double>(10, 0))); + CALL_SUBTEST_2((test_dplr_solve<double>(10, 1))); + CALL_SUBTEST_2((test_dplr_solve<double>(30, 4))); + CALL_SUBTEST_2((test_dplr_solve<double>(8, 8))); + CALL_SUBTEST_2((test_dplr_solve<float>(12, 2))); + CALL_SUBTEST_2((test_dplr_solve<std::complex<double>>(14, 3))); + CALL_SUBTEST_2((test_dplr_solve<std::complex<float>>(8, 2))); + + // Inverse closure, determinant lemma, matrix-free GMRES, fixed sizes. + CALL_SUBTEST_3((test_dplr_inverse_determinant<double>(10, 2))); + CALL_SUBTEST_3((test_dplr_inverse_determinant<double>(6, 0))); + CALL_SUBTEST_3((test_dplr_inverse_determinant<std::complex<double>>(8, 3))); + CALL_SUBTEST_3((test_dplr_matrix_free_gmres<double>(24, 3))); + CALL_SUBTEST_3((test_dplr_fixed<double, 6, 2>())); + CALL_SUBTEST_3((test_dplr_fixed<std::complex<float>, 5, 1>())); + + // Review regressions: aliased and value-nested (owning) delayed products, + // mixed-scalar products, the balanced determinant accumulation, and the + // fixed rank-0 instantiation. + CALL_SUBTEST_4((test_dplr_aliased_product<double>(11, 3))); + CALL_SUBTEST_4((test_dplr_aliased_product<std::complex<double>>(9, 2))); + CALL_SUBTEST_4((test_dplr_delayed_product<double>(13, 2))); + CALL_SUBTEST_4((test_dplr_delayed_product<std::complex<double>>(10, 3))); + CALL_SUBTEST_4((test_dplr_mixed_scalar<double>(12, 3))); + CALL_SUBTEST_4((test_dplr_mixed_scalar<float>(9, 2))); + CALL_SUBTEST_4(test_dplr_determinant_scaled()); + CALL_SUBTEST_4(test_dplr_determinant_capacitance()); + CALL_SUBTEST_4(test_dplr_capacitance_overflow()); + CALL_SUBTEST_4((test_dplr_capacitance_fast_path<double>(11, 3))); + CALL_SUBTEST_4((test_dplr_capacitance_fast_path<std::complex<double>>(9, 2))); + CALL_SUBTEST_4((test_dplr_fixed_rank0<double, 5>())); + CALL_SUBTEST_4((test_dplr_fixed_rank0<std::complex<float>, 4>())); + } +}