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] 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`