From 29277be12bac96caf7396022ac4a7beee815208b Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Wed, 8 Jul 2026 07:26:11 -0500 Subject: [PATCH 1/3] ENH: Add Eigen-backed itk::Math::SolveSymmetric and InverseSymmetric Add itk::Math::SolveSymmetric and itk::Math::InverseSymmetric for symmetric (indefinite) linear systems and symmetric-matrix inversion, backed by Eigen's LDLT (Bunch-Kaufman) factorization. Overloads accept itk::Array2D/itk::Array and map the contiguous block through Eigen::Map with no copy; a matrix-RHS overload solves A X = B with a single factorization. ITK_MATH_HAS_SOLVE_SYMMETRIC advertises the capability so consumers can feature-test the header. itkMathLDLTGTest checks SPD and indefinite systems, the matrix-RHS solve and the symmetric inverse against a vnl reference. --- Modules/Core/Common/include/itkMathLDLT.h | 274 ++++++++++++++++++ Modules/Core/Common/test/CMakeLists.txt | 1 + Modules/Core/Common/test/itkMathLDLTGTest.cxx | 242 ++++++++++++++++ 3 files changed, 517 insertions(+) create mode 100644 Modules/Core/Common/include/itkMathLDLT.h create mode 100644 Modules/Core/Common/test/itkMathLDLTGTest.cxx diff --git a/Modules/Core/Common/include/itkMathLDLT.h b/Modules/Core/Common/include/itkMathLDLT.h new file mode 100644 index 00000000000..2d37372b348 --- /dev/null +++ b/Modules/Core/Common/include/itkMathLDLT.h @@ -0,0 +1,274 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMathLDLT_h +#define itkMathLDLT_h + +#include "itkMacro.h" +#include "itkArray.h" +#include "itkArray2D.h" +#include "itkMatrix.h" +#include "itkVector.h" +#include "vnl/vnl_matrix.h" +#include "vnl/vnl_matrix_fixed.h" +#include "vnl/vnl_vector.h" +#include "vnl/vnl_vector_fixed.h" + +#include "itk_eigen.h" +#include ITK_EIGEN(Dense) + +#include + +/** Capability macro. Downstream code selects the Eigen-backed symmetric solve + * with, e.g.: + * \code + * #if __has_include() + * # include + * #endif + * #ifdef ITK_MATH_HAS_SOLVE_SYMMETRIC + * x = itk::Math::SolveSymmetric(A, b); // itk::Array2D, itk::Array + * #else + * x = vnl_cholesky(A.as_ref()).solve(b); // legacy vnl fallback + * #endif + * \endcode + */ +#define ITK_MATH_HAS_SOLVE_SYMMETRIC 1 + +namespace itk +{ +namespace Math +{ +namespace detail +{ +// Solve the symmetric system A x = b via Eigen's robust LDLT (Bunch-Kaufman +// pivoting). A is assumed symmetric; column-major mapping of ITK's row-major +// storage is exact for a symmetric matrix (A == A^T). Throws on a non-finite +// input. +template +void +SolveSymmetricLDLTEigen(const TReal * aData, const TReal * bData, unsigned int n, TReal * xData) +{ + using ColMajor = Eigen::Matrix; + using Vector = Eigen::Matrix; + + const Eigen::Map aMap(aData, n, n); + const Eigen::Map bMap(bData, n); + + const Eigen::LDLT ldlt(aMap); + if (ldlt.info() != Eigen::Success) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric failed; input is likely non-finite (NaN/Inf)."); + } + Eigen::Map(xData, n) = ldlt.solve(bMap); +} + +// Solve A X = B for symmetric A (n x n) and B (n x m) with a SINGLE LDLT +// factorization (all columns solved together). A is symmetric so its col-major +// map of ITK's row-major storage is exact; B and X are mapped row-major. Throws +// on a non-finite input. +template +void +SolveSymmetricMatrixLDLTEigen(const TReal * aData, const TReal * bData, unsigned int n, unsigned int m, TReal * xData) +{ + using ColMajor = Eigen::Matrix; + using RowMajor = Eigen::Matrix; + + const Eigen::Map aMap(aData, n, n); + const Eigen::Map bMap(bData, n, m); + + const Eigen::LDLT ldlt(aMap); + if (ldlt.info() != Eigen::Success) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric failed; input is likely non-finite (NaN/Inf)."); + } + Eigen::Map(xData, n, m) = ldlt.solve(bMap); +} + +// Invert symmetric A via LDLT (solves A X = I). Eigen's LDLT::solve silently +// pseudo-inverts zero pivots, so singularity is rejected here from vectorD(). +template +void +InverseSymmetricLDLTEigen(const TReal * aData, unsigned int n, TReal * invData) +{ + using ColMajor = Eigen::Matrix; + using RowMajor = Eigen::Matrix; + + const Eigen::Map aMap(aData, n, n); + + const Eigen::LDLT ldlt(aMap); + if (ldlt.info() != Eigen::Success) + { + itkGenericExceptionMacro("itk::Math::InverseSymmetric failed; input is likely non-finite (NaN/Inf)."); + } + const auto dAbs = ldlt.vectorD().cwiseAbs().eval(); + const TReal dMax = dAbs.maxCoeff(); + if (dMax == TReal{ 0 } || dAbs.minCoeff() <= dMax * static_cast(n) * std::numeric_limits::epsilon()) + { + itkGenericExceptionMacro("itk::Math::InverseSymmetric failed; input matrix is singular."); + } + Eigen::Map(invData, n, n) = ldlt.solve(ColMajor::Identity(n, n)); +} +} // namespace detail + +/** \brief Solve the symmetric linear system \c A x = b via LDLT, backed by Eigen. + * + * Eigen-backed alternative to vnl_cholesky for a symmetric \a A (positive + * definite or indefinite). Uses Eigen's robust LDLT (Bunch-Kaufman symmetric + * pivoting), which handles near-singular and indefinite symmetric systems + * gracefully at Cholesky-class cost -- so a single call replaces the common + * cholesky-plus-SVD-fallback idiom. Only \a A's symmetry is assumed; the caller + * supplies a symmetric matrix. + * + * The interface uses ITK container types (no vnl or Eigen type appears in the + * signature): a runtime-sized \c itk::Array2D / \c itk::Array pair, or a + * fixed-size \c itk::Matrix / \c itk::Vector pair. + * + * A non-finite input throws an itk::ExceptionObject. + * + * \ingroup ITKCommon + */ +template +Array +SolveSymmetric(const Array2D & A, const Array & b) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n || b.size() != n) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric requires a non-empty square A and a matching b."); + } + Array x(n); + detail::SolveSymmetricLDLTEigen(A.data_block(), b.data_block(), n, x.data_block()); + return x; +} + +/** Fixed-size symmetric solve A x = b via LDLT. */ +template +Vector +SolveSymmetric(const Matrix & A, const Vector & b) +{ + Vector x; + detail::SolveSymmetricLDLTEigen(A.GetVnlMatrix().data_block(), b.GetDataPointer(), VDim, x.GetDataPointer()); + return x; +} + +// --- vnl convenience overloads --------------------------------------------- +// Provided so consumers that still hold vnl types (e.g. legacy filters) can call +// directly without wrapping. The ITK-typed overloads above are the forward- +// looking interface; the very-long-term goal is to eliminate vnl from the ITK +// API, so prefer the ITK-typed signatures in new code. + +/** Runtime-sized symmetric solve on vnl types (convenience). */ +template +vnl_vector +SolveSymmetric(const vnl_matrix & A, const vnl_vector & b) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n || b.size() != n) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric requires a non-empty square A and a matching b."); + } + vnl_vector x(n); + detail::SolveSymmetricLDLTEigen(A.data_block(), b.data_block(), n, x.data_block()); + return x; +} + +/** Fixed-size symmetric solve on vnl types (convenience). */ +template +vnl_vector_fixed +SolveSymmetric(const vnl_matrix_fixed & A, const vnl_vector_fixed & b) +{ + vnl_vector_fixed x; + detail::SolveSymmetricLDLTEigen(A.data_block(), b.data_block(), VDim, x.data_block()); + return x; +} + +// --- multi-RHS solve and symmetric inverse --------------------------------- + +/** Solve the symmetric system \c A X = B for a matrix right-hand side \a B via a + * single LDLT factorization (all columns solved together -- O(n^3), not one + * factorization per column). \a A must be symmetric. */ +template +Array2D +SolveSymmetric(const Array2D & A, const Array2D & B) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n || B.rows() != n) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric requires a non-empty square A and a matching B."); + } + Array2D X(n, B.cols()); + detail::SolveSymmetricMatrixLDLTEigen(A.data_block(), B.data_block(), n, B.cols(), X.data_block()); + return X; +} + +/** \brief Inverse of a symmetric matrix, backed by Eigen LDLT. + * + * Returns \c A^-1 for a symmetric \a A via a single LDLT factorization (solves + * A X = I). Prefer the solve overloads when only \c A^-1 b or \c A^-1 B is + * needed -- forming the full inverse is rarely necessary. A non-finite or + * singular input throws an itk::ExceptionObject. + * + * \ingroup ITKCommon + */ +template +Array2D +InverseSymmetric(const Array2D & A) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n) + { + itkGenericExceptionMacro("itk::Math::InverseSymmetric requires a non-empty square A."); + } + Array2D inv(n, n); + detail::InverseSymmetricLDLTEigen(A.data_block(), n, inv.data_block()); + return inv; +} + +/** Multi-RHS symmetric solve on vnl types (convenience). */ +template +vnl_matrix +SolveSymmetric(const vnl_matrix & A, const vnl_matrix & B) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n || B.rows() != n) + { + itkGenericExceptionMacro("itk::Math::SolveSymmetric requires a non-empty square A and a matching B."); + } + vnl_matrix X(n, B.cols()); + detail::SolveSymmetricMatrixLDLTEigen(A.data_block(), B.data_block(), n, B.cols(), X.data_block()); + return X; +} + +/** Symmetric inverse on vnl types (convenience). */ +template +vnl_matrix +InverseSymmetric(const vnl_matrix & A) +{ + const unsigned int n = A.rows(); + if (n == 0 || A.cols() != n) + { + itkGenericExceptionMacro("itk::Math::InverseSymmetric requires a non-empty square A."); + } + vnl_matrix inv(n, n); + detail::InverseSymmetricLDLTEigen(A.data_block(), n, inv.data_block()); + return inv; +} + +} // namespace Math +} // namespace itk + +#endif // itkMathLDLT_h diff --git a/Modules/Core/Common/test/CMakeLists.txt b/Modules/Core/Common/test/CMakeLists.txt index 53be9fce787..e13484091f1 100644 --- a/Modules/Core/Common/test/CMakeLists.txt +++ b/Modules/Core/Common/test/CMakeLists.txt @@ -1485,6 +1485,7 @@ set( itkMathRoundGTest.cxx itkMathVnlParityGTest.cxx itkMathSVDGTest.cxx + itkMathLDLTGTest.cxx itkVnlSVDEngineGTest.cxx itkMatrixExponentialGTest.cxx itkMatrixGTest.cxx diff --git a/Modules/Core/Common/test/itkMathLDLTGTest.cxx b/Modules/Core/Common/test/itkMathLDLTGTest.cxx new file mode 100644 index 00000000000..b54f18ce5eb --- /dev/null +++ b/Modules/Core/Common/test/itkMathLDLTGTest.cxx @@ -0,0 +1,242 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +// First include the header file to be tested: +#include "itkMathLDLT.h" + +// Exercise the deprecated VNL engines as cross-check oracles. +// They are unavailable under ITK_FUTURE_LEGACY_REMOVE. +#ifndef ITK_FUTURE_LEGACY_REMOVE +# define ITK_LEGACY_TEST +# include "vnl/algo/vnl_cholesky.h" +# include "vnl/algo/vnl_matrix_inverse.h" +#endif + +#include +#include + +namespace +{ +// Symmetric positive-definite matrix G^T G + shift*I, as an itk::Array2D. +template +itk::Array2D +MakeSPD(unsigned int n, T shift) +{ + vnl_matrix G(n + 2, n); + for (unsigned int i = 0; i < n + 2; ++i) + { + for (unsigned int j = 0; j < n; ++j) + { + G(i, j) = static_cast(std::sin(0.7 * i + 1.3 * j)); + } + } + vnl_matrix A = G.transpose() * G; + for (unsigned int i = 0; i < n; ++i) + { + A(i, i) += shift; + } + return itk::Array2D(A); +} +} // namespace + +// The capability macro must be defined by the header. +#ifndef ITK_MATH_HAS_SOLVE_SYMMETRIC +# error "itkMathLDLT.h must define ITK_MATH_HAS_SOLVE_SYMMETRIC" +#endif + +TEST(MathLDLT, DynamicMatchesVnlCholeskyOnSPD) +{ + for (unsigned int n : { 1u, 2u, 3u, 6u, 12u, 20u }) + { + const itk::Array2D A = MakeSPD(n, 0.5); + itk::Array b(n); + for (unsigned int i = 0; i < n; ++i) + { + b[i] = std::cos(0.3 * i + 1.0); + } + + const itk::Array x = itk::Math::SolveSymmetric(A, b); + ASSERT_EQ(x.size(), n); +#ifndef ITK_FUTURE_LEGACY_REMOVE + // Reference: itk::Array2D / itk::Array upcast to their vnl bases. + const vnl_vector ref = vnl_cholesky(A).solve(b); + for (unsigned int i = 0; i < n; ++i) + { + EXPECT_NEAR(x[i], ref[i], 1e-10) << "n=" << n << " i=" << i; + } +#endif + // Residual A x - b is near zero (vnl ops via the base classes). + const vnl_vector r = A * x - b; + EXPECT_LT(r.inf_norm(), 1e-9) << "n=" << n; + } +} + +TEST(MathLDLT, FixedMatchesDynamic) +{ + constexpr unsigned int VDim = 4; + const itk::Array2D Ad = MakeSPD(VDim, 0.3); + + itk::Matrix Af; + itk::Vector bf; + itk::Array bd(VDim); + for (unsigned int i = 0; i < VDim; ++i) + { + for (unsigned int j = 0; j < VDim; ++j) + { + Af(i, j) = Ad(i, j); + } + bf[i] = bd[i] = std::cos(0.3 * i + 1.0); + } + + const itk::Vector xf = itk::Math::SolveSymmetric(Af, bf); + const itk::Array xd = itk::Math::SolveSymmetric(Ad, bd); + for (unsigned int i = 0; i < VDim; ++i) + { + EXPECT_NEAR(xf[i], xd[i], 1e-12) << "i=" << i; + } +} + +// LDLT solves indefinite symmetric systems where a plain Cholesky (SPD-only) +// would fail -- the key robustness property motivating the JLF adoption. +TEST(MathLDLT, HandlesIndefiniteSymmetric) +{ + constexpr unsigned int n = 3; + itk::Matrix A; + A.Fill(0.0); + A(0, 0) = 2.0; + A(1, 1) = -1.0; // negative eigenvalue -> indefinite + A(2, 2) = 3.0; + A(0, 1) = A(1, 0) = 0.5; + itk::Vector b; + b[0] = 1.0; + b[1] = 2.0; + b[2] = -1.0; + + const itk::Vector x = itk::Math::SolveSymmetric(A, b); + const itk::Vector r = A * x - b; + EXPECT_LT(r.GetNorm(), 1e-10); +} + +// The vnl convenience overload (for consumers still holding vnl types, e.g. the +// ANTs joint-label-fusion MxBar) must agree with the ITK-typed path. +TEST(MathLDLT, VnlConvenienceMatchesItkTyped) +{ + constexpr unsigned int n = 8; + const itk::Array2D A = MakeSPD(n, 0.5); + itk::Array b(n); + for (unsigned int i = 0; i < n; ++i) + { + b[i] = std::cos(0.3 * i + 1.0); + } + + // vnl overload: A / b upcast to their vnl bases. + const vnl_matrix & Avnl = A; + const vnl_vector & bvnl = b; + const vnl_vector xvnl = itk::Math::SolveSymmetric(Avnl, bvnl); + const itk::Array xitk = itk::Math::SolveSymmetric(A, b); + + ASSERT_EQ(xvnl.size(), n); + for (unsigned int i = 0; i < n; ++i) + { + EXPECT_NEAR(xvnl[i], xitk[i], 1e-12) << "i=" << i; + } +} + +TEST(MathLDLT, InverseSymmetricMatchesVnl) +{ + for (unsigned int n : { 1u, 3u, 6u, 12u }) + { + const itk::Array2D A = MakeSPD(n, 0.5); + const itk::Array2D inv = itk::Math::InverseSymmetric(A); + ASSERT_EQ(inv.rows(), n); +#ifndef ITK_FUTURE_LEGACY_REMOVE + const vnl_matrix ref = vnl_matrix_inverse(A).inverse(); + for (unsigned int i = 0; i < n; ++i) + { + for (unsigned int j = 0; j < n; ++j) + { + EXPECT_NEAR(inv(i, j), ref(i, j), 1e-9) << "n=" << n << " (" << i << "," << j << ")"; + } + } +#endif + // A * inv == I + const vnl_matrix prod = A * inv; + for (unsigned int i = 0; i < n; ++i) + { + for (unsigned int j = 0; j < n; ++j) + { + EXPECT_NEAR(prod(i, j), (i == j ? 1.0 : 0.0), 1e-9); + } + } + } +} + +TEST(MathLDLT, InverseSymmetricThrowsOnSingular) +{ + // Rank-1 symmetric matrix [[1,2],[2,4]] is exactly singular. + itk::Array2D A(2, 2); + A(0, 0) = 1.0; + A(0, 1) = A(1, 0) = 2.0; + A(1, 1) = 4.0; + EXPECT_THROW(itk::Math::InverseSymmetric(A), itk::ExceptionObject); + + const vnl_matrix & Avnl = A; + EXPECT_THROW(itk::Math::InverseSymmetric(Avnl), itk::ExceptionObject); + + itk::Array2D zero(3, 3); + zero.fill(0.0); + EXPECT_THROW(itk::Math::InverseSymmetric(zero), itk::ExceptionObject); +} + +TEST(MathLDLT, MatrixRhsSolveMatchesColumnwise) +{ + constexpr unsigned int n = 6; + const itk::Array2D A = MakeSPD(n, 0.4); + itk::Array2D B(n, 3); + for (unsigned int i = 0; i < n; ++i) + { + for (unsigned int j = 0; j < 3; ++j) + { + B(i, j) = std::cos(0.3 * i + 0.7 * j); + } + } + + const itk::Array2D X = itk::Math::SolveSymmetric(A, B); + // A X == B + const vnl_matrix prod = A * X; + for (unsigned int i = 0; i < n; ++i) + { + for (unsigned int j = 0; j < 3; ++j) + { + EXPECT_NEAR(prod(i, j), B(i, j), 1e-9) << "(" << i << "," << j << ")"; + } + } +} + +TEST(MathLDLT, FloatPrecision) +{ + const itk::Array2D A = MakeSPD(5, 0.5f); + itk::Array b(5); + for (unsigned int i = 0; i < 5; ++i) + { + b[i] = static_cast(std::cos(0.3 * i)); + } + const itk::Array x = itk::Math::SolveSymmetric(A, b); + const vnl_vector r = A * x - b; + EXPECT_LT(r.inf_norm(), 1e-4f); +} From 4810e1c855795898ead602ab70618417e4a4abf8 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Wed, 8 Jul 2026 20:15:12 -0500 Subject: [PATCH 2/3] ENH: Add vnl_cholesky engine coverage as itkVnlCholeskyEngineGTest Adapt vxl core/vnl/algo/tests/test_cholesky.cxx into an ITK GoogleTest so the coverage runs in ITK CI and guards the native Cholesky engine that replaced the netlib LINPACK path (98f7498a3a9). Beyond the scavenged cases (inverse vs direct inverse, two-sided identity in both modes, known-solution solve), adds determinant, triangular-factor, rank-deficiency/rcond, and an equivalence check against the Eigen-backed itk::Math::SolveSymmetricPositiveDefinite. The suite defines ITK_LEGACY_TEST and compiles out under ITK_FUTURE_LEGACY_REMOVE, matching the vnl_cholesky deprecation guard. Follows the itkVnlSVDEngineGTest precedent (PR #6546). Tracking: #6403. --- Modules/Core/Common/test/CMakeLists.txt | 1 + .../Common/test/itkVnlCholeskyEngineGTest.cxx | 189 ++++++++++++++++++ 2 files changed, 190 insertions(+) create mode 100644 Modules/Core/Common/test/itkVnlCholeskyEngineGTest.cxx diff --git a/Modules/Core/Common/test/CMakeLists.txt b/Modules/Core/Common/test/CMakeLists.txt index e13484091f1..1593d8456f6 100644 --- a/Modules/Core/Common/test/CMakeLists.txt +++ b/Modules/Core/Common/test/CMakeLists.txt @@ -1486,6 +1486,7 @@ set( itkMathVnlParityGTest.cxx itkMathSVDGTest.cxx itkMathLDLTGTest.cxx + itkVnlCholeskyEngineGTest.cxx itkVnlSVDEngineGTest.cxx itkMatrixExponentialGTest.cxx itkMatrixGTest.cxx diff --git a/Modules/Core/Common/test/itkVnlCholeskyEngineGTest.cxx b/Modules/Core/Common/test/itkVnlCholeskyEngineGTest.cxx new file mode 100644 index 00000000000..908323d9a43 --- /dev/null +++ b/Modules/Core/Common/test/itkVnlCholeskyEngineGTest.cxx @@ -0,0 +1,189 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +// Exercises the vnl_cholesky engine end-to-end (inverse, solve, determinant, +// triangular factors, condition estimate), adapted from vxl +// core/vnl/algo/tests/test_cholesky.cxx so the coverage runs in ITK CI and +// guards any future change of the underlying Cholesky engine. + +#include "itkCholeskySolve.h" + +// The deprecated VNL engine under test is unavailable under ITK_FUTURE_LEGACY_REMOVE. +#ifndef ITK_FUTURE_LEGACY_REMOVE +# define ITK_LEGACY_TEST +# include "vnl/algo/vnl_cholesky.h" +# include "vnl/vnl_inverse.h" +# include "vnl/vnl_random.h" + +# include +# include +# include + +namespace +{ +// Deterministic symmetric positive-definite matrix A = R R^T. +vnl_matrix +MakeSPD(unsigned int n, unsigned int seed) +{ + vnl_random rng(seed); + vnl_matrix R(n, n); + for (unsigned int i = 0; i < n; ++i) + { + for (unsigned int j = 0; j < n; ++j) + { + R(i, j) = rng.drand32(-1.0, 1.0); + } + } + return R * R.transpose(); +} +} // namespace + + +// Scavenged from vxl test_cholesky: cholesky inverse matches the direct inverse. +TEST(VnlCholeskyEngine, InverseMatchesDirectInverse) +{ + const vnl_matrix A = MakeSPD(3, 1000); + const vnl_cholesky chol(A); + EXPECT_NEAR((chol.inverse() - vnl_inverse(A)).fro_norm(), 0.0, 1e-10); +} + + +// Scavenged from vxl test_cholesky: inverse is a two-sided identity, in both +// the default and estimate_condition modes. +TEST(VnlCholeskyEngine, InverseIsTwoSidedIdentity) +{ + const vnl_matrix A = MakeSPD(3, 1000); + vnl_matrix identity(3, 3); + identity.set_identity(); + + { + const vnl_cholesky chol(A); + EXPECT_NEAR((chol.inverse() * A - identity).fro_norm(), 0.0, 1e-10); + EXPECT_NEAR((A * chol.inverse() - identity).fro_norm(), 0.0, 1e-10); + } + { + const vnl_cholesky chol(A, vnl_cholesky::estimate_condition); + EXPECT_NEAR((chol.inverse() * A - identity).fro_norm(), 0.0, 1e-10); + EXPECT_NEAR((A * chol.inverse() - identity).fro_norm(), 0.0, 1e-10); + } +} + + +// Scavenged from vxl test_cholesky: solve recovers a known solution; the +// out-parameter overload agrees with the returning overload. +TEST(VnlCholeskyEngine, SolveRecoversKnownSolution) +{ + const vnl_matrix A = MakeSPD(3, 1000); + vnl_random rng(2000); + vnl_vector x0(3); + for (unsigned int i = 0; i < 3; ++i) + { + x0[i] = rng.drand32(-1.0, 1.0); + } + const vnl_vector b = A * x0; + + const vnl_cholesky chol(A); + const vnl_vector x = chol.solve(b); + EXPECT_NEAR((x - x0).one_norm(), 0.0, 1e-6); + + vnl_vector xOut(3); + chol.solve(b, &xOut); + EXPECT_NEAR((xOut - x).one_norm(), 0.0, 1e-14); +} + + +// determinant() equals the analytic determinant of L0 L0^T for a known lower +// factor L0: det = (prod diag(L0))^2. +TEST(VnlCholeskyEngine, DeterminantMatchesAnalytic) +{ + constexpr unsigned int n = 4; + vnl_matrix L0(n, n, 0.0); + double expected = 1.0; + for (unsigned int i = 0; i < n; ++i) + { + L0(i, i) = static_cast(i + 1); + expected *= L0(i, i) * L0(i, i); + for (unsigned int j = 0; j < i; ++j) + { + L0(i, j) = 0.25 * static_cast(i + j + 1); + } + } + const vnl_matrix A = L0 * L0.transpose(); + + const vnl_cholesky chol(A); + EXPECT_NEAR(chol.determinant(), expected, 1e-9 * expected); +} + + +// lower_triangle()/upper_triangle() reconstruct A and are transposes. +TEST(VnlCholeskyEngine, TriangularFactorsReconstruct) +{ + const vnl_matrix A = MakeSPD(5, 3000); + const vnl_cholesky chol(A); + const vnl_matrix L = chol.lower_triangle(); + const vnl_matrix U = chol.upper_triangle(); + + EXPECT_NEAR((L * L.transpose() - A).fro_norm() / A.fro_norm(), 0.0, 1e-12); + EXPECT_NEAR((U - L.transpose()).fro_norm(), 0.0, 1e-14); + for (unsigned int i = 0; i < 5; ++i) + { + for (unsigned int j = i + 1; j < 5; ++j) + { + EXPECT_EQ(L(i, j), 0.0); + } + } +} + + +// SPD input factors fully (rank_deficiency 0) with a usable condition estimate; +// an indefinite input reports non-positive-definiteness. +TEST(VnlCholeskyEngine, RankDeficiencyAndConditionEstimate) +{ + const vnl_matrix A = MakeSPD(4, 4000); + const vnl_cholesky chol(A, vnl_cholesky::estimate_condition); + EXPECT_EQ(chol.rank_deficiency(), 0); + EXPECT_GT(chol.rcond(), std::sqrt(std::numeric_limits::epsilon())); + EXPECT_LE(chol.rcond(), 1.0); + + vnl_matrix indefinite(2, 2); + indefinite(0, 0) = 1.0; + indefinite(0, 1) = 2.0; + indefinite(1, 0) = 2.0; + indefinite(1, 1) = 1.0; // eigenvalues 3, -1 + const vnl_cholesky badChol(indefinite, vnl_cholesky::quiet); + EXPECT_GT(badChol.rank_deficiency(), 0); +} + + +// The engine agrees with the Eigen-backed replacement API it is deprecated for. +TEST(VnlCholeskyEngine, EquivalentToItkCholeskySolve) +{ + const vnl_matrix A = MakeSPD(6, 5000); + vnl_vector b(6); + for (unsigned int i = 0; i < 6; ++i) + { + b[i] = std::cos(0.5 * (i + 1)); + } + + const vnl_cholesky chol(A, vnl_cholesky::quiet); + const vnl_vector xVnl = chol.solve(b); + const vnl_vector xItk = itk::Math::SolveSymmetricPositiveDefinite(A, b); + EXPECT_LT((xItk - xVnl).two_norm() / xVnl.two_norm(), 1e-10); +} + +#endif // ITK_FUTURE_LEGACY_REMOVE From 210fb5c1b9b79c1567dfeb4f52d1a1ae0cc5ddee Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Wed, 8 Jul 2026 20:15:14 -0500 Subject: [PATCH 3/3] STYLE: Brace single-line loop bodies in itkCholeskySolveGTest --- .../Core/Common/test/itkCholeskySolveGTest.cxx | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/Modules/Core/Common/test/itkCholeskySolveGTest.cxx b/Modules/Core/Common/test/itkCholeskySolveGTest.cxx index f91267a5406..d9ea71892ec 100644 --- a/Modules/Core/Common/test/itkCholeskySolveGTest.cxx +++ b/Modules/Core/Common/test/itkCholeskySolveGTest.cxx @@ -38,11 +38,17 @@ MakeSPD(unsigned int n) { vnl_matrix R(n, n); for (unsigned int i = 0; i < n; ++i) + { for (unsigned int j = 0; j < n; ++j) + { R(i, j) = static_cast(std::sin(0.7 * (i + 1) * (j + 2)) + 0.3 * (i + 1)); + } + } vnl_matrix A = R * R.transpose(); for (unsigned int i = 0; i < n; ++i) + { A(i, i) += static_cast(n); + } return A; } } // namespace @@ -55,7 +61,9 @@ TEST(CholeskySolve, SolveResidual) const vnl_matrix A = MakeSPD(n); vnl_vector b(n); for (unsigned int i = 0; i < n; ++i) + { b[i] = static_cast(i) - 1.5; + } const vnl_vector x = itk::Math::SolveSymmetricPositiveDefinite(A, b); const vnl_vector residual = A * x - b; @@ -72,8 +80,12 @@ TEST(CholeskySolve, LowerTriangleReconstructsMatrix) // L is lower triangular. for (unsigned int i = 0; i < n; ++i) + { for (unsigned int j = i + 1; j < n; ++j) + { EXPECT_NEAR(L(i, j), 0.0, 1e-12); + } + } const vnl_matrix reconstructed = L * L.transpose(); EXPECT_LT((reconstructed - A).fro_norm() / A.fro_norm(), 1e-12); @@ -88,7 +100,9 @@ TEST(CholeskySolve, EquivalentToVnlCholesky) const vnl_matrix A = MakeSPD(n); vnl_vector b(n); for (unsigned int i = 0; i < n; ++i) + { b[i] = std::cos(0.5 * (i + 1)); + } const vnl_vector xItk = itk::Math::SolveSymmetricPositiveDefinite(A, b); @@ -107,7 +121,9 @@ TEST(CholeskySolve, FloatResidual) const vnl_matrix A = MakeSPD(n); vnl_vector b(n); for (unsigned int i = 0; i < n; ++i) + { b[i] = static_cast(i) + 0.25f; + } const vnl_vector x = itk::Math::SolveSymmetricPositiveDefinite(A, b); const vnl_vector residual = A * x - b;