From f10239d53fc067358749f08cd17e54dbaa42b2af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?G=C3=B6ran=20B=C3=A4cklund?= Date: Mon, 27 Jul 2026 20:44:12 +0200 Subject: [PATCH 1/2] feat(numerics): add LU and Cholesky decompositions - LuDecomposition with partial pivoting (P*A = L*U): cached factorization with Solve for vector/list/matrix right-hand sides, Determinant, Inverse, IsSingular, Lower/Upper/PermutationMatrix - CholeskyDecomposition (A = L*L^T) with IsPositiveDefinite, Solve, Determinant, Inverse - Extension facade matrix.Lu() / matrix.Cholesky() - Matrix.Inverse now uses LU internally (O(n^3) instead of adjugate O(n!)), same public API and error messages - LinearSystemSolver overloads solve via LU substitution instead of computing the full inverse - 25 new unit tests; full suite green (1486 passed) - Check off LinearAlgebra roadmap Phase 1 Co-Authored-By: Claude Fable 5 --- .../NumericTest/CholeskyDecompositionTest.cs | 147 +++++++++ Numerics/NumericTest/LuDecompositionTest.cs | 213 ++++++++++++++ .../DifferentialEquationExtensions.cs | 27 +- .../Decompositions/CholeskyDecomposition.cs | 218 ++++++++++++++ .../Decompositions/LuDecomposition.cs | 278 ++++++++++++++++++ .../MatrixDecompositionExtensions.cs | 19 ++ Numerics/Numerics/Numerics/Objects/Matrix.cs | 9 +- docs/LinearAlgebraRoadMap.md | 10 +- 8 files changed, 899 insertions(+), 22 deletions(-) create mode 100644 Numerics/NumericTest/CholeskyDecompositionTest.cs create mode 100644 Numerics/NumericTest/LuDecompositionTest.cs create mode 100644 Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/CholeskyDecomposition.cs create mode 100644 Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs create mode 100644 Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs diff --git a/Numerics/NumericTest/CholeskyDecompositionTest.cs b/Numerics/NumericTest/CholeskyDecompositionTest.cs new file mode 100644 index 0000000..a885456 --- /dev/null +++ b/Numerics/NumericTest/CholeskyDecompositionTest.cs @@ -0,0 +1,147 @@ +using CSharpNumerics.Numerics.LinearAlgebra; +using CSharpNumerics.Numerics.Objects; + +namespace NumericsTests +{ + [TestClass] + public class CholeskyDecompositionTest + { + private const double Tolerance = 1e-10; + + private static readonly double[,] SpdValues = + { + { 4, 12, -16 }, + { 12, 37, -43 }, + { -16, -43, 98 } + }; + + private static void AssertMatricesEqual(Matrix expected, Matrix actual, double tolerance = Tolerance) + { + Assert.AreEqual(expected.rowLength, actual.rowLength); + Assert.AreEqual(expected.columnLength, actual.columnLength); + + for (var i = 0; i < expected.rowLength; i++) + { + for (var j = 0; j < expected.columnLength; j++) + { + Assert.IsTrue(Math.Abs(expected.values[i, j] - actual.values[i, j]) < tolerance, + $"Mismatch at ({i},{j}): expected {expected.values[i, j]}, actual {actual.values[i, j]}"); + } + } + } + + [TestMethod] + public void TestKnownFactorization() + { + var matrix = new Matrix(SpdValues); + var cholesky = matrix.Cholesky(); + + Assert.IsTrue(cholesky.IsPositiveDefinite); + + var expectedLower = new Matrix(new double[,] + { + { 2, 0, 0 }, + { 6, 1, 0 }, + { -8, 5, 3 } + }); + + AssertMatricesEqual(expectedLower, cholesky.Lower); + } + + [TestMethod] + public void TestRoundTrip() + { + var matrix = new Matrix(SpdValues); + var lower = matrix.Cholesky().Lower; + + AssertMatricesEqual(matrix, lower * lower.Transpose()); + } + + [TestMethod] + public void TestSolveKnownSystem() + { + var matrix = new Matrix(SpdValues); + var x = new VectorN(new double[] { 1, 2, 3 }); + var b = matrix * x; + + var solution = matrix.Cholesky().Solve(b); + + Assert.IsTrue(Math.Abs(solution[0] - 1) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution[2] - 3) < Tolerance); + } + + [TestMethod] + public void TestSolveList() + { + var matrix = new Matrix(SpdValues); + var b = matrix * new List { 1, 2, 3 }; + + var solution = matrix.Cholesky().Solve(b); + + Assert.IsTrue(Math.Abs(solution[0] - 1) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution[2] - 3) < Tolerance); + } + + [TestMethod] + public void TestSolveMatchesLu() + { + var matrix = new Matrix(SpdValues); + var b = new VectorN(new double[] { 1, -2, 5 }); + + var choleskySolution = matrix.Cholesky().Solve(b); + var luSolution = matrix.Lu().Solve(b); + + for (var i = 0; i < 3; i++) + { + Assert.IsTrue(Math.Abs(choleskySolution[i] - luSolution[i]) < Tolerance); + } + } + + [TestMethod] + public void TestDeterminant() + { + var matrix = new Matrix(SpdValues); + + Assert.IsTrue(Math.Abs(matrix.Cholesky().Determinant() - 36) < Tolerance); + } + + [TestMethod] + public void TestInverseRoundTrip() + { + var matrix = new Matrix(SpdValues); + var inverse = matrix.Cholesky().Inverse(); + + AssertMatricesEqual(new Matrix(matrix.identity), matrix * inverse, 1e-9); + } + + [TestMethod] + public void TestSymmetricIndefiniteIsRejected() + { + var matrix = new Matrix(new double[,] { { 1, 2 }, { 2, 1 } }); + var cholesky = matrix.Cholesky(); + + Assert.IsFalse(cholesky.IsPositiveDefinite); + Assert.ThrowsException(() => cholesky.Solve(new VectorN(new double[] { 1, 1 }))); + } + + [TestMethod] + public void TestNonSymmetricIsRejected() + { + var matrix = new Matrix(new double[,] { { 4, 1 }, { 2, 3 } }); + var cholesky = matrix.Cholesky(); + + Assert.IsFalse(cholesky.IsPositiveDefinite); + Assert.ThrowsException(() => cholesky.Lower); + } + + [TestMethod] + public void TestNonSquareMatrixThrows() + { + var matrix = new Matrix(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } }); + + Assert.ThrowsException(() => matrix.Cholesky()); + } + } +} diff --git a/Numerics/NumericTest/LuDecompositionTest.cs b/Numerics/NumericTest/LuDecompositionTest.cs new file mode 100644 index 0000000..c9f5054 --- /dev/null +++ b/Numerics/NumericTest/LuDecompositionTest.cs @@ -0,0 +1,213 @@ +using CSharpNumerics.Numerics; +using CSharpNumerics.Numerics.LinearAlgebra; +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; +using CSharpNumerics.Numerics.Objects; + +namespace NumericsTests +{ + [TestClass] + public class LuDecompositionTest + { + private const double Tolerance = 1e-10; + + private static void AssertMatricesEqual(Matrix expected, Matrix actual, double tolerance = Tolerance) + { + Assert.AreEqual(expected.rowLength, actual.rowLength); + Assert.AreEqual(expected.columnLength, actual.columnLength); + + for (var i = 0; i < expected.rowLength; i++) + { + for (var j = 0; j < expected.columnLength; j++) + { + Assert.IsTrue(Math.Abs(expected.values[i, j] - actual.values[i, j]) < tolerance, + $"Mismatch at ({i},{j}): expected {expected.values[i, j]}, actual {actual.values[i, j]}"); + } + } + } + + [TestMethod] + public void TestKnownFactorization() + { + var matrix = new Matrix(new double[,] { { 4, 3 }, { 6, 3 } }); + var lu = matrix.Lu(); + + var lower = lu.Lower; + var upper = lu.Upper; + + Assert.IsTrue(Math.Abs(lower.values[0, 0] - 1) < Tolerance); + Assert.IsTrue(Math.Abs(lower.values[1, 0] - 2.0 / 3.0) < Tolerance); + Assert.IsTrue(Math.Abs(lower.values[1, 1] - 1) < Tolerance); + + Assert.IsTrue(Math.Abs(upper.values[0, 0] - 6) < Tolerance); + Assert.IsTrue(Math.Abs(upper.values[0, 1] - 3) < Tolerance); + Assert.IsTrue(Math.Abs(upper.values[1, 1] - 1) < Tolerance); + } + + [TestMethod] + public void TestRoundTrip() + { + var matrix = new Matrix(new double[,] + { + { 2, 1, 1, 0 }, + { 4, 3, 3, 1 }, + { 8, 7, 9, 5 }, + { 6, 7, 9, 8 } + }); + + var lu = matrix.Lu(); + AssertMatricesEqual(lu.PermutationMatrix * matrix, lu.Lower * lu.Upper); + } + + [TestMethod] + public void TestLowerIsUnitLowerAndUpperIsUpper() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + var lu = matrix.Lu(); + + for (var i = 0; i < 3; i++) + { + Assert.IsTrue(lu.Lower.values[i, i] == 1); + for (var j = i + 1; j < 3; j++) + { + Assert.IsTrue(lu.Lower.values[i, j] == 0); + Assert.IsTrue(lu.Upper.values[j, i] == 0); + } + } + } + + [TestMethod] + public void TestSolveKnownSystem() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var solution = matrix.Lu().Solve(new VectorN(new double[] { 8, -11, -3 })); + + Assert.IsTrue(Math.Abs(solution[0] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 3) < Tolerance); + Assert.IsTrue(Math.Abs(solution[2] - (-1)) < Tolerance); + } + + [TestMethod] + public void TestSolveList() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var solution = matrix.Lu().Solve(new List { 8, -11, -3 }); + + Assert.IsTrue(Math.Abs(solution[0] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 3) < Tolerance); + Assert.IsTrue(Math.Abs(solution[2] - (-1)) < Tolerance); + } + + [TestMethod] + public void TestSolveReusesFactorization() + { + var matrix = new Matrix(new double[,] { { 4, 3 }, { 6, 3 } }); + var lu = matrix.Lu(); + + var x1 = lu.Solve(new VectorN(new double[] { 10, 12 })); + var x2 = lu.Solve(new VectorN(new double[] { 7, 9 })); + + var b1 = matrix * x1; + var b2 = matrix * x2; + + Assert.IsTrue(Math.Abs(b1[0] - 10) < Tolerance); + Assert.IsTrue(Math.Abs(b1[1] - 12) < Tolerance); + Assert.IsTrue(Math.Abs(b2[0] - 7) < Tolerance); + Assert.IsTrue(Math.Abs(b2[1] - 9) < Tolerance); + } + + [TestMethod] + public void TestSolveMultipleRightHandSides() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var rhs = new Matrix(new double[,] { { 8, 1 }, { -11, 0 }, { -3, 2 } }); + + var x = matrix.Lu().Solve(rhs); + + AssertMatricesEqual(rhs, matrix * x); + } + + [TestMethod] + public void TestDeterminantMatchesCofactorExpansion() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + + Assert.IsTrue(Math.Abs(matrix.Lu().Determinant() - 77) < Tolerance); + Assert.IsTrue(Math.Abs(matrix.Lu().Determinant() - matrix.Determinant()) < Tolerance); + } + + [TestMethod] + public void TestInverseRoundTrip() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + var inverse = matrix.Lu().Inverse(); + + AssertMatricesEqual(new Matrix(matrix.identity), matrix * inverse); + AssertMatricesEqual(new Matrix(matrix.identity), inverse * matrix); + } + + [TestMethod] + public void TestMatrixInverseUsesLuAndMatchesAdjugate() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + var inverse = matrix.Inverse(); + var expected = matrix.Adjugate() / matrix.Determinant(); + + AssertMatricesEqual(expected, inverse); + } + + [TestMethod] + public void TestSingularMatrixIsDetected() + { + var matrix = new Matrix(new double[,] { { 1, 2 }, { 2, 4 } }); + var lu = matrix.Lu(); + + Assert.IsTrue(lu.IsSingular); + Assert.ThrowsException(() => lu.Solve(new VectorN(new double[] { 1, 1 }))); + Assert.ThrowsException(() => matrix.Inverse()); + } + + [TestMethod] + public void TestNonSquareMatrixThrows() + { + var matrix = new Matrix(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } }); + + Assert.ThrowsException(() => matrix.Lu()); + } + + [TestMethod] + public void TestPivotingHandlesZeroLeadingElement() + { + var matrix = new Matrix(new double[,] { { 0, 1 }, { 1, 0 } }); + var lu = matrix.Lu(); + + Assert.IsFalse(lu.IsSingular); + + var solution = lu.Solve(new VectorN(new double[] { 3, 5 })); + + Assert.IsTrue(Math.Abs(solution[0] - 5) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 3) < Tolerance); + } + + [TestMethod] + public void TestLinearSystemSolverStillWorks() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var solution = matrix.LinearSystemSolver(new VectorN(new double[] { 8, -11, -3 })); + + Assert.IsTrue(Math.Abs(solution[0] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - 3) < Tolerance); + Assert.IsTrue(Math.Abs(solution[2] - (-1)) < Tolerance); + } + + [TestMethod] + public void TestLinearSystemSolverVector() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var solution = matrix.LinearSystemSolver(new Vector(8, -11, -3)); + + Assert.IsTrue(Math.Abs(solution.x - 2) < Tolerance); + Assert.IsTrue(Math.Abs(solution.y - 3) < Tolerance); + Assert.IsTrue(Math.Abs(solution.z - (-1)) < Tolerance); + } + } +} diff --git a/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs b/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs index 9dcd812..a07ef38 100644 --- a/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs +++ b/Numerics/Numerics/Numerics/DifferentialEquationExtensions.cs @@ -1,4 +1,5 @@ -using CSharpNumerics.Numerics.Objects; +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; +using CSharpNumerics.Numerics.Objects; using System.Collections.Generic; using System.Linq; using System; @@ -526,42 +527,42 @@ public static double[] VelocityVerlet( #endregion /// - /// Solves a linear system A x = b by computing x = A^{-1} b. + /// Solves a linear system A x = b via LU decomposition with partial pivoting. /// - /// Coefficient matrix A. + /// Coefficient matrix A (2x2 or 3x3). /// Right-hand side vector b. /// Solution vector x. public static Vector LinearSystemSolver(this Matrix matrix, Vector vector) { - var values = matrix.Inverse() * vector; + var b = matrix.columnLength == 2 + ? new List { vector.x, vector.y } + : new List { vector.x, vector.y, vector.z }; + + var x = new LuDecomposition(matrix).Solve(b); - return values; + return new Vector(x[0], x[1], x.Count > 2 ? x[2] : 0); } /// - /// Solves a linear system A x = b by computing x = A^{-1} b. + /// Solves a linear system A x = b via LU decomposition with partial pivoting. /// /// Coefficient matrix A. /// Right-hand side vector b. /// Solution vector x. public static VectorN LinearSystemSolver(this Matrix matrix, VectorN vector) { - var values = matrix.Inverse() * vector; - - return values; + return new LuDecomposition(matrix).Solve(vector); } /// - /// Solves a linear system A x = b by computing x = A^{-1} b. + /// Solves a linear system A x = b via LU decomposition with partial pivoting. /// /// Coefficient matrix A. /// Right-hand side vector b. /// Solution vector x as a list. public static List LinearSystemSolver(this Matrix matrix, List vector) { - var values = matrix.Inverse() * vector; - - return values; + return new LuDecomposition(matrix).Solve(vector); } /// diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/CholeskyDecomposition.cs b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/CholeskyDecomposition.cs new file mode 100644 index 0000000..7c24774 --- /dev/null +++ b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/CholeskyDecomposition.cs @@ -0,0 +1,218 @@ +using System; +using System.Collections.Generic; +using CSharpNumerics.Numerics.Objects; + +namespace CSharpNumerics.Numerics.LinearAlgebra.Decompositions; + +/// +/// Cholesky decomposition A = L·Lᵀ for symmetric positive definite matrices, +/// where L is lower triangular. Roughly twice as fast as LU for such systems. +/// The factorization is computed once and can be reused for multiple solves. +/// +public sealed class CholeskyDecomposition +{ + private const double SymmetryTolerance = 1e-10; + + private readonly double[,] l; + private readonly int n; + + public CholeskyDecomposition(Matrix matrix) + { + if (matrix.rowLength != matrix.columnLength) + { + throw new Exception("Is not a NxN matrix"); + } + + n = matrix.rowLength; + l = new double[n, n]; + IsPositiveDefinite = IsSymmetric(matrix); + + for (var j = 0; j < n; j++) + { + var d = matrix.values[j, j]; + for (var k = 0; k < j; k++) + { + d -= l[j, k] * l[j, k]; + } + + if (d <= 0.0) + { + IsPositiveDefinite = false; + return; + } + + l[j, j] = Math.Sqrt(d); + + for (var i = j + 1; i < n; i++) + { + var sum = matrix.values[i, j]; + for (var k = 0; k < j; k++) + { + sum -= l[i, k] * l[j, k]; + } + l[i, j] = sum / l[j, j]; + } + } + } + + /// + /// True if the matrix is symmetric and all pivots are positive, i.e. the factorization A = L·Lᵀ exists. + /// + public bool IsPositiveDefinite { get; } + + /// + /// The lower triangular factor L. + /// + public Matrix Lower + { + get + { + EnsurePositiveDefinite(); + return new Matrix((double[,])l.Clone()); + } + } + + /// + /// Computes the determinant from the factorization as ∏ L[i,i]². + /// + public double Determinant() + { + EnsurePositiveDefinite(); + var determinant = 1.0; + for (var j = 0; j < n; j++) + { + determinant *= l[j, j] * l[j, j]; + } + return determinant; + } + + /// + /// Solves A x = b using the cached factorization (forward and back substitution). + /// + public VectorN Solve(VectorN b) + { + if (b.Length != n) + { + throw new Exception("The vector length must match the matrix dimension"); + } + + var x = new double[n]; + for (var i = 0; i < n; i++) + { + x[i] = b[i]; + } + + SolveInPlace(x); + return new VectorN(x); + } + + /// + /// Solves A x = b using the cached factorization. + /// + public List Solve(List b) + { + if (b.Count != n) + { + throw new Exception("The vector length must match the matrix dimension"); + } + + var x = b.ToArray(); + SolveInPlace(x); + return new List(x); + } + + /// + /// Solves A X = B for multiple right-hand sides (one per column of B). + /// + public Matrix Solve(Matrix b) + { + if (b.rowLength != n) + { + throw new Exception("The row length of B must match the matrix dimension"); + } + + var columns = b.columnLength; + var x = new double[n, columns]; + + for (var c = 0; c < columns; c++) + { + var column = new double[n]; + for (var i = 0; i < n; i++) + { + column[i] = b.values[i, c]; + } + + SolveInPlace(column); + + for (var i = 0; i < n; i++) + { + x[i, c] = column[i]; + } + } + + return new Matrix(x); + } + + /// + /// Computes the inverse by solving A X = I. Reuses the cached factorization. + /// + public Matrix Inverse() + { + var identity = new double[n, n]; + for (var i = 0; i < n; i++) + { + identity[i, i] = 1.0; + } + return Solve(new Matrix(identity)); + } + + private void SolveInPlace(double[] x) + { + EnsurePositiveDefinite(); + + for (var i = 0; i < n; i++) + { + for (var k = 0; k < i; k++) + { + x[i] -= l[i, k] * x[k]; + } + x[i] /= l[i, i]; + } + + for (var i = n - 1; i >= 0; i--) + { + for (var k = i + 1; k < n; k++) + { + x[i] -= l[k, i] * x[k]; + } + x[i] /= l[i, i]; + } + } + + private void EnsurePositiveDefinite() + { + if (!IsPositiveDefinite) + { + throw new Exception("The matrix is not symmetric positive definite"); + } + } + + private static bool IsSymmetric(Matrix matrix) + { + for (var i = 0; i < matrix.rowLength; i++) + { + for (var j = i + 1; j < matrix.columnLength; j++) + { + var a = matrix.values[i, j]; + var b = matrix.values[j, i]; + var scale = Math.Max(1.0, Math.Max(Math.Abs(a), Math.Abs(b))); + + if (Math.Abs(a - b) > SymmetryTolerance * scale) + { + return false; + } + } + } + return true; + } +} diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs new file mode 100644 index 0000000..0cce65f --- /dev/null +++ b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/LuDecomposition.cs @@ -0,0 +1,278 @@ +using System; +using System.Collections.Generic; +using CSharpNumerics.Numerics.Objects; + +namespace CSharpNumerics.Numerics.LinearAlgebra.Decompositions; + +/// +/// LU decomposition with partial (row) pivoting: P·A = L·U where L is unit lower triangular, +/// U is upper triangular and P is a permutation matrix. +/// The factorization is computed once and can be reused for multiple solves. +/// +public sealed class LuDecomposition +{ + private readonly double[,] lu; + private readonly int[] pivot; + private readonly int pivotSign; + private readonly int n; + + public LuDecomposition(Matrix matrix) + { + if (matrix.rowLength != matrix.columnLength) + { + throw new Exception("Is not a NxN matrix"); + } + + n = matrix.rowLength; + lu = (double[,])matrix.values.Clone(); + pivot = new int[n]; + pivotSign = 1; + + for (var i = 0; i < n; i++) + { + pivot[i] = i; + } + + for (var j = 0; j < n; j++) + { + for (var i = 0; i < n; i++) + { + var kMax = Math.Min(i, j); + var sum = 0.0; + for (var k = 0; k < kMax; k++) + { + sum += lu[i, k] * lu[k, j]; + } + lu[i, j] -= sum; + } + + var p = j; + for (var i = j + 1; i < n; i++) + { + if (Math.Abs(lu[i, j]) > Math.Abs(lu[p, j])) + { + p = i; + } + } + + if (p != j) + { + for (var k = 0; k < n; k++) + { + (lu[p, k], lu[j, k]) = (lu[j, k], lu[p, k]); + } + + (pivot[p], pivot[j]) = (pivot[j], pivot[p]); + pivotSign = -pivotSign; + } + + if (lu[j, j] != 0.0) + { + for (var i = j + 1; i < n; i++) + { + lu[i, j] /= lu[j, j]; + } + } + } + } + + /// + /// True if the matrix is singular (a zero pivot was encountered) and cannot be solved. + /// + public bool IsSingular + { + get + { + for (var j = 0; j < n; j++) + { + if (lu[j, j] == 0.0) + { + return true; + } + } + return false; + } + } + + /// + /// The unit lower triangular factor L. + /// + public Matrix Lower + { + get + { + var l = new double[n, n]; + for (var i = 0; i < n; i++) + { + for (var j = 0; j < n; j++) + { + l[i, j] = i > j ? lu[i, j] : (i == j ? 1.0 : 0.0); + } + } + return new Matrix(l); + } + } + + /// + /// The upper triangular factor U. + /// + public Matrix Upper + { + get + { + var u = new double[n, n]; + for (var i = 0; i < n; i++) + { + for (var j = i; j < n; j++) + { + u[i, j] = lu[i, j]; + } + } + return new Matrix(u); + } + } + + /// + /// The row permutation applied by pivoting: row i of P·A is row Permutation[i] of A. + /// + public int[] Permutation => (int[])pivot.Clone(); + + /// + /// The permutation matrix P such that P·A = L·U. + /// + public Matrix PermutationMatrix + { + get + { + var p = new double[n, n]; + for (var i = 0; i < n; i++) + { + p[i, pivot[i]] = 1.0; + } + return new Matrix(p); + } + } + + /// + /// Computes the determinant from the factorization as sign(P) · ∏ U[i,i]. + /// + public double Determinant() + { + var determinant = (double)pivotSign; + for (var j = 0; j < n; j++) + { + determinant *= lu[j, j]; + } + return determinant; + } + + /// + /// Solves A x = b using the cached factorization (forward and back substitution). + /// + public VectorN Solve(VectorN b) + { + if (b.Length != n) + { + throw new Exception("The vector length must match the matrix dimension"); + } + + var x = new double[n]; + for (var i = 0; i < n; i++) + { + x[i] = b[pivot[i]]; + } + + SolveInPlace(x); + return new VectorN(x); + } + + /// + /// Solves A x = b using the cached factorization. + /// + public List Solve(List b) + { + if (b.Count != n) + { + throw new Exception("The vector length must match the matrix dimension"); + } + + var x = new double[n]; + for (var i = 0; i < n; i++) + { + x[i] = b[pivot[i]]; + } + + SolveInPlace(x); + return new List(x); + } + + /// + /// Solves A X = B for multiple right-hand sides (one per column of B). + /// + public Matrix Solve(Matrix b) + { + if (b.rowLength != n) + { + throw new Exception("The row length of B must match the matrix dimension"); + } + + var columns = b.columnLength; + var x = new double[n, columns]; + + for (var c = 0; c < columns; c++) + { + var column = new double[n]; + for (var i = 0; i < n; i++) + { + column[i] = b.values[pivot[i], c]; + } + + SolveInPlace(column); + + for (var i = 0; i < n; i++) + { + x[i, c] = column[i]; + } + } + + return new Matrix(x); + } + + /// + /// Computes the inverse by solving A X = I. Reuses the cached factorization. + /// + public Matrix Inverse() + { + var identity = new double[n, n]; + for (var i = 0; i < n; i++) + { + identity[i, i] = 1.0; + } + return Solve(new Matrix(identity)); + } + + private void SolveInPlace(double[] x) + { + if (IsSingular) + { + throw new Exception("This matrix is not invertible"); + } + + for (var i = 1; i < n; i++) + { + for (var k = 0; k < i; k++) + { + x[i] -= lu[i, k] * x[k]; + } + } + + for (var i = n - 1; i >= 0; i--) + { + for (var k = i + 1; k < n; k++) + { + x[i] -= lu[i, k] * x[k]; + } + x[i] /= lu[i, i]; + } + } +} diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs new file mode 100644 index 0000000..5689667 --- /dev/null +++ b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs @@ -0,0 +1,19 @@ +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; +using CSharpNumerics.Numerics.Objects; + +namespace CSharpNumerics.Numerics.LinearAlgebra; + +public static class MatrixDecompositionExtensions +{ + /// + /// Computes the LU decomposition with partial pivoting: P·A = L·U. + /// The returned object caches the factorization for reuse across multiple solves. + /// + public static LuDecomposition Lu(this Matrix matrix) => new LuDecomposition(matrix); + + /// + /// Computes the Cholesky decomposition A = L·Lᵀ for a symmetric positive definite matrix. + /// The returned object caches the factorization for reuse across multiple solves. + /// + public static CholeskyDecomposition Cholesky(this Matrix matrix) => new CholeskyDecomposition(matrix); +} diff --git a/Numerics/Numerics/Numerics/Objects/Matrix.cs b/Numerics/Numerics/Numerics/Objects/Matrix.cs index 83ad3b5..76dea01 100644 --- a/Numerics/Numerics/Numerics/Objects/Matrix.cs +++ b/Numerics/Numerics/Numerics/Objects/Matrix.cs @@ -1,5 +1,6 @@ using System; using System.Collections.Generic; +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; namespace CSharpNumerics.Numerics.Objects; @@ -54,12 +55,12 @@ public Matrix(int rows, int cols) } public Matrix Inverse() { - var determinant =Determinant(); - if (determinant == 0) { + var lu = new LuDecomposition(this); + if (lu.IsSingular) + { throw new Exception("This matrix is not invertible"); } - var adj = Adjugate(); - return adj / determinant; + return lu.Inverse(); } diff --git a/docs/LinearAlgebraRoadMap.md b/docs/LinearAlgebraRoadMap.md index 03fc653..bd7165c 100644 --- a/docs/LinearAlgebraRoadMap.md +++ b/docs/LinearAlgebraRoadMap.md @@ -105,11 +105,11 @@ Placeras i `Numerics/RootFinding/`. Direkt användbart i: `KeplerOrbit` (Keplers ## Implementationsplan — Faser ### Phase 1 — Dekompositionsgrund -- [ ] Skapa `Numerics/LinearAlgebra/Decompositions/`-struktur -- [ ] Implementera `LuDecomposition` med partiell pivotering + `Solve`/`Determinant`/`Inverse` -- [ ] Implementera `CholeskyDecomposition` + `IsPositiveDefinite` -- [ ] Refaktorera `Matrix.Inverse` och `LinearSystemSolver` till LU internt (inga API-ändringar) -- [ ] Enhetstester: kända faktoriseringar, singulära matriser, round-trip `A ≈ P·L·U` +- [x] Skapa `Numerics/LinearAlgebra/Decompositions/`-struktur +- [x] Implementera `LuDecomposition` med partiell pivotering + `Solve`/`Determinant`/`Inverse` +- [x] Implementera `CholeskyDecomposition` + `IsPositiveDefinite` +- [x] Refaktorera `Matrix.Inverse` och `LinearSystemSolver` till LU internt (inga API-ändringar) +- [x] Enhetstester: kända faktoriseringar, singulära matriser, round-trip `A ≈ P·L·U` ### Phase 2 — QR & egendekomposition - [ ] Implementera `QrDecomposition` (Householder) + minsta kvadrat-`Solve` From 8f627dc8cbb09bb3c98cb052041885df4b867bed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?G=C3=B6ran=20B=C3=A4cklund?= Date: Fri, 7 Aug 2026 16:19:33 +0200 Subject: [PATCH 2/2] feat(numerics): add QR and eigenvalue decompositions - QrDecomposition via Householder reflections: economy Q/R, IsFullRank with relative tolerance, least squares Solve for overdetermined systems - EigenDecomposition (EISPACK/JAMA port): symmetric path via Householder tridiagonalization + implicit QL with shifts (ascending eigenvalues, orthonormal eigenvectors); non-symmetric path via Hessenberg reduction + shifted QR iteration with complex conjugate pair support (A*V = V*D with block diagonal D) - Extension facade matrix.Qr() / matrix.Eigen() - Quantum module migrated: SymmetricEigenSolver (Jacobi) removed, SchrodingerExtensions now uses the general EigenDecomposition - 24 new unit tests: Q orthonormality, eigenpair residuals, least squares vs normal equations, complex pairs, analytic Laplacian eigenvalues - Check off LinearAlgebra roadmap Phase 2 Co-Authored-By: Claude Fable 5 --- .../NumericTest/EigenDecompositionTest.cs | 211 ++++ Numerics/NumericTest/QrDecompositionTest.cs | 139 +++ .../Decompositions/EigenDecomposition.cs | 939 ++++++++++++++++++ .../Decompositions/QrDecomposition.cs | 283 ++++++ .../MatrixDecompositionExtensions.cs | 12 + .../Physics/Quantum/SchrodingerExtensions.cs | 13 +- .../Physics/Quantum/SymmetricEigenSolver.cs | 102 -- docs/LinearAlgebraRoadMap.md | 10 +- 8 files changed, 1600 insertions(+), 109 deletions(-) create mode 100644 Numerics/NumericTest/EigenDecompositionTest.cs create mode 100644 Numerics/NumericTest/QrDecompositionTest.cs create mode 100644 Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/EigenDecomposition.cs create mode 100644 Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/QrDecomposition.cs delete mode 100644 Numerics/Numerics/Physics/Quantum/SymmetricEigenSolver.cs diff --git a/Numerics/NumericTest/EigenDecompositionTest.cs b/Numerics/NumericTest/EigenDecompositionTest.cs new file mode 100644 index 0000000..e6074b6 --- /dev/null +++ b/Numerics/NumericTest/EigenDecompositionTest.cs @@ -0,0 +1,211 @@ +using CSharpNumerics.Numerics.LinearAlgebra; +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; +using CSharpNumerics.Numerics.Objects; + +namespace NumericsTests +{ + [TestClass] + public class EigenDecompositionTest + { + private const double Tolerance = 1e-9; + + private static void AssertMatricesEqual(Matrix expected, Matrix actual, double tolerance = Tolerance) + { + Assert.AreEqual(expected.rowLength, actual.rowLength); + Assert.AreEqual(expected.columnLength, actual.columnLength); + + for (var i = 0; i < expected.rowLength; i++) + { + for (var j = 0; j < expected.columnLength; j++) + { + Assert.IsTrue(Math.Abs(expected.values[i, j] - actual.values[i, j]) < tolerance, + $"Mismatch at ({i},{j}): expected {expected.values[i, j]}, actual {actual.values[i, j]}"); + } + } + } + + [TestMethod] + public void TestSymmetricKnownEigenvalues() + { + var matrix = new Matrix(new double[,] { { 2, 1 }, { 1, 2 } }); + var eigen = matrix.Eigen(); + + Assert.IsTrue(eigen.IsSymmetric); + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[0] - 1) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[1] - 3) < Tolerance); + Assert.IsTrue(eigen.ImaginaryEigenvalues[0] == 0); + Assert.IsTrue(eigen.ImaginaryEigenvalues[1] == 0); + } + + [TestMethod] + public void TestSymmetricEigenvaluesAreAscending() + { + var matrix = new Matrix(new double[,] + { + { 4, 1, -2, 2 }, + { 1, 2, 0, 1 }, + { -2, 0, 3, -2 }, + { 2, 1, -2, -1 } + }); + var eigen = matrix.Eigen(); + + for (var i = 1; i < 4; i++) + { + Assert.IsTrue(eigen.RealEigenvalues[i] >= eigen.RealEigenvalues[i - 1]); + } + } + + [TestMethod] + public void TestSymmetricEigenpairResiduals() + { + var matrix = new Matrix(new double[,] + { + { 4, 1, -2, 2 }, + { 1, 2, 0, 1 }, + { -2, 0, 3, -2 }, + { 2, 1, -2, -1 } + }); + var eigen = matrix.Eigen(); + var v = eigen.EigenVectors; + + for (var k = 0; k < 4; k++) + { + var vector = v.ColumnSlice(k); + var av = matrix * vector; + + for (var i = 0; i < 4; i++) + { + Assert.IsTrue(Math.Abs(av[i] - eigen.RealEigenvalues[k] * vector[i]) < Tolerance, + $"Residual too large for eigenpair {k}, component {i}"); + } + } + } + + [TestMethod] + public void TestSymmetricEigenvectorsAreOrthonormal() + { + var matrix = new Matrix(new double[,] + { + { 4, 1, -2, 2 }, + { 1, 2, 0, 1 }, + { -2, 0, 3, -2 }, + { 2, 1, -2, -1 } + }); + var v = matrix.Eigen().EigenVectors; + var identity = v.Transpose() * v; + + for (var i = 0; i < 4; i++) + { + for (var j = 0; j < 4; j++) + { + Assert.IsTrue(Math.Abs(identity.values[i, j] - (i == j ? 1 : 0)) < Tolerance); + } + } + } + + [TestMethod] + public void TestSymmetricRoundTrip() + { + var matrix = new Matrix(new double[,] { { 4, 12, -16 }, { 12, 37, -43 }, { -16, -43, 98 } }); + var eigen = matrix.Eigen(); + + AssertMatricesEqual(matrix * eigen.EigenVectors, eigen.EigenVectors * eigen.DiagonalMatrix, 1e-8); + } + + [TestMethod] + public void TestNonSymmetricRealEigenvalues() + { + var matrix = new Matrix(new double[,] { { 4, 1 }, { 2, 3 } }); + var eigen = matrix.Eigen(); + + Assert.IsFalse(eigen.IsSymmetric); + + var values = new List(eigen.RealEigenvalues); + values.Sort(); + + Assert.IsTrue(Math.Abs(values[0] - 2) < Tolerance); + Assert.IsTrue(Math.Abs(values[1] - 5) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.ImaginaryEigenvalues[0]) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.ImaginaryEigenvalues[1]) < Tolerance); + } + + [TestMethod] + public void TestNonSymmetricRoundTrip() + { + var matrix = new Matrix(new double[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 10 } }); + var eigen = matrix.Eigen(); + + AssertMatricesEqual(matrix * eigen.EigenVectors, eigen.EigenVectors * eigen.DiagonalMatrix, 1e-8); + } + + [TestMethod] + public void TestComplexConjugatePair() + { + var matrix = new Matrix(new double[,] { { 0, -1 }, { 1, 0 } }); + var eigen = matrix.Eigen(); + + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[0]) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[1]) < Tolerance); + Assert.IsTrue(Math.Abs(Math.Abs(eigen.ImaginaryEigenvalues[0]) - 1) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.ImaginaryEigenvalues[0] + eigen.ImaginaryEigenvalues[1]) < Tolerance); + + AssertMatricesEqual(matrix * eigen.EigenVectors, eigen.EigenVectors * eigen.DiagonalMatrix); + } + + [TestMethod] + public void TestDeterminantMatchesEigenvalueProduct() + { + var matrix = new Matrix(new double[,] { { 4, 12, -16 }, { 12, 37, -43 }, { -16, -43, 98 } }); + var eigen = matrix.Eigen(); + + var product = 1.0; + foreach (var value in eigen.RealEigenvalues) + { + product *= value; + } + + Assert.IsTrue(Math.Abs(product - 36) < 1e-6); + } + + [TestMethod] + public void TestDiagonalMatrixIsTrivial() + { + var matrix = new Matrix(new double[,] { { 3, 0 }, { 0, 7 } }); + var eigen = matrix.Eigen(); + + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[0] - 3) < Tolerance); + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[1] - 7) < Tolerance); + } + + [TestMethod] + public void TestNonSquareMatrixThrows() + { + var matrix = new Matrix(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } }); + + Assert.ThrowsException(() => matrix.Eigen()); + } + + [TestMethod] + public void TestLargeSymmetricTridiagonal() + { + // Discrete 1D Laplacian: eigenvalues 2 - 2cos(k*pi/(n+1)) are known analytically. + const int n = 20; + var values = new double[n, n]; + for (var i = 0; i < n; i++) + { + values[i, i] = 2; + if (i > 0) values[i, i - 1] = -1; + if (i < n - 1) values[i, i + 1] = -1; + } + + var eigen = new EigenDecomposition(new Matrix(values)); + + for (var k = 0; k < n; k++) + { + var expected = 2 - 2 * Math.Cos((k + 1) * Math.PI / (n + 1)); + Assert.IsTrue(Math.Abs(eigen.RealEigenvalues[k] - expected) < 1e-8, + $"Eigenvalue {k}: expected {expected}, actual {eigen.RealEigenvalues[k]}"); + } + } + } +} diff --git a/Numerics/NumericTest/QrDecompositionTest.cs b/Numerics/NumericTest/QrDecompositionTest.cs new file mode 100644 index 0000000..3d940e8 --- /dev/null +++ b/Numerics/NumericTest/QrDecompositionTest.cs @@ -0,0 +1,139 @@ +using CSharpNumerics.Numerics.LinearAlgebra; +using CSharpNumerics.Numerics.Objects; + +namespace NumericsTests +{ + [TestClass] + public class QrDecompositionTest + { + private const double Tolerance = 1e-10; + + private static void AssertMatricesEqual(Matrix expected, Matrix actual, double tolerance = Tolerance) + { + Assert.AreEqual(expected.rowLength, actual.rowLength); + Assert.AreEqual(expected.columnLength, actual.columnLength); + + for (var i = 0; i < expected.rowLength; i++) + { + for (var j = 0; j < expected.columnLength; j++) + { + Assert.IsTrue(Math.Abs(expected.values[i, j] - actual.values[i, j]) < tolerance, + $"Mismatch at ({i},{j}): expected {expected.values[i, j]}, actual {actual.values[i, j]}"); + } + } + } + + [TestMethod] + public void TestRoundTripSquare() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + var qr = matrix.Qr(); + + AssertMatricesEqual(matrix, qr.Q * qr.R); + } + + [TestMethod] + public void TestRoundTripOverdetermined() + { + var matrix = new Matrix(new double[,] { { 1, 1 }, { 1, 2 }, { 1, 3 }, { 1, 4 } }); + var qr = matrix.Qr(); + + AssertMatricesEqual(matrix, qr.Q * qr.R); + } + + [TestMethod] + public void TestQHasOrthonormalColumns() + { + var matrix = new Matrix(new double[,] { { 1, 1 }, { 1, 2 }, { 1, 3 }, { 1, 4 } }); + var q = matrix.Qr().Q; + var identity = q.Transpose() * q; + + for (var i = 0; i < 2; i++) + { + for (var j = 0; j < 2; j++) + { + Assert.IsTrue(Math.Abs(identity.values[i, j] - (i == j ? 1 : 0)) < Tolerance); + } + } + } + + [TestMethod] + public void TestRIsUpperTriangular() + { + var matrix = new Matrix(new double[,] { { 1, 5, 2 }, { 0, 3, 7 }, { 2, -1, 4 } }); + var r = matrix.Qr().R; + + for (var i = 1; i < 3; i++) + { + for (var j = 0; j < i; j++) + { + Assert.IsTrue(r.values[i, j] == 0); + } + } + } + + [TestMethod] + public void TestLeastSquaresSolution() + { + var matrix = new Matrix(new double[,] { { 1, 1 }, { 1, 2 }, { 1, 3 } }); + var solution = matrix.Qr().Solve(new VectorN(new double[] { 6, 0, 0 })); + + Assert.IsTrue(Math.Abs(solution[0] - 8) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - (-3)) < Tolerance); + } + + [TestMethod] + public void TestSquareSolveMatchesLu() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var b = new VectorN(new double[] { 8, -11, -3 }); + + var qrSolution = matrix.Qr().Solve(b); + var luSolution = matrix.Lu().Solve(b); + + for (var i = 0; i < 3; i++) + { + Assert.IsTrue(Math.Abs(qrSolution[i] - luSolution[i]) < Tolerance); + } + } + + [TestMethod] + public void TestSolveList() + { + var matrix = new Matrix(new double[,] { { 1, 1 }, { 1, 2 }, { 1, 3 } }); + var solution = matrix.Qr().Solve(new List { 6, 0, 0 }); + + Assert.IsTrue(Math.Abs(solution[0] - 8) < Tolerance); + Assert.IsTrue(Math.Abs(solution[1] - (-3)) < Tolerance); + } + + [TestMethod] + public void TestSolveMultipleRightHandSides() + { + var matrix = new Matrix(new double[,] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } }); + var rhs = new Matrix(new double[,] { { 8, 1 }, { -11, 0 }, { -3, 2 } }); + + var x = matrix.Qr().Solve(rhs); + + AssertMatricesEqual(rhs, matrix * x, 1e-9); + } + + [TestMethod] + public void TestRankDeficientIsDetected() + { + var matrix = new Matrix(new double[,] { { 1, 2 }, { 2, 4 }, { 3, 6 } }); + var qr = matrix.Qr(); + + Assert.IsFalse(qr.IsFullRank); + Assert.ThrowsException(() => qr.Solve(new VectorN(new double[] { 1, 1, 1 }))); + } + + [TestMethod] + public void TestUnderdeterminedThrows() + { + var matrix = new Matrix(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } }); + + Assert.ThrowsException(() => matrix.Qr()); + } + } +} diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/EigenDecomposition.cs b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/EigenDecomposition.cs new file mode 100644 index 0000000..509708d --- /dev/null +++ b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/EigenDecomposition.cs @@ -0,0 +1,939 @@ +using System; +using CSharpNumerics.Numerics.Objects; + +namespace CSharpNumerics.Numerics.LinearAlgebra.Decompositions; + +/// +/// Eigenvalue decomposition A·V = V·D. +/// For symmetric matrices: Householder tridiagonalization followed by the implicit QL algorithm — +/// all eigenvalues are real, returned in ascending order with orthonormal eigenvectors. +/// For non-symmetric matrices: Hessenberg reduction followed by the shifted QR algorithm — +/// eigenvalues may be complex conjugate pairs (real parts in , +/// imaginary parts in ) and D is block diagonal. +/// Based on the EISPACK/JAMA reference implementation. +/// +public sealed class EigenDecomposition +{ + private const double SymmetryTolerance = 1e-10; + + private readonly int n; + private readonly double[] d; + private readonly double[] e; + private readonly double[,] v; + private double[,] hess; + private double[] ort; + private double cdivr, cdivi; + + public EigenDecomposition(Matrix matrix) + { + if (matrix.rowLength != matrix.columnLength) + { + throw new Exception("Is not a NxN matrix"); + } + + n = matrix.rowLength; + d = new double[n]; + e = new double[n]; + v = new double[n, n]; + IsSymmetric = CheckSymmetric(matrix); + + if (IsSymmetric) + { + for (var i = 0; i < n; i++) + { + for (var j = 0; j < n; j++) + { + v[i, j] = matrix.values[i, j]; + } + } + + Tridiagonalize(); + TridiagonalQl(); + } + else + { + hess = (double[,])matrix.values.Clone(); + ort = new double[n]; + + ReduceToHessenberg(); + HessenbergQr(); + } + } + + /// + /// True if the matrix was detected as symmetric (all eigenvalues real, eigenvectors orthonormal). + /// + public bool IsSymmetric { get; } + + /// + /// Real parts of the eigenvalues. For symmetric matrices these are all eigenvalues, in ascending order. + /// + public double[] RealEigenvalues => (double[])d.Clone(); + + /// + /// Imaginary parts of the eigenvalues (all zero for symmetric matrices). + /// + public double[] ImaginaryEigenvalues => (double[])e.Clone(); + + /// + /// The eigenvector matrix V (eigenvector k in column k, matching eigenvalue k). + /// + public Matrix EigenVectors => new Matrix((double[,])v.Clone()); + + /// + /// The block diagonal eigenvalue matrix D such that A·V = V·D. + /// Complex conjugate pairs appear as 2×2 blocks [[re, im], [-im, re]]. + /// + public Matrix DiagonalMatrix + { + get + { + var diagonal = new double[n, n]; + for (var i = 0; i < n; i++) + { + diagonal[i, i] = d[i]; + if (e[i] > 0) + { + diagonal[i, i + 1] = e[i]; + } + else if (e[i] < 0) + { + diagonal[i, i - 1] = e[i]; + } + } + return new Matrix(diagonal); + } + } + + private bool CheckSymmetric(Matrix matrix) + { + for (var i = 0; i < n; i++) + { + for (var j = i + 1; j < n; j++) + { + var a = matrix.values[i, j]; + var b = matrix.values[j, i]; + var scale = Math.Max(1.0, Math.Max(Math.Abs(a), Math.Abs(b))); + + if (Math.Abs(a - b) > SymmetryTolerance * scale) + { + return false; + } + } + } + return true; + } + + // Householder reduction of a symmetric matrix to tridiagonal form (EISPACK tred2). + private void Tridiagonalize() + { + for (var j = 0; j < n; j++) + { + d[j] = v[n - 1, j]; + } + + for (var i = n - 1; i > 0; i--) + { + var scale = 0.0; + var h = 0.0; + for (var k = 0; k < i; k++) + { + scale += Math.Abs(d[k]); + } + + if (scale == 0.0) + { + e[i] = d[i - 1]; + for (var j = 0; j < i; j++) + { + d[j] = v[i - 1, j]; + v[i, j] = 0.0; + v[j, i] = 0.0; + } + } + else + { + for (var k = 0; k < i; k++) + { + d[k] /= scale; + h += d[k] * d[k]; + } + + var f = d[i - 1]; + var g = Math.Sqrt(h); + if (f > 0) + { + g = -g; + } + + e[i] = scale * g; + h -= f * g; + d[i - 1] = f - g; + for (var j = 0; j < i; j++) + { + e[j] = 0.0; + } + + for (var j = 0; j < i; j++) + { + f = d[j]; + v[j, i] = f; + g = e[j] + v[j, j] * f; + for (var k = j + 1; k <= i - 1; k++) + { + g += v[k, j] * d[k]; + e[k] += v[k, j] * f; + } + e[j] = g; + } + + f = 0.0; + for (var j = 0; j < i; j++) + { + e[j] /= h; + f += e[j] * d[j]; + } + + var hh = f / (h + h); + for (var j = 0; j < i; j++) + { + e[j] -= hh * d[j]; + } + + for (var j = 0; j < i; j++) + { + f = d[j]; + g = e[j]; + for (var k = j; k <= i - 1; k++) + { + v[k, j] -= f * e[k] + g * d[k]; + } + d[j] = v[i - 1, j]; + v[i, j] = 0.0; + } + } + + d[i] = h; + } + + for (var i = 0; i < n - 1; i++) + { + v[n - 1, i] = v[i, i]; + v[i, i] = 1.0; + var h = d[i + 1]; + + if (h != 0.0) + { + for (var k = 0; k <= i; k++) + { + d[k] = v[k, i + 1] / h; + } + for (var j = 0; j <= i; j++) + { + var g = 0.0; + for (var k = 0; k <= i; k++) + { + g += v[k, i + 1] * v[k, j]; + } + for (var k = 0; k <= i; k++) + { + v[k, j] -= g * d[k]; + } + } + } + + for (var k = 0; k <= i; k++) + { + v[k, i + 1] = 0.0; + } + } + + for (var j = 0; j < n; j++) + { + d[j] = v[n - 1, j]; + v[n - 1, j] = 0.0; + } + v[n - 1, n - 1] = 1.0; + e[0] = 0.0; + } + + // Implicit QL algorithm with shifts for a symmetric tridiagonal matrix (EISPACK tql2). + // Eigenvalues end up sorted ascending with matching eigenvector columns. + private void TridiagonalQl() + { + for (var i = 1; i < n; i++) + { + e[i - 1] = e[i]; + } + e[n - 1] = 0.0; + + var f = 0.0; + var tst1 = 0.0; + var eps = Math.Pow(2.0, -52.0); + + for (var l = 0; l < n; l++) + { + tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l])); + var m = l; + while (m < n) + { + if (Math.Abs(e[m]) <= eps * tst1) + { + break; + } + m++; + } + + if (m > l) + { + do + { + var g = d[l]; + var p = (d[l + 1] - g) / (2.0 * e[l]); + var r = Hypot(p, 1.0); + if (p < 0) + { + r = -r; + } + + d[l] = e[l] / (p + r); + d[l + 1] = e[l] * (p + r); + var dl1 = d[l + 1]; + var h = g - d[l]; + for (var i = l + 2; i < n; i++) + { + d[i] -= h; + } + f += h; + + p = d[m]; + var c = 1.0; + var c2 = c; + var c3 = c; + var el1 = e[l + 1]; + var s = 0.0; + var s2 = 0.0; + + for (var i = m - 1; i >= l; i--) + { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = Hypot(p, e[i]); + e[i + 1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i + 1] = h + s * (c * g + s * d[i]); + + for (var k = 0; k < n; k++) + { + h = v[k, i + 1]; + v[k, i + 1] = s * v[k, i] + c * h; + v[k, i] = c * v[k, i] - s * h; + } + } + + p = -s * s2 * c3 * el1 * e[l] / dl1; + e[l] = s * p; + d[l] = c * p; + } while (Math.Abs(e[l]) > eps * tst1); + } + + d[l] += f; + e[l] = 0.0; + } + + for (var i = 0; i < n - 1; i++) + { + var k = i; + var p = d[i]; + for (var j = i + 1; j < n; j++) + { + if (d[j] < p) + { + k = j; + p = d[j]; + } + } + + if (k != i) + { + d[k] = d[i]; + d[i] = p; + for (var j = 0; j < n; j++) + { + (v[j, i], v[j, k]) = (v[j, k], v[j, i]); + } + } + } + } + + // Householder reduction of a general matrix to Hessenberg form (EISPACK orthes). + private void ReduceToHessenberg() + { + var low = 0; + var high = n - 1; + + for (var m = low + 1; m <= high - 1; m++) + { + var scale = 0.0; + for (var i = m; i <= high; i++) + { + scale += Math.Abs(hess[i, m - 1]); + } + + if (scale != 0.0) + { + var h = 0.0; + for (var i = high; i >= m; i--) + { + ort[i] = hess[i, m - 1] / scale; + h += ort[i] * ort[i]; + } + + var g = Math.Sqrt(h); + if (ort[m] > 0) + { + g = -g; + } + h -= ort[m] * g; + ort[m] -= g; + + for (var j = m; j < n; j++) + { + var f = 0.0; + for (var i = high; i >= m; i--) + { + f += ort[i] * hess[i, j]; + } + f /= h; + for (var i = m; i <= high; i++) + { + hess[i, j] -= f * ort[i]; + } + } + + for (var i = 0; i <= high; i++) + { + var f = 0.0; + for (var j = high; j >= m; j--) + { + f += ort[j] * hess[i, j]; + } + f /= h; + for (var j = m; j <= high; j++) + { + hess[i, j] -= f * ort[j]; + } + } + + ort[m] *= scale; + hess[m, m - 1] = scale * g; + } + } + + for (var i = 0; i < n; i++) + { + for (var j = 0; j < n; j++) + { + v[i, j] = i == j ? 1.0 : 0.0; + } + } + + for (var m = high - 1; m >= low + 1; m--) + { + if (hess[m, m - 1] != 0.0) + { + for (var i = m + 1; i <= high; i++) + { + ort[i] = hess[i, m - 1]; + } + + for (var j = m; j <= high; j++) + { + var g = 0.0; + for (var i = m; i <= high; i++) + { + g += ort[i] * v[i, j]; + } + g = (g / ort[m]) / hess[m, m - 1]; + for (var i = m; i <= high; i++) + { + v[i, j] += g * ort[i]; + } + } + } + } + } + + private void ComplexDivide(double xr, double xi, double yr, double yi) + { + if (Math.Abs(yr) > Math.Abs(yi)) + { + var r = yi / yr; + var dd = yr + r * yi; + cdivr = (xr + r * xi) / dd; + cdivi = (xi - r * xr) / dd; + } + else + { + var r = yr / yi; + var dd = yi + r * yr; + cdivr = (r * xr + xi) / dd; + cdivi = (r * xi - xr) / dd; + } + } + + // Shifted QR iteration on a Hessenberg matrix with eigenvector back-substitution (EISPACK hqr2). + private void HessenbergQr() + { + var nn = n; + var current = nn - 1; + var low = 0; + var high = nn - 1; + var eps = Math.Pow(2.0, -52.0); + var exshift = 0.0; + double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y; + + var norm = 0.0; + for (var i = 0; i < nn; i++) + { + if (i < low || i > high) + { + d[i] = hess[i, i]; + e[i] = 0.0; + } + for (var j = Math.Max(i - 1, 0); j < nn; j++) + { + norm += Math.Abs(hess[i, j]); + } + } + + var iter = 0; + while (current >= low) + { + var l = current; + while (l > low) + { + s = Math.Abs(hess[l - 1, l - 1]) + Math.Abs(hess[l, l]); + if (s == 0.0) + { + s = norm; + } + if (Math.Abs(hess[l, l - 1]) < eps * s) + { + break; + } + l--; + } + + if (l == current) + { + hess[current, current] += exshift; + d[current] = hess[current, current]; + e[current] = 0.0; + current--; + iter = 0; + } + else if (l == current - 1) + { + w = hess[current, current - 1] * hess[current - 1, current]; + p = (hess[current - 1, current - 1] - hess[current, current]) / 2.0; + q = p * p + w; + z = Math.Sqrt(Math.Abs(q)); + hess[current, current] += exshift; + hess[current - 1, current - 1] += exshift; + x = hess[current, current]; + + if (q >= 0) + { + z = p >= 0 ? p + z : p - z; + d[current - 1] = x + z; + d[current] = d[current - 1]; + if (z != 0.0) + { + d[current] = x - w / z; + } + e[current - 1] = 0.0; + e[current] = 0.0; + x = hess[current, current - 1]; + s = Math.Abs(x) + Math.Abs(z); + p = x / s; + q = z / s; + r = Math.Sqrt(p * p + q * q); + p /= r; + q /= r; + + for (var j = current - 1; j < nn; j++) + { + z = hess[current - 1, j]; + hess[current - 1, j] = q * z + p * hess[current, j]; + hess[current, j] = q * hess[current, j] - p * z; + } + + for (var i = 0; i <= current; i++) + { + z = hess[i, current - 1]; + hess[i, current - 1] = q * z + p * hess[i, current]; + hess[i, current] = q * hess[i, current] - p * z; + } + + for (var i = low; i <= high; i++) + { + z = v[i, current - 1]; + v[i, current - 1] = q * z + p * v[i, current]; + v[i, current] = q * v[i, current] - p * z; + } + } + else + { + d[current - 1] = x + p; + d[current] = x + p; + e[current - 1] = z; + e[current] = -z; + } + + current -= 2; + iter = 0; + } + else + { + x = hess[current, current]; + y = 0.0; + w = 0.0; + if (l < current) + { + y = hess[current - 1, current - 1]; + w = hess[current, current - 1] * hess[current - 1, current]; + } + + if (iter == 10 || iter == 20) + { + exshift += x; + for (var i = low; i <= current; i++) + { + hess[i, i] -= x; + } + s = Math.Abs(hess[current, current - 1]) + Math.Abs(hess[current - 1, current - 2]); + x = y = 0.75 * s; + w = -0.4375 * s * s; + } + + if (iter == 30) + { + s = (y - x) / 2.0; + s = s * s + w; + if (s > 0) + { + s = Math.Sqrt(s); + if (y < x) + { + s = -s; + } + s = x - w / ((y - x) / 2.0 + s); + for (var i = low; i <= current; i++) + { + hess[i, i] -= s; + } + exshift += s; + x = y = w = 0.964; + } + } + + iter++; + if (iter > 250) + { + throw new Exception("Eigenvalue iteration did not converge"); + } + + var m = current - 2; + while (m >= l) + { + z = hess[m, m]; + r = x - z; + s = y - z; + p = (r * s - w) / hess[m + 1, m] + hess[m, m + 1]; + q = hess[m + 1, m + 1] - z - r - s; + r = hess[m + 2, m + 1]; + s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + p /= s; + q /= s; + r /= s; + if (m == l) + { + break; + } + if (Math.Abs(hess[m, m - 1]) * (Math.Abs(q) + Math.Abs(r)) < + eps * (Math.Abs(p) * (Math.Abs(hess[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(hess[m + 1, m + 1])))) + { + break; + } + m--; + } + + for (var i = m + 2; i <= current; i++) + { + hess[i, i - 2] = 0.0; + if (i > m + 2) + { + hess[i, i - 3] = 0.0; + } + } + + for (var k = m; k <= current - 1; k++) + { + var notLast = k != current - 1; + if (k != m) + { + p = hess[k, k - 1]; + q = hess[k + 1, k - 1]; + r = notLast ? hess[k + 2, k - 1] : 0.0; + x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r); + if (x == 0.0) + { + continue; + } + p /= x; + q /= x; + r /= x; + } + + s = Math.Sqrt(p * p + q * q + r * r); + if (p < 0) + { + s = -s; + } + + if (s != 0) + { + if (k != m) + { + hess[k, k - 1] = -s * x; + } + else if (l != m) + { + hess[k, k - 1] = -hess[k, k - 1]; + } + + p += s; + x = p / s; + y = q / s; + z = r / s; + q /= p; + r /= p; + + for (var j = k; j < nn; j++) + { + p = hess[k, j] + q * hess[k + 1, j]; + if (notLast) + { + p += r * hess[k + 2, j]; + hess[k + 2, j] -= p * z; + } + hess[k, j] -= p * x; + hess[k + 1, j] -= p * y; + } + + for (var i = 0; i <= Math.Min(current, k + 3); i++) + { + p = x * hess[i, k] + y * hess[i, k + 1]; + if (notLast) + { + p += z * hess[i, k + 2]; + hess[i, k + 2] -= p * r; + } + hess[i, k] -= p; + hess[i, k + 1] -= p * q; + } + + for (var i = low; i <= high; i++) + { + p = x * v[i, k] + y * v[i, k + 1]; + if (notLast) + { + p += z * v[i, k + 2]; + v[i, k + 2] -= p * r; + } + v[i, k] -= p; + v[i, k + 1] -= p * q; + } + } + } + } + } + + if (norm == 0.0) + { + return; + } + + for (current = nn - 1; current >= 0; current--) + { + p = d[current]; + q = e[current]; + + if (q == 0) + { + var l = current; + hess[current, current] = 1.0; + for (var i = current - 1; i >= 0; i--) + { + w = hess[i, i] - p; + r = 0.0; + for (var j = l; j <= current; j++) + { + r += hess[i, j] * hess[j, current]; + } + + if (e[i] < 0.0) + { + z = w; + s = r; + } + else + { + l = i; + if (e[i] == 0.0) + { + hess[i, current] = w != 0.0 ? -r / w : -r / (eps * norm); + } + else + { + x = hess[i, i + 1]; + y = hess[i + 1, i]; + q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; + t = (x * s - z * r) / q; + hess[i, current] = t; + hess[i + 1, current] = Math.Abs(x) > Math.Abs(z) + ? (-r - w * t) / x + : (-s - y * t) / z; + } + + t = Math.Abs(hess[i, current]); + if (eps * t * t > 1) + { + for (var j = i; j <= current; j++) + { + hess[j, current] /= t; + } + } + } + } + } + else if (q < 0) + { + var l = current - 1; + + if (Math.Abs(hess[current, current - 1]) > Math.Abs(hess[current - 1, current])) + { + hess[current - 1, current - 1] = q / hess[current, current - 1]; + hess[current - 1, current] = -(hess[current, current] - p) / hess[current, current - 1]; + } + else + { + ComplexDivide(0.0, -hess[current - 1, current], hess[current - 1, current - 1] - p, q); + hess[current - 1, current - 1] = cdivr; + hess[current - 1, current] = cdivi; + } + + hess[current, current - 1] = 0.0; + hess[current, current] = 1.0; + for (var i = current - 2; i >= 0; i--) + { + var ra = 0.0; + var sa = 0.0; + for (var j = l; j <= current; j++) + { + ra += hess[i, j] * hess[j, current - 1]; + sa += hess[i, j] * hess[j, current]; + } + w = hess[i, i] - p; + + if (e[i] < 0.0) + { + z = w; + r = ra; + s = sa; + } + else + { + l = i; + if (e[i] == 0) + { + ComplexDivide(-ra, -sa, w, q); + hess[i, current - 1] = cdivr; + hess[i, current] = cdivi; + } + else + { + x = hess[i, i + 1]; + y = hess[i + 1, i]; + var vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; + var vi = (d[i] - p) * 2.0 * q; + if (vr == 0.0 && vi == 0.0) + { + vr = eps * norm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z)); + } + ComplexDivide(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi); + hess[i, current - 1] = cdivr; + hess[i, current] = cdivi; + + if (Math.Abs(x) > Math.Abs(z) + Math.Abs(q)) + { + hess[i + 1, current - 1] = (-ra - w * hess[i, current - 1] + q * hess[i, current]) / x; + hess[i + 1, current] = (-sa - w * hess[i, current] - q * hess[i, current - 1]) / x; + } + else + { + ComplexDivide(-r - y * hess[i, current - 1], -s - y * hess[i, current], z, q); + hess[i + 1, current - 1] = cdivr; + hess[i + 1, current] = cdivi; + } + } + + t = Math.Max(Math.Abs(hess[i, current - 1]), Math.Abs(hess[i, current])); + if (eps * t * t > 1) + { + for (var j = i; j <= current; j++) + { + hess[j, current - 1] /= t; + hess[j, current] /= t; + } + } + } + } + } + } + + for (var j = nn - 1; j >= low; j--) + { + for (var i = low; i <= high; i++) + { + z = 0.0; + for (var k = low; k <= Math.Min(j, high); k++) + { + z += v[i, k] * hess[k, j]; + } + v[i, j] = z; + } + } + } + + private static double Hypot(double a, double b) + { + if (Math.Abs(a) > Math.Abs(b)) + { + var r = b / a; + return Math.Abs(a) * Math.Sqrt(1 + r * r); + } + if (b != 0) + { + var r = a / b; + return Math.Abs(b) * Math.Sqrt(1 + r * r); + } + return 0.0; + } +} diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/QrDecomposition.cs b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/QrDecomposition.cs new file mode 100644 index 0000000..3a913ef --- /dev/null +++ b/Numerics/Numerics/Numerics/LinearAlgebra/Decompositions/QrDecomposition.cs @@ -0,0 +1,283 @@ +using System; +using System.Collections.Generic; +using CSharpNumerics.Numerics.Objects; + +namespace CSharpNumerics.Numerics.LinearAlgebra.Decompositions; + +/// +/// QR decomposition via Householder reflections: A = Q·R where Q has orthonormal columns +/// and R is upper triangular. For an m×n matrix with m ≥ n, Solve computes the least squares +/// solution min ‖A x − b‖, which makes QR the numerically stable choice for overdetermined +/// systems (curve fitting, linear regression) — no normal equations needed. +/// The factorization is computed once and can be reused for multiple solves. +/// +public sealed class QrDecomposition +{ + private readonly double[,] qr; + private readonly double[] rDiag; + private readonly int m; + private readonly int n; + + public QrDecomposition(Matrix matrix) + { + if (matrix.rowLength < matrix.columnLength) + { + throw new Exception("The matrix must have at least as many rows as columns"); + } + + m = matrix.rowLength; + n = matrix.columnLength; + qr = (double[,])matrix.values.Clone(); + rDiag = new double[n]; + + for (var k = 0; k < n; k++) + { + var norm = 0.0; + for (var i = k; i < m; i++) + { + norm = Hypot(norm, qr[i, k]); + } + + if (norm != 0.0) + { + if (qr[k, k] < 0) + { + norm = -norm; + } + + for (var i = k; i < m; i++) + { + qr[i, k] /= norm; + } + qr[k, k] += 1.0; + + for (var j = k + 1; j < n; j++) + { + var s = 0.0; + for (var i = k; i < m; i++) + { + s += qr[i, k] * qr[i, j]; + } + s = -s / qr[k, k]; + for (var i = k; i < m; i++) + { + qr[i, j] += s * qr[i, k]; + } + } + } + + rDiag[k] = -norm; + } + } + + /// + /// True if the columns of A are linearly independent. Diagonal entries of R are compared + /// against a tolerance relative to the largest one, since Householder rounding leaves + /// rank-deficient columns at ~1e-15 rather than exactly zero. + /// + public bool IsFullRank + { + get + { + var maxDiag = 0.0; + for (var j = 0; j < n; j++) + { + maxDiag = Math.Max(maxDiag, Math.Abs(rDiag[j])); + } + + if (maxDiag == 0.0) + { + return false; + } + + var threshold = maxDiag * Math.Max(m, n) * Math.Pow(2.0, -52.0); + for (var j = 0; j < n; j++) + { + if (Math.Abs(rDiag[j]) <= threshold) + { + return false; + } + } + return true; + } + } + + /// + /// The economy-size orthogonal factor Q (m×n) with orthonormal columns. + /// + public Matrix Q + { + get + { + var q = new double[m, n]; + for (var k = n - 1; k >= 0; k--) + { + for (var i = 0; i < m; i++) + { + q[i, k] = 0.0; + } + q[k, k] = 1.0; + + for (var j = k; j < n; j++) + { + if (qr[k, k] != 0.0) + { + var s = 0.0; + for (var i = k; i < m; i++) + { + s += qr[i, k] * q[i, j]; + } + s = -s / qr[k, k]; + for (var i = k; i < m; i++) + { + q[i, j] += s * qr[i, k]; + } + } + } + } + return new Matrix(q); + } + } + + /// + /// The upper triangular factor R (n×n). + /// + public Matrix R + { + get + { + var r = new double[n, n]; + for (var i = 0; i < n; i++) + { + for (var j = i; j < n; j++) + { + r[i, j] = i == j ? rDiag[i] : qr[i, j]; + } + } + return new Matrix(r); + } + } + + /// + /// Solves A x = b in the least squares sense: returns x minimizing ‖A x − b‖. + /// For square full-rank systems this is the exact solution. + /// + public VectorN Solve(VectorN b) + { + if (b.Length != m) + { + throw new Exception("The vector length must match the matrix row length"); + } + + var x = new double[m]; + for (var i = 0; i < m; i++) + { + x[i] = b[i]; + } + + SolveInPlace(x); + + var result = new double[n]; + Array.Copy(x, result, n); + return new VectorN(result); + } + + /// + /// Solves A x = b in the least squares sense. + /// + public List Solve(List b) + { + if (b.Count != m) + { + throw new Exception("The vector length must match the matrix row length"); + } + + var x = b.ToArray(); + SolveInPlace(x); + + var result = new List(n); + for (var i = 0; i < n; i++) + { + result.Add(x[i]); + } + return result; + } + + /// + /// Solves A X = B in the least squares sense for multiple right-hand sides (one per column of B). + /// + public Matrix Solve(Matrix b) + { + if (b.rowLength != m) + { + throw new Exception("The row length of B must match the matrix row length"); + } + + var columns = b.columnLength; + var x = new double[n, columns]; + + for (var c = 0; c < columns; c++) + { + var column = new double[m]; + for (var i = 0; i < m; i++) + { + column[i] = b.values[i, c]; + } + + SolveInPlace(column); + + for (var i = 0; i < n; i++) + { + x[i, c] = column[i]; + } + } + + return new Matrix(x); + } + + private void SolveInPlace(double[] x) + { + if (!IsFullRank) + { + throw new Exception("The matrix is rank deficient"); + } + + for (var k = 0; k < n; k++) + { + var s = 0.0; + for (var i = k; i < m; i++) + { + s += qr[i, k] * x[i]; + } + s = -s / qr[k, k]; + for (var i = k; i < m; i++) + { + x[i] += s * qr[i, k]; + } + } + + for (var k = n - 1; k >= 0; k--) + { + x[k] /= rDiag[k]; + for (var i = 0; i < k; i++) + { + x[i] -= x[k] * qr[i, k]; + } + } + } + + private static double Hypot(double a, double b) + { + if (Math.Abs(a) > Math.Abs(b)) + { + var r = b / a; + return Math.Abs(a) * Math.Sqrt(1 + r * r); + } + if (b != 0) + { + var r = a / b; + return Math.Abs(b) * Math.Sqrt(1 + r * r); + } + return 0.0; + } +} diff --git a/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs index 5689667..554b032 100644 --- a/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs +++ b/Numerics/Numerics/Numerics/LinearAlgebra/MatrixDecompositionExtensions.cs @@ -16,4 +16,16 @@ public static class MatrixDecompositionExtensions /// The returned object caches the factorization for reuse across multiple solves. /// public static CholeskyDecomposition Cholesky(this Matrix matrix) => new CholeskyDecomposition(matrix); + + /// + /// Computes the QR decomposition via Householder reflections: A = Q·R. + /// For overdetermined systems (rows > columns) Solve gives the least squares solution. + /// + public static QrDecomposition Qr(this Matrix matrix) => new QrDecomposition(matrix); + + /// + /// Computes the eigenvalue decomposition A·V = V·D. + /// Symmetric matrices give real, ascending eigenvalues with orthonormal eigenvectors. + /// + public static EigenDecomposition Eigen(this Matrix matrix) => new EigenDecomposition(matrix); } diff --git a/Numerics/Numerics/Physics/Quantum/SchrodingerExtensions.cs b/Numerics/Numerics/Physics/Quantum/SchrodingerExtensions.cs index fa804a0..4eac00d 100644 --- a/Numerics/Numerics/Physics/Quantum/SchrodingerExtensions.cs +++ b/Numerics/Numerics/Physics/Quantum/SchrodingerExtensions.cs @@ -1,4 +1,5 @@ using CSharpNumerics.Physics.Constants; +using CSharpNumerics.Numerics.LinearAlgebra.Decompositions; using CSharpNumerics.Numerics.Objects; using System; @@ -62,7 +63,10 @@ public static StationaryStates SolveStationaryStates( if (i < points - 1) h[i, i + 1] = -kinetic; } - var (energies, vectors) = SymmetricEigenSolver.Solve(h); + // Symmetric eigensolver: ascending eigenvalues, eigenvector k in column k. + var eigen = new EigenDecomposition(new Matrix(h)); + var energies = eigen.RealEigenvalues; + var vectors = eigen.EigenVectors; int count = states < 0 ? points : Math.Min(states, points); var selectedEnergies = new double[count]; @@ -71,7 +75,12 @@ public static StationaryStates SolveStationaryStates( for (int k = 0; k < count; k++) { selectedEnergies[k] = energies[k]; - waveFunctions[k] = NormaliseReal(vectors[k], dx); + + var vector = new double[points]; + for (int i = 0; i < points; i++) + vector[i] = vectors.values[i, k]; + + waveFunctions[k] = NormaliseReal(vector, dx); } return new StationaryStates(selectedEnergies, waveFunctions, grid, dx); diff --git a/Numerics/Numerics/Physics/Quantum/SymmetricEigenSolver.cs b/Numerics/Numerics/Physics/Quantum/SymmetricEigenSolver.cs deleted file mode 100644 index 10daacd..0000000 --- a/Numerics/Numerics/Physics/Quantum/SymmetricEigenSolver.cs +++ /dev/null @@ -1,102 +0,0 @@ -using System; - -namespace CSharpNumerics.Physics.Quantum; - -/// -/// Cyclic Jacobi eigenvalue solver for real symmetric matrices. Unlike power iteration with -/// deflation, it returns all eigenvalues and orthonormal eigenvectors accurately, -/// including the smallest ones — which is exactly what the Schrödinger Hamiltonian needs (the -/// low-lying energy levels are the smallest eigenvalues). Eigenvalues are returned in ascending -/// order with matching eigenvectors. -/// -internal static class SymmetricEigenSolver -{ - /// - /// Diagonalises a symmetric matrix. Returns ascending eigenvalues and the corresponding - /// eigenvectors as vectors[k] (length n). - /// - public static (double[] values, double[][] vectors) Solve(double[,] matrix, int maxSweeps = 100, double tolerance = 1e-12) - { - int n = matrix.GetLength(0); - if (matrix.GetLength(1) != n) - throw new ArgumentException("Matrix must be square.", nameof(matrix)); - - var a = (double[,])matrix.Clone(); - var v = new double[n, n]; - for (int i = 0; i < n; i++) v[i, i] = 1.0; - - for (int sweep = 0; sweep < maxSweeps; sweep++) - { - double off = OffDiagonalNorm(a, n); - if (off < tolerance) - break; - - for (int p = 0; p < n - 1; p++) - { - for (int q = p + 1; q < n; q++) - { - if (Math.Abs(a[p, q]) < 1e-300) - continue; - - // Rotation angle that zeroes the (p,q) entry of JᵀAJ: - // tan(2φ) = 2·a_pq / (a_qq − a_pp). - double phi = 0.5 * Math.Atan2(2.0 * a[p, q], a[q, q] - a[p, p]); - double c = Math.Cos(phi); - double s = Math.Sin(phi); - - // A := A·J (rotate columns p, q) - for (int i = 0; i < n; i++) - { - double aip = a[i, p], aiq = a[i, q]; - a[i, p] = (c * aip) - (s * aiq); - a[i, q] = (s * aip) + (c * aiq); - } - // A := Jᵀ·A (rotate rows p, q) - for (int j = 0; j < n; j++) - { - double apj = a[p, j], aqj = a[q, j]; - a[p, j] = (c * apj) - (s * aqj); - a[q, j] = (s * apj) + (c * aqj); - } - // Accumulate eigenvectors V := V·J - for (int i = 0; i < n; i++) - { - double vip = v[i, p], viq = v[i, q]; - v[i, p] = (c * vip) - (s * viq); - v[i, q] = (s * vip) + (c * viq); - } - } - } - } - - // Extract eigenvalues (diagonal) and sort ascending. - var values = new double[n]; - for (int i = 0; i < n; i++) values[i] = a[i, i]; - - var order = new int[n]; - for (int i = 0; i < n; i++) order[i] = i; - Array.Sort(order, (x, y) => values[x].CompareTo(values[y])); - - var sortedValues = new double[n]; - var vectors = new double[n][]; - for (int k = 0; k < n; k++) - { - int col = order[k]; - sortedValues[k] = values[col]; - var vec = new double[n]; - for (int i = 0; i < n; i++) vec[i] = v[i, col]; - vectors[k] = vec; - } - - return (sortedValues, vectors); - } - - private static double OffDiagonalNorm(double[,] a, int n) - { - double sum = 0.0; - for (int i = 0; i < n; i++) - for (int j = i + 1; j < n; j++) - sum += a[i, j] * a[i, j]; - return Math.Sqrt(2.0 * sum); - } -} diff --git a/docs/LinearAlgebraRoadMap.md b/docs/LinearAlgebraRoadMap.md index bd7165c..4573dab 100644 --- a/docs/LinearAlgebraRoadMap.md +++ b/docs/LinearAlgebraRoadMap.md @@ -112,11 +112,11 @@ Placeras i `Numerics/RootFinding/`. Direkt användbart i: `KeplerOrbit` (Keplers - [x] Enhetstester: kända faktoriseringar, singulära matriser, round-trip `A ≈ P·L·U` ### Phase 2 — QR & egendekomposition -- [ ] Implementera `QrDecomposition` (Householder) + minsta kvadrat-`Solve` -- [ ] Implementera symmetrisk `EigenDecomposition` (QR med shift) -- [ ] Implementera osymmetrisk egenlösare (Hessenberg + QR-iteration) -- [ ] Migrera kvantmodulens symmetriska egenlösare till den nya -- [ ] Enhetstester: ortogonalitet `QᵀQ = I`, egenpar-residualer `‖Av − λv‖` +- [x] Implementera `QrDecomposition` (Householder) + minsta kvadrat-`Solve` +- [x] Implementera symmetrisk `EigenDecomposition` (QR med shift) +- [x] Implementera osymmetrisk egenlösare (Hessenberg + QR-iteration) +- [x] Migrera kvantmodulens symmetriska egenlösare till den nya +- [x] Enhetstester: ortogonalitet `QᵀQ = I`, egenpar-residualer `‖Av − λv‖` ### Phase 3 — SVD - [ ] Implementera `SvdDecomposition` (Golub–Kahan)