feat(maths): add LU decomposition algorithm using Doolittle method

MD-Mushfiqur123 · MD-Mushfiqur123 · commit 3f8bb42c5b6e · 2026-05-19T07:29:20.000+06:00
diff --git a/src/main/java/com/thealgorithms/maths/LUDecomposition.java b/src/main/java/com/thealgorithms/maths/LUDecomposition.java
@@ -0,0 +1,147 @@
+package com.thealgorithms.maths;
+
+/**
+ * LU Decomposition (Doolittle Algorithm)
+ *
+ * Decomposes a square matrix A into a lower triangular matrix L and an upper
+ * triangular matrix U such that A = L * U. The diagonal of L contains all ones.
+ *
+ * This decomposition is useful for solving systems of linear equations, computing
+ * matrix determinants, and finding matrix inverses efficiently.
+ *
+ * The algorithm proceeds by computing elements of L and U iteratively:
+ * - For U: U[k][j] = A[k][j] - sum(L[k][s] * U[s][j]) for s &lt; k
+ * - For L: L[i][k] = (A[i][k] - sum(L[i][s] * U[s][k])) / U[k][k] for s &lt; k
+ *
+ * Time Complexity: O(n^3) for decomposition, O(n^2) for solving a system.
+ * Space Complexity: O(n^2) for storing L and U matrices.
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/LU_decomposition">LU Decomposition</a>
+ */
+public final class LUDecomposition {
+
+    private LUDecomposition() {
+    }
+
+    /**
+     * Performs LU decomposition on a square matrix using the Doolittle algorithm.
+     *
+     * @param matrix a square matrix represented as a 2D array
+     * @return a 2D array where the lower triangle (excluding diagonal) contains L
+     *         elements (with implicit 1s on the diagonal) and the upper triangle
+     *         (including diagonal) contains U elements
+     * @throws IllegalArgumentException if the matrix is not square
+     * @throws ArithmeticException      if a zero pivot is encountered
+     */
+    public static double[][] decompose(double[][] matrix) {
+        int n = matrix.length;
+        for (double[] row : matrix) {
+            if (row.length != n) {
+                throw new IllegalArgumentException("Matrix must be square.");
+            }
+        }
+
+        double[][] lu = new double[n][n];
+        for (int i = 0; i < n; i++) {
+            for (int j = 0; j < n; j++) {
+                lu[i][j] = matrix[i][j];
+            }
+        }
+
+        for (int k = 0; k < n; k++) {
+            // Compute U elements for row k
+            for (int j = k; j < n; j++) {
+                double sum = 0;
+                for (int s = 0; s < k; s++) {
+                    sum += lu[k][s] * lu[s][j];
+                }
+                lu[k][j] -= sum;
+            }
+
+            // Check for zero pivot
+            if (lu[k][k] == 0) {
+                throw new ArithmeticException("Zero pivot encountered. Matrix may be singular.");
+            }
+
+            // Compute L elements for column k
+            for (int i = k + 1; i < n; i++) {
+                double sum = 0;
+                for (int s = 0; s < k; s++) {
+                    sum += lu[i][s] * lu[s][k];
+                }
+                lu[i][k] = (lu[i][k] - sum) / lu[k][k];
+            }
+        }
+
+        return lu;
+    }
+
+    /**
+     * Extracts the lower triangular matrix L from the combined LU matrix.
+     * The diagonal of L is set to 1 (Doolittle convention).
+     *
+     * @param lu the combined LU matrix from {@link #decompose(double[][])}
+     * @return the lower triangular matrix L with 1s on the diagonal
+     */
+    public static double[][] getLowerMatrix(double[][] lu) {
+        int n = lu.length;
+        double[][] lower = new double[n][n];
+        for (int i = 0; i < n; i++) {
+            lower[i][i] = 1.0;
+            for (int j = 0; j < i; j++) {
+                lower[i][j] = lu[i][j];
+            }
+        }
+        return lower;
+    }
+
+    /**
+     * Extracts the upper triangular matrix U from the combined LU matrix.
+     *
+     * @param lu the combined LU matrix from {@link #decompose(double[][])}
+     * @return the upper triangular matrix U
+     */
+    public static double[][] getUpperMatrix(double[][] lu) {
+        int n = lu.length;
+        double[][] upper = new double[n][n];
+        for (int i = 0; i < n; i++) {
+            for (int j = i; j < n; j++) {
+                upper[i][j] = lu[i][j];
+            }
+        }
+        return upper;
+    }
+
+    /**
+     * Solves a system of linear equations Ax = b using LU decomposition.
+     *
+     * @param lu the combined LU matrix from {@link #decompose(double[][])}
+     * @param b  the right-hand side vector
+     * @return the solution vector x
+     */
+    public static double[] solve(double[][] lu, double[] b) {
+        int n = lu.length;
+        double[] y = new double[n];
+        double[] x = new double[n];
+
+        // Forward substitution: solve Ly = b
+        for (int i = 0; i < n; i++) {
+            double sum = 0;
+            for (int j = 0; j < i; j++) {
+                sum += lu[i][j] * y[j];
+            }
+            y[i] = b[i] - sum;
+        }
+
+        // Back substitution: solve Ux = y
+        for (int i = n - 1; i >= 0; i--) {
+            double sum = 0;
+            for (int j = i + 1; j < n; j++) {
+                sum += lu[i][j] * x[j];
+            }
+            x[i] = (y[i] - sum) / lu[i][i];
+        }
+
+        return x;
+    }
+}
diff --git a/src/test/java/com/thealgorithms/maths/LUDecompositionTest.java b/src/test/java/com/thealgorithms/maths/LUDecompositionTest.java
@@ -0,0 +1,116 @@
+package com.thealgorithms.maths;
+
+import static org.junit.jupiter.api.Assertions.assertArrayEquals;
+import static org.junit.jupiter.api.Assertions.assertThrows;
+
+import org.junit.jupiter.api.Test;
+
+public class LUDecompositionTest {
+
+    private static final double DELTA = 1e-9;
+
+    @Test
+    void testDecomposeSimpleMatrix() {
+        double[][] matrix = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}};
+        double[][] lu = LUDecomposition.decompose(matrix);
+
+        double[][] lower = LUDecomposition.getLowerMatrix(lu);
+        double[][] upper = LUDecomposition.getUpperMatrix(lu);
+
+        // Verify L has 1s on diagonal
+        assertArrayEquals(new double[]{1, 1, 1}, new double[]{lower[0][0], lower[1][1], lower[2][2]}, DELTA);
+
+        // Verify L * U = A
+        double[][] product = multiply(lower, upper);
+        assertArrayEquals(new double[]{2, 1, 1, 4, 3, 3, 8, 7, 9}, flatten(product), DELTA);
+    }
+
+    @Test
+    void testDecomposeTwoByTwo() {
+        double[][] matrix = {{1, 2}, {3, 4}};
+        double[][] lu = LUDecomposition.decompose(matrix);
+
+        double[][] lower = LUDecomposition.getLowerMatrix(lu);
+        double[][] upper = LUDecomposition.getUpperMatrix(lu);
+
+        double[][] product = multiply(lower, upper);
+        assertArrayEquals(new double[]{1, 2, 3, 4}, flatten(product), DELTA);
+    }
+
+    @Test
+    void testDecomposeIdentityMatrix() {
+        double[][] matrix = {{1, 0}, {0, 1}};
+        double[][] lu = LUDecomposition.decompose(matrix);
+
+        double[][] lower = LUDecomposition.getLowerMatrix(lu);
+        double[][] upper = LUDecomposition.getUpperMatrix(lu);
+
+        // For identity matrix, L = I and U = I
+        assertArrayEquals(new double[]{1, 0, 0, 1}, flatten(lower), DELTA);
+        assertArrayEquals(new double[]{1, 0, 0, 1}, flatten(upper), DELTA);
+    }
+
+    @Test
+    void testDecomposeNonSquareMatrixThrows() {
+        double[][] matrix = {{1, 2, 3}, {4, 5, 6}};
+        assertThrows(IllegalArgumentException.class, () -> LUDecomposition.decompose(matrix));
+    }
+
+    @Test
+    void testDecomposeSingularMatrixThrows() {
+        double[][] matrix = {{0, 1}, {1, 0}};
+        assertThrows(ArithmeticException.class, () -> LUDecomposition.decompose(matrix));
+    }
+
+    @Test
+    void testSolveLinearSystem() {
+        // 2x + y + z = 8
+        // 4x + 3y + 3z = 20
+        // 8x + 7y + 9z = 46
+        // Solution: x=1, y=3, z=3
+        double[][] matrix = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}};
+        double[] b = {8, 20, 46};
+        double[][] lu = LUDecomposition.decompose(matrix);
+        double[] solution = LUDecomposition.solve(lu, b);
+
+        assertArrayEquals(new double[]{1, 3, 3}, solution, DELTA);
+    }
+
+    @Test
+    void testSolveTwoByTwoSystem() {
+        // 2x + y = 5
+        // x + 3y = 7
+        // Solution: x=1.6, y=1.8
+        double[][] matrix = {{2, 1}, {1, 3}};
+        double[] b = {5, 7};
+        double[][] lu = LUDecomposition.decompose(matrix);
+        double[] solution = LUDecomposition.solve(lu, b);
+
+        assertArrayEquals(new double[]{1.6, 1.8}, solution, DELTA);
+    }
+
+    private static double[][] multiply(double[][] a, double[][] b) {
+        int n = a.length;
+        double[][] result = new double[n][n];
+        for (int i = 0; i < n; i++) {
+            for (int j = 0; j < n; j++) {
+                for (int k = 0; k < n; k++) {
+                    result[i][j] += a[i][k] * b[k][j];
+                }
+            }
+        }
+        return result;
+    }
+
+    private static double[] flatten(double[][] matrix) {
+        int n = matrix.length;
+        double[] result = new double[n * n];
+        int idx = 0;
+        for (double[] row : matrix) {
+            for (double val : row) {
+                result[idx++] = val;
+            }
+        }
+        return result;
+    }
+}
