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 <limits>
+
+/** Capability macro. Downstream code selects the Eigen-backed symmetric solve
+ * with, e.g.:
+ * \code
+ *   #if __has_include(<itkMathLDLT.h>)
+ *   #  include <itkMathLDLT.h>
+ *   #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 <typename TReal>
+void
+SolveSymmetricLDLTEigen(const TReal * aData, const TReal * bData, unsigned int n, TReal * xData)
+{
+  using ColMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic>;
+  using Vector = Eigen::Matrix<TReal, Eigen::Dynamic, 1>;
+
+  const Eigen::Map<const ColMajor> aMap(aData, n, n);
+  const Eigen::Map<const Vector>   bMap(bData, n);
+
+  const Eigen::LDLT<ColMajor> ldlt(aMap);
+  if (ldlt.info() != Eigen::Success)
+  {
+    itkGenericExceptionMacro("itk::Math::SolveSymmetric failed; input is likely non-finite (NaN/Inf).");
+  }
+  Eigen::Map<Vector>(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 <typename TReal>
+void
+SolveSymmetricMatrixLDLTEigen(const TReal * aData, const TReal * bData, unsigned int n, unsigned int m, TReal * xData)
+{
+  using ColMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic>;
+  using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+  const Eigen::Map<const ColMajor> aMap(aData, n, n);
+  const Eigen::Map<const RowMajor> bMap(bData, n, m);
+
+  const Eigen::LDLT<ColMajor> ldlt(aMap);
+  if (ldlt.info() != Eigen::Success)
+  {
+    itkGenericExceptionMacro("itk::Math::SolveSymmetric failed; input is likely non-finite (NaN/Inf).");
+  }
+  Eigen::Map<RowMajor>(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 <typename TReal>
+void
+InverseSymmetricLDLTEigen(const TReal * aData, unsigned int n, TReal * invData)
+{
+  using ColMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic>;
+  using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+  const Eigen::Map<const ColMajor> aMap(aData, n, n);
+
+  const Eigen::LDLT<ColMajor> 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<TReal>(n) * std::numeric_limits<TReal>::epsilon())
+  {
+    itkGenericExceptionMacro("itk::Math::InverseSymmetric failed; input matrix is singular.");
+  }
+  Eigen::Map<RowMajor>(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 <typename TReal>
+Array<TReal>
+SolveSymmetric(const Array2D<TReal> & A, const Array<TReal> & 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<TReal> x(n);
+  detail::SolveSymmetricLDLTEigen<TReal>(A.data_block(), b.data_block(), n, x.data_block());
+  return x;
+}
+
+/** Fixed-size symmetric solve A x = b via LDLT. */
+template <typename TReal, unsigned int VDim>
+Vector<TReal, VDim>
+SolveSymmetric(const Matrix<TReal, VDim, VDim> & A, const Vector<TReal, VDim> & b)
+{
+  Vector<TReal, VDim> x;
+  detail::SolveSymmetricLDLTEigen<TReal>(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 <typename TReal>
+vnl_vector<TReal>
+SolveSymmetric(const vnl_matrix<TReal> & A, const vnl_vector<TReal> & 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<TReal> x(n);
+  detail::SolveSymmetricLDLTEigen<TReal>(A.data_block(), b.data_block(), n, x.data_block());
+  return x;
+}
+
+/** Fixed-size symmetric solve on vnl types (convenience). */
+template <typename TReal, unsigned int VDim>
+vnl_vector_fixed<TReal, VDim>
+SolveSymmetric(const vnl_matrix_fixed<TReal, VDim, VDim> & A, const vnl_vector_fixed<TReal, VDim> & b)
+{
+  vnl_vector_fixed<TReal, VDim> x;
+  detail::SolveSymmetricLDLTEigen<TReal>(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 <typename TReal>
+Array2D<TReal>
+SolveSymmetric(const Array2D<TReal> & A, const Array2D<TReal> & 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<TReal> X(n, B.cols());
+  detail::SolveSymmetricMatrixLDLTEigen<TReal>(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 <typename TReal>
+Array2D<TReal>
+InverseSymmetric(const Array2D<TReal> & A)
+{
+  const unsigned int n = A.rows();
+  if (n == 0 || A.cols() != n)
+  {
+    itkGenericExceptionMacro("itk::Math::InverseSymmetric requires a non-empty square A.");
+  }
+  Array2D<TReal> inv(n, n);
+  detail::InverseSymmetricLDLTEigen<TReal>(A.data_block(), n, inv.data_block());
+  return inv;
+}
+
+/** Multi-RHS symmetric solve on vnl types (convenience). */
+template <typename TReal>
+vnl_matrix<TReal>
+SolveSymmetric(const vnl_matrix<TReal> & A, const vnl_matrix<TReal> & 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<TReal> X(n, B.cols());
+  detail::SolveSymmetricMatrixLDLTEigen<TReal>(A.data_block(), B.data_block(), n, B.cols(), X.data_block());
+  return X;
+}
+
+/** Symmetric inverse on vnl types (convenience). */
+template <typename TReal>
+vnl_matrix<TReal>
+InverseSymmetric(const vnl_matrix<TReal> & 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<TReal> inv(n, n);
+  detail::InverseSymmetricLDLTEigen<TReal>(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..1593d8456f6 100644
--- a/Modules/Core/Common/test/CMakeLists.txt
+++ b/Modules/Core/Common/test/CMakeLists.txt
@@ -1485,6 +1485,8 @@ set(
   itkMathRoundGTest.cxx
   itkMathVnlParityGTest.cxx
   itkMathSVDGTest.cxx
+  itkMathLDLTGTest.cxx
+  itkVnlCholeskyEngineGTest.cxx
   itkVnlSVDEngineGTest.cxx
   itkMatrixExponentialGTest.cxx
   itkMatrixGTest.cxx
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<T> R(n, n);
   for (unsigned int i = 0; i < n; ++i)
+  {
     for (unsigned int j = 0; j < n; ++j)
+    {
       R(i, j) = static_cast<T>(std::sin(0.7 * (i + 1) * (j + 2)) + 0.3 * (i + 1));
+    }
+  }
   vnl_matrix<T> A = R * R.transpose();
   for (unsigned int i = 0; i < n; ++i)
+  {
     A(i, i) += static_cast<T>(n);
+  }
   return A;
 }
 } // namespace
@@ -55,7 +61,9 @@ TEST(CholeskySolve, SolveResidual)
   const vnl_matrix<double> A = MakeSPD<double>(n);
   vnl_vector<double>       b(n);
   for (unsigned int i = 0; i < n; ++i)
+  {
     b[i] = static_cast<double>(i) - 1.5;
+  }
 
   const vnl_vector<double> x = itk::Math::SolveSymmetricPositiveDefinite(A, b);
   const vnl_vector<double> 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<double> 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<double> A = MakeSPD<double>(n);
   vnl_vector<double>       b(n);
   for (unsigned int i = 0; i < n; ++i)
+  {
     b[i] = std::cos(0.5 * (i + 1));
+  }
 
   const vnl_vector<double> xItk = itk::Math::SolveSymmetricPositiveDefinite(A, b);
 
@@ -107,7 +121,9 @@ TEST(CholeskySolve, FloatResidual)
   const vnl_matrix<float> A = MakeSPD<float>(n);
   vnl_vector<float>       b(n);
   for (unsigned int i = 0; i < n; ++i)
+  {
     b[i] = static_cast<float>(i) + 0.25f;
+  }
 
   const vnl_vector<float> x = itk::Math::SolveSymmetricPositiveDefinite(A, b);
   const vnl_vector<float> residual = A * x - b;
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 <gtest/gtest.h>
+#include <cmath>
+
+namespace
+{
+// Symmetric positive-definite matrix G^T G + shift*I, as an itk::Array2D.
+template <typename T>
+itk::Array2D<T>
+MakeSPD(unsigned int n, T shift)
+{
+  vnl_matrix<T> 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<T>(std::sin(0.7 * i + 1.3 * j));
+    }
+  }
+  vnl_matrix<T> A = G.transpose() * G;
+  for (unsigned int i = 0; i < n; ++i)
+  {
+    A(i, i) += shift;
+  }
+  return itk::Array2D<T>(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<double> A = MakeSPD<double>(n, 0.5);
+    itk::Array<double>         b(n);
+    for (unsigned int i = 0; i < n; ++i)
+    {
+      b[i] = std::cos(0.3 * i + 1.0);
+    }
+
+    const itk::Array<double> 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<double> 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<double> r = A * x - b;
+    EXPECT_LT(r.inf_norm(), 1e-9) << "n=" << n;
+  }
+}
+
+TEST(MathLDLT, FixedMatchesDynamic)
+{
+  constexpr unsigned int     VDim = 4;
+  const itk::Array2D<double> Ad = MakeSPD<double>(VDim, 0.3);
+
+  itk::Matrix<double, VDim, VDim> Af;
+  itk::Vector<double, VDim>       bf;
+  itk::Array<double>              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<double, VDim> xf = itk::Math::SolveSymmetric(Af, bf);
+  const itk::Array<double>        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<double, n, n> 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<double, n> b;
+  b[0] = 1.0;
+  b[1] = 2.0;
+  b[2] = -1.0;
+
+  const itk::Vector<double, n> x = itk::Math::SolveSymmetric(A, b);
+  const itk::Vector<double, n> 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<double> A = MakeSPD<double>(n, 0.5);
+  itk::Array<double>         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<double> & Avnl = A;
+  const vnl_vector<double> & bvnl = b;
+  const vnl_vector<double>   xvnl = itk::Math::SolveSymmetric(Avnl, bvnl);
+  const itk::Array<double>   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<double> A = MakeSPD<double>(n, 0.5);
+    const itk::Array2D<double> inv = itk::Math::InverseSymmetric(A);
+    ASSERT_EQ(inv.rows(), n);
+#ifndef ITK_FUTURE_LEGACY_REMOVE
+    const vnl_matrix<double> ref = vnl_matrix_inverse<double>(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<double> 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<double> 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<double> & Avnl = A;
+  EXPECT_THROW(itk::Math::InverseSymmetric(Avnl), itk::ExceptionObject);
+
+  itk::Array2D<double> 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<double> A = MakeSPD<double>(n, 0.4);
+  itk::Array2D<double>       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<double> X = itk::Math::SolveSymmetric(A, B);
+  // A X == B
+  const vnl_matrix<double> 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<float> A = MakeSPD<float>(5, 0.5f);
+  itk::Array<float>         b(5);
+  for (unsigned int i = 0; i < 5; ++i)
+  {
+    b[i] = static_cast<float>(std::cos(0.3 * i));
+  }
+  const itk::Array<float> x = itk::Math::SolveSymmetric(A, b);
+  const vnl_vector<float> r = A * x - b;
+  EXPECT_LT(r.inf_norm(), 1e-4f);
+}
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 <gtest/gtest.h>
+#  include <cmath>
+#  include <limits>
+
+namespace
+{
+// Deterministic symmetric positive-definite matrix A = R R^T.
+vnl_matrix<double>
+MakeSPD(unsigned int n, unsigned int seed)
+{
+  vnl_random         rng(seed);
+  vnl_matrix<double> 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<double> 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<double> A = MakeSPD(3, 1000);
+  vnl_matrix<double>       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<double> A = MakeSPD(3, 1000);
+  vnl_random               rng(2000);
+  vnl_vector<double>       x0(3);
+  for (unsigned int i = 0; i < 3; ++i)
+  {
+    x0[i] = rng.drand32(-1.0, 1.0);
+  }
+  const vnl_vector<double> b = A * x0;
+
+  const vnl_cholesky       chol(A);
+  const vnl_vector<double> x = chol.solve(b);
+  EXPECT_NEAR((x - x0).one_norm(), 0.0, 1e-6);
+
+  vnl_vector<double> 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<double>     L0(n, n, 0.0);
+  double                 expected = 1.0;
+  for (unsigned int i = 0; i < n; ++i)
+  {
+    L0(i, i) = static_cast<double>(i + 1);
+    expected *= L0(i, i) * L0(i, i);
+    for (unsigned int j = 0; j < i; ++j)
+    {
+      L0(i, j) = 0.25 * static_cast<double>(i + j + 1);
+    }
+  }
+  const vnl_matrix<double> 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<double> A = MakeSPD(5, 3000);
+  const vnl_cholesky       chol(A);
+  const vnl_matrix<double> L = chol.lower_triangle();
+  const vnl_matrix<double> 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<double> 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<double>::epsilon()));
+  EXPECT_LE(chol.rcond(), 1.0);
+
+  vnl_matrix<double> 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<double> A = MakeSPD(6, 5000);
+  vnl_vector<double>       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<double> xVnl = chol.solve(b);
+  const vnl_vector<double> xItk = itk::Math::SolveSymmetricPositiveDefinite(A, b);
+  EXPECT_LT((xItk - xVnl).two_norm() / xVnl.two_norm(), 1e-10);
+}
+
+#endif // ITK_FUTURE_LEGACY_REMOVE