Skip to content

Commit b657b05

Browse files
ooplesclaude
ooples
and
claude
authored
test(decomposition): add integration tests for 7 matrix decompositions (#674)
* test(decomposition): add integration tests for 7 matrix decompositions (#620) Add comprehensive integration tests for previously untested decompositions: - GramSchmidtDecomposition (Classical and Modified variants) - NmfDecomposition (Non-negative Matrix Factorization) - IcaDecomposition (Independent Component Analysis) - CramerDecomposition (Cramer's rule for linear systems) - NormalDecomposition (Normal equation method) - ComplexMatrixDecomposition (Complex number wrapper) - TakagiDecomposition (Takagi factorization with multiple algorithms) Tests verify: decomposition correctness, Solve/Invert functionality, algorithm variants, error handling, and numerical properties. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> * fix(decomposition): fix bugs found during integration testing - CholeskyDecomposition: use tolerance-based symmetry validation instead of exact floating-point comparison to handle A^T*A numerical errors - CramerDecomposition: add base case for 0×0 matrix determinant (returns 1 by empty product convention) to fix 1×1 matrix inversion - TakagiDecomposition: fix singular value calculation to use |eigenvalue| instead of sqrt(|eigenvalue|), and fix Solve method to use correct matrix operations x = U * S^(-1) * U^H * b 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> * fix: correct takagi decomposition jacobi algorithm - Fix tau formula: use diff = (app - aqq) instead of (aqq - app) The correct formula for zeroing A'[p,q] is tan(2θ) = 2*apq / (app - aqq) - Fix t formula: use t = tau / (sqrt(1 + tau²) + 1) [stable form] Previous formula t = 1 / (|tau| + sqrt(1 + tau²)) was incorrect - Add proper handling for equal diagonal elements (45-degree rotation) - Add detailed comments explaining the Jacobi rotation formulas - Update V' eigenvector accumulation with correct formulas - Build U matrix with proper phase handling for negative eigenvalues - Update test to expect correct singular values (eigenvalue magnitudes, not square roots) The Jacobi algorithm now correctly preserves eigenvalues and produces accurate reconstruction (error ~1e-14 machine precision). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> * docs: fix comment in cholesky validation Correct comment from "relative tolerance" to "absolute tolerance" to accurately describe the 1e-10 tolerance used for symmetry checking. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> * test: improve decomposition integration tests - Add 1x1 matrix inversion test for Cramer (enabled by 0x0 determinant fix) - Fix unused parameters in Theory test signatures (NMF, Normal) - Improve test names for GramSchmidt and ICA to better describe behavior - Use relative tolerance assertions in NMF convergence tests 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> * test: improve takagi solve and invert tests to verify correctness - Solve test now verifies A * x ≈ b, not just finite values - Invert test now verifies A * A^-1 ≈ I, not just finite values Addresses PR review comments about proper solution verification. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <[email protected]> --------- Co-authored-by: Claude Opus 4.5 <[email protected]>
1 parent 59d3be8 commit b657b05

10 files changed

+3205
-32
lines changed

‎src/DecompositionMethods/MatrixDecomposition/CholeskyDecomposition.cs‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -127,11 +127,16 @@ private Matrix<T> ComputeDecomposition(Matrix<T> matrix, CholeskyAlgorithmType a
127127
private void ValidateSymmetric(Matrix<T> matrix)
128128
{
129129
int n = matrix.Rows;
130+
// Use an absolute tolerance for floating-point comparison
131+
// This is necessary because computations like A^T*A may have tiny numerical errors
132+
T tolerance = NumOps.FromDouble(1e-10);
133+
130134
for (int i = 0; i < n; i++)
131135
{
132136
for (int j = i + 1; j < n; j++)
133137
{
134-
if (!NumOps.Equals(matrix[i, j], matrix[j, i]))
138+
T diff = NumOps.Abs(NumOps.Subtract(matrix[i, j], matrix[j, i]));
139+
if (NumOps.GreaterThan(diff, tolerance))
135140
{
136141
throw new ArgumentException("Matrix must be symmetric for Cholesky decomposition.");
137142
}
@@ -168,9 +173,14 @@ private Matrix<T> ComputeCholeskyDefault(Matrix<T> matrix)
168173
{
169174
for (int j = 0; j <= i; j++)
170175
{
171-
if (i != j && !NumOps.Equals(matrix[i, j], matrix[j, i]))
176+
if (i != j)
172177
{
173-
throw new ArgumentException("Matrix must be symmetric for Cholesky decomposition.");
178+
T diff = NumOps.Abs(NumOps.Subtract(matrix[i, j], matrix[j, i]));
179+
T tolerance = NumOps.FromDouble(1e-10);
180+
if (NumOps.GreaterThan(diff, tolerance))
181+
{
182+
throw new ArgumentException("Matrix must be symmetric for Cholesky decomposition.");
183+
}
174184
}
175185

176186
if (j == i) // Diagonal elements

‎src/DecompositionMethods/MatrixDecomposition/CramerDecomposition.cs‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,13 @@ public override Matrix<T> Invert()
137137
/// <returns>The determinant value.</returns>
138138
private T Determinant(Matrix<T> matrix)
139139
{
140+
// Base case: 0x0 matrix (empty matrix) has determinant 1 by convention
141+
// This is needed for cofactor expansion of 1x1 matrices to work correctly
142+
if (matrix.Rows == 0)
143+
{
144+
return NumOps.One;
145+
}
146+
140147
if (matrix.Rows == 1)
141148
{
142149
return matrix[0, 0];

‎src/DecompositionMethods/MatrixDecomposition/TakagiDecomposition.cs‎

Lines changed: 167 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -121,7 +121,8 @@ protected override void Decompose()
121121
{
122122
int n = matrix.Rows;
123123
var S = new Matrix<T>(n, n);
124-
var U = Matrix<Complex<T>>.CreateIdentity(n);
124+
// Start with real identity for eigenvector accumulation
125+
var V = Matrix<T>.CreateIdentityMatrix(n);
125126
var A = matrix.Clone();
126127

127128
const int maxIterations = 100;
@@ -152,20 +153,53 @@ protected override void Decompose()
152153
break;
153154
}
154155

155-
// Compute the Jacobi rotation
156+
// Compute the Jacobi rotation angle
156157
T app = A[p, p];
157158
T aqq = A[q, q];
158159
T apq = A[p, q];
159-
T theta = NumOps.Divide(NumOps.Subtract(app, aqq), NumOps.Multiply(NumOps.FromDouble(2), apq));
160-
T t = NumOps.Divide(NumOps.FromDouble(1), NumOps.Add(NumOps.Abs(theta), NumOps.Sqrt(NumOps.Add(NumOps.Square(theta), NumOps.One))));
161-
if (NumOps.LessThan(theta, NumOps.Zero))
160+
161+
// Handle case when app == aqq
162+
T diff = NumOps.Subtract(app, aqq); // Note: (app - aqq), not (aqq - app)
163+
T c, s;
164+
if (NumOps.LessThan(NumOps.Abs(diff), tolerance))
162165
{
163-
t = NumOps.Negate(t);
166+
// When diagonal elements are equal, use 45-degree rotation
167+
c = NumOps.FromDouble(Math.Sqrt(0.5));
168+
s = NumOps.FromDouble(Math.Sqrt(0.5));
169+
if (NumOps.LessThan(apq, NumOps.Zero))
170+
{
171+
s = NumOps.Negate(s);
172+
}
173+
}
174+
else
175+
{
176+
// For A' = G^T * A * G to zero out A'[p,q], we need:
177+
// A'[p,q] = (c² - s²)*apq + cs*(aqq - app) = 0
178+
// This gives tan(2θ) = 2*apq / (app - aqq)
179+
// For Jacobi eigenvalue algorithm, t = tan(θ) where:
180+
// t = tau / (sqrt(1 + tau²) + 1) [numerically stable form]
181+
T tau = NumOps.Divide(NumOps.Multiply(NumOps.FromDouble(2), apq), diff);
182+
T sqrtTerm = NumOps.Sqrt(NumOps.Add(NumOps.One, NumOps.Square(tau)));
183+
T t = NumOps.Divide(tau, NumOps.Add(sqrtTerm, NumOps.One));
184+
c = NumOps.Divide(NumOps.One, NumOps.Sqrt(NumOps.Add(NumOps.One, NumOps.Square(t))));
185+
s = NumOps.Multiply(t, c);
164186
}
165-
T c = NumOps.Divide(NumOps.FromDouble(1), NumOps.Sqrt(NumOps.Add(NumOps.Square(t), NumOps.One)));
166-
T s = NumOps.Multiply(t, c);
167187

168-
// Update A
188+
// Update A with correct Jacobi rotation formulas
189+
// For G^T * A * G where G has G[p,q] = -s, G[q,p] = s:
190+
// A'[p,p] = c²*app + 2*c*s*apq + s²*aqq
191+
// A'[q,q] = s²*app - 2*c*s*apq + c²*aqq
192+
T c2 = NumOps.Square(c);
193+
T s2 = NumOps.Square(s);
194+
T cs2 = NumOps.Multiply(NumOps.Multiply(NumOps.FromDouble(2), c), s);
195+
T csapq = NumOps.Multiply(cs2, apq);
196+
197+
T newApp = NumOps.Add(NumOps.Add(NumOps.Multiply(c2, app), NumOps.Multiply(s2, aqq)), csapq);
198+
T newAqq = NumOps.Subtract(NumOps.Add(NumOps.Multiply(s2, app), NumOps.Multiply(c2, aqq)), csapq);
199+
200+
// Update off-diagonal elements: A' = G^T * A * G
201+
// A'[i,p] = c*A[i,p] + s*A[i,q] (for i != p,q)
202+
// A'[i,q] = -s*A[i,p] + c*A[i,q]
169203
for (int i = 0; i < n; i++)
170204
{
171205
if (i != p && i != q)
@@ -178,27 +212,49 @@ protected override void Decompose()
178212
A[i, q] = A[q, i];
179213
}
180214
}
181-
A[p, p] = NumOps.Add(NumOps.Multiply(NumOps.Square(c), app), NumOps.Multiply(NumOps.Square(s), aqq));
182-
A[q, q] = NumOps.Add(NumOps.Multiply(NumOps.Square(s), app), NumOps.Multiply(NumOps.Square(c), aqq));
215+
A[p, p] = newApp;
216+
A[q, q] = newAqq;
183217
A[p, q] = NumOps.Zero;
184218
A[q, p] = NumOps.Zero;
185219

186-
// Update U
220+
// Update eigenvectors V: V' = V * G
221+
// V'[:,p] = c*V[:,p] + s*V[:,q]
222+
// V'[:,q] = -s*V[:,p] + c*V[:,q]
187223
for (int i = 0; i < n; i++)
188224
{
189-
Complex<T> uip = U[i, p];
190-
Complex<T> uiq = U[i, q];
191-
U[i, p] = new Complex<T>(NumOps.Add(NumOps.Multiply(c, uip.Real), NumOps.Multiply(s, uiq.Real)),
192-
NumOps.Add(NumOps.Multiply(c, uip.Imaginary), NumOps.Multiply(s, uiq.Imaginary)));
193-
U[i, q] = new Complex<T>(NumOps.Subtract(NumOps.Multiply(c, uiq.Real), NumOps.Multiply(s, uip.Real)),
194-
NumOps.Subtract(NumOps.Multiply(c, uiq.Imaginary), NumOps.Multiply(s, uip.Imaginary)));
225+
T vip = V[i, p];
226+
T viq = V[i, q];
227+
V[i, p] = NumOps.Add(NumOps.Multiply(c, vip), NumOps.Multiply(s, viq));
228+
V[i, q] = NumOps.Subtract(NumOps.Multiply(c, viq), NumOps.Multiply(s, vip));
195229
}
196230
}
197231

198-
// Extract singular values
199-
for (int i = 0; i < n; i++)
232+
// Build Takagi U matrix with proper phase handling for negative eigenvalues
233+
// For Takagi: A = U * S * U^T, where S has non-negative singular values
234+
// If eigenvalue λ_i < 0, we need U column to have phase factor i (imaginary unit)
235+
var U = new Matrix<Complex<T>>(n, n);
236+
for (int j = 0; j < n; j++)
200237
{
201-
S[i, i] = NumOps.Sqrt(NumOps.Abs(A[i, i]));
238+
T eigenvalue = A[j, j];
239+
S[j, j] = NumOps.Abs(eigenvalue);
240+
241+
if (NumOps.GreaterThanOrEquals(eigenvalue, NumOps.Zero))
242+
{
243+
// Positive eigenvalue: U column is real
244+
for (int i = 0; i < n; i++)
245+
{
246+
U[i, j] = new Complex<T>(V[i, j], NumOps.Zero);
247+
}
248+
}
249+
else
250+
{
251+
// Negative eigenvalue: U column is purely imaginary (multiply by i)
252+
// This ensures: (i*v) * |λ| * (i*v)^T = -|λ| * v*v^T = λ * v*v^T
253+
for (int i = 0; i < n; i++)
254+
{
255+
U[i, j] = new Complex<T>(NumOps.Zero, V[i, j]);
256+
}
257+
}
202258
}
203259

204260
return (S, U);
@@ -358,7 +414,9 @@ private T CalculateMagnitude(Complex<T> complex)
358414
// VECTORIZED: Process each column of eigenvectors as a vector operation
359415
for (int i = 0; i < rows; i++)
360416
{
361-
S[i, i] = NumOps.Sqrt(NumOps.Abs(eigenValues[i]));
417+
// For symmetric matrices, singular values = |eigenvalues|
418+
// (not sqrt(|eigenvalues|) since singular values = sqrt(eigenvalues of A^T A) = sqrt(λ²) = |λ|)
419+
S[i, i] = NumOps.Abs(eigenValues[i]);
362420

363421
Vector<T> eigenVector = eigenVectors.GetColumn(i);
364422
Vector<Complex<T>> complexVector = new Vector<Complex<T>>(eigenVector.Select(val => new Complex<T>(val, NumOps.Zero)));
@@ -521,10 +579,36 @@ private T CalculateMagnitude(Complex<T> complex)
521579
/// </remarks>
522580
public override Matrix<T> Invert()
523581
{
582+
// For Takagi decomposition A = U * S * U^T:
583+
// A^(-1) = (U^T)^(-1) * S^(-1) * U^(-1) = conj(U) * S^(-1) * U^H
584+
int n = SigmaMatrix.Rows;
524585
var invSigma = SigmaMatrix.InvertDiagonalMatrix();
525-
var invU = UnitaryMatrix.InvertUnitaryMatrix();
526586
var invSigmaComplex = invSigma.ToComplexMatrix();
527-
var inv = invU.Multiply(invSigmaComplex).Multiply(invU.Transpose());
587+
588+
// Compute U^H (conjugate transpose of U)
589+
var uConjTranspose = new Matrix<Complex<T>>(n, n);
590+
for (int i = 0; i < n; i++)
591+
{
592+
for (int j = 0; j < n; j++)
593+
{
594+
var u = UnitaryMatrix[j, i];
595+
uConjTranspose[i, j] = new Complex<T>(u.Real, NumOps.Negate(u.Imaginary));
596+
}
597+
}
598+
599+
// Compute conj(U) (element-wise conjugate)
600+
var uConj = new Matrix<Complex<T>>(n, n);
601+
for (int i = 0; i < n; i++)
602+
{
603+
for (int j = 0; j < n; j++)
604+
{
605+
var u = UnitaryMatrix[i, j];
606+
uConj[i, j] = new Complex<T>(u.Real, NumOps.Negate(u.Imaginary));
607+
}
608+
}
609+
610+
// Compute conj(U) * S^(-1) * U^H
611+
var inv = uConj.Multiply(invSigmaComplex).Multiply(uConjTranspose);
528612

529613
return inv.ToRealMatrix();
530614
}
@@ -541,15 +625,69 @@ public override Matrix<T> Invert()
541625
/// but with matrices instead of simple numbers. Using the decomposition makes
542626
/// this process much more efficient than directly inverting the matrix.
543627
/// </para>
628+
/// <para>
629+
/// For Takagi decomposition A = U * S * U^T, solve x = conj(U) * S^(-1) * U^H * b
630+
/// where U^H is the conjugate transpose and conj(U) is element-wise conjugate.
631+
/// </para>
544632
/// </remarks>
545633
public override Vector<T> Solve(Vector<T> bVector)
546634
{
547-
// VECTORIZED: Convert to complex vector using Select
548-
var bComplex = new Vector<Complex<T>>(bVector.Select(val => new Complex<T>(val, NumOps.Zero)));
635+
int n = bVector.Length;
636+
637+
// For Takagi decomposition A = U * S * U^T, solve x = conj(U) * S^(-1) * U^H * b
638+
// Step 1: Compute y = U^H * b (conjugate transpose of U times b)
639+
var yVector = new Vector<Complex<T>>(n);
640+
for (int i = 0; i < n; i++)
641+
{
642+
Complex<T> sum = new Complex<T>(NumOps.Zero, NumOps.Zero);
643+
for (int j = 0; j < n; j++)
644+
{
645+
// U^H[i,j] = conj(U[j,i])
646+
var uConjugate = new Complex<T>(UnitaryMatrix[j, i].Real, NumOps.Negate(UnitaryMatrix[j, i].Imaginary));
647+
var bVal = new Complex<T>(bVector[j], NumOps.Zero);
648+
var product = new Complex<T>(
649+
NumOps.Subtract(NumOps.Multiply(uConjugate.Real, bVal.Real), NumOps.Multiply(uConjugate.Imaginary, bVal.Imaginary)),
650+
NumOps.Add(NumOps.Multiply(uConjugate.Real, bVal.Imaginary), NumOps.Multiply(uConjugate.Imaginary, bVal.Real)));
651+
sum = new Complex<T>(NumOps.Add(sum.Real, product.Real), NumOps.Add(sum.Imaginary, product.Imaginary));
652+
}
653+
yVector[i] = sum;
654+
}
655+
656+
// Step 2: Compute z = S^(-1) * y (divide by diagonal singular values)
657+
var zVector = new Vector<Complex<T>>(n);
658+
T tolerance = NumOps.FromDouble(1e-10);
659+
for (int i = 0; i < n; i++)
660+
{
661+
T sigma = SigmaMatrix[i, i];
662+
// Use tolerance-based check for near-zero singular values (pseudoinverse behavior)
663+
zVector[i] = NumOps.LessThan(NumOps.Abs(sigma), tolerance)
664+
? new Complex<T>(NumOps.Zero, NumOps.Zero)
665+
: new Complex<T>(
666+
NumOps.Divide(yVector[i].Real, sigma),
667+
NumOps.Divide(yVector[i].Imaginary, sigma));
668+
}
549669

550-
var yVector = UnitaryMatrix.ForwardSubstitution(bComplex);
670+
// Step 3: Compute x = conj(U) * z
671+
// Note: For Takagi A = U*S*U^T, we need (U^T)^(-1) = conj(U) for unitary U
672+
var xVector = new Vector<T>(n);
673+
for (int i = 0; i < n; i++)
674+
{
675+
Complex<T> sum = new Complex<T>(NumOps.Zero, NumOps.Zero);
676+
for (int j = 0; j < n; j++)
677+
{
678+
// Use conjugate of U[i,j]
679+
var u = UnitaryMatrix[i, j];
680+
var uConj = new Complex<T>(u.Real, NumOps.Negate(u.Imaginary));
681+
var z = zVector[j];
682+
var product = new Complex<T>(
683+
NumOps.Subtract(NumOps.Multiply(uConj.Real, z.Real), NumOps.Multiply(uConj.Imaginary, z.Imaginary)),
684+
NumOps.Add(NumOps.Multiply(uConj.Real, z.Imaginary), NumOps.Multiply(uConj.Imaginary, z.Real)));
685+
sum = new Complex<T>(NumOps.Add(sum.Real, product.Real), NumOps.Add(sum.Imaginary, product.Imaginary));
686+
}
687+
// Take real part of result
688+
xVector[i] = sum.Real;
689+
}
551690

552-
var result = SigmaMatrix.ToComplexMatrix().BackwardSubstitution(yVector);
553-
return result.ToRealVector();
691+
return xVector;
554692
}
555693
}

0 commit comments

Comments
 (0)