feat: Add LU decomposition algorithm for matrices using Doolittle's algorithm

src/main/java/com/thealgorithms/matrix/LUDecomposition.java

@@ -1,88 +1,186 @@
 package com.thealgorithms.matrix;
 
 /**
- * LU Decomposition algorithm
- * --------------------------
- * Decomposes a square matrix a into a product of two matrices:
- * a = l * u
- * where:
- * - l is a lower triangular matrix with 1s on its diagonal
- * - u is an upper triangular matrix
+ * LU Decomposition algorithm for square matrices.
+ * Decomposes a matrix A into L (lower triangular) and U (upper triangular)
+ * such that A = L * U
  *
- * Reference:
- * https://en.wikipedia.org/wiki/lu_decomposition
+ * <p>Time Complexity: O(n^3)
+ * <p>Space Complexity: O(n^2)
+ *
+ * @author Raghu0703
+ * @see <a href="https://en.wikipedia.org/wiki/LU_decomposition">LU Decomposition</a>
  */
 public final class LUDecomposition {
-
     private LUDecomposition() {
     }
 
     /**
-     * A helper class to store both l and u matrices
+     * Performs LU decomposition on a square matrix using Doolittle's method.
+     *
+     * @param matrix The input square matrix
+     * @return A Result object containing L and U matrices
+     * @throws IllegalArgumentException if matrix is not square or singular
      */
-    public static class LU {
-        double[][] l;
-        double[][] u;
+    public static Result decompose(double[][] matrix) {
+        int n = matrix.length;
 
-        LU(double[][] l, double[][] u) {
-            this.l = l;
-            this.u = u;
+        // Validate input
+        if (n == 0) {
+            throw new IllegalArgumentException("Matrix cannot be empty");
+        }
+        for (double[] row : matrix) {
+            if (row.length != n) {
+                throw new IllegalArgumentException("Matrix must be square");
+            }
         }
-    }
 
-    /**
-     * Performs LU Decomposition on a square matrix a
-     *
-     * @param a input square matrix
-     * @return LU object containing l and u matrices
-     */
-    public static LU decompose(double[][] a) {
-        int n = a.length;
         double[][] l = new double[n][n];
         double[][] u = new double[n][n];
 
+        // Initialize L with identity matrix
         for (int i = 0; i < n; i++) {
-            // upper triangular matrix
-            for (int k = i; k < n; k++) {
-                double sum = 0;
-                for (int j = 0; j < i; j++) {
-                    sum += l[i][j] * u[j][k];
+            l[i][i] = 1.0;
+        }
+
+        // Perform LU decomposition using Doolittle's method
+        for (int j = 0; j < n; j++) {
+            // Calculate U matrix elements
+            for (int i = 0; i <= j; i++) {
+                double sum = 0.0;
+                for (int k = 0; k < i; k++) {
+                    sum += l[i][k] * u[k][j];
                 }
-                u[i][k] = a[i][k] - sum;
+                u[i][j] = matrix[i][j] - sum;
             }
 
-            // lower triangular matrix
-            for (int k = i; k < n; k++) {
-                if (i == k) {
-                    l[i][i] = 1; // diagonal as 1
-                } else {
-                    double sum = 0;
-                    for (int j = 0; j < i; j++) {
-                        sum += l[k][j] * u[j][i];
-                    }
-                    l[k][i] = (a[k][i] - sum) / u[i][i];
+            // Calculate L matrix elements
+            for (int i = j + 1; i < n; i++) {
+                double sum = 0.0;
+                for (int k = 0; k < j; k++) {
+                    sum += l[i][k] * u[k][j];
                 }
+
+                if (Math.abs(u[j][j]) < 1e-10) {
+                    throw new IllegalArgumentException(
+                        "Matrix is singular or nearly singular"
+                    );
+                }
+
+                l[i][j] = (matrix[i][j] - sum) / u[j][j];
             }
         }
 
-        return new LU(l, u);
+        return new Result(l, u);
+    }
+
+    /**
+     * Main method for demonstration.
+     *
+     * @param args command line arguments (not used)
+     */
+    public static void main(String[] args) {
+        // Example from the issue
+        double[][] matrix = {
+            {2, -1, -2},
+            {-4, 6, 3},
+            {-4, -2, 8}
+        };
+
+        System.out.println("LU Decomposition Example");
+        System.out.println("========================\n");
+
+        Result result = decompose(matrix);
+
+        System.out.println("Original Matrix A:");
+        printMatrix(matrix);
+
+        System.out.println("\nLower Triangular Matrix L:");
+        printMatrix(result.getL());
+
+        System.out.println("\nUpper Triangular Matrix U:");
+        printMatrix(result.getU());
+
+        // Verify that L * U = A
+        System.out.println("\nVerification: L * U = A");
+        double[][] product = multiplyMatrices(result.getL(), result.getU());
+        printMatrix(product);
     }
 
     /**
-     * Utility function to print a matrix
+     * Helper method to print a matrix.
      *
-     * @param m matrix to print
+     * @param matrix the matrix to print
      */
-    public static void printMatrix(double[][] m) {
-        for (double[] row : m) {
-            System.out.print("[");
-            for (int j = 0; j < row.length; j++) {
-                System.out.printf("%7.3f", row[j]);
-                if (j < row.length - 1) {
+    private static void printMatrix(double[][] matrix) {
+        for (double[] row : matrix) {
+            System.out.print("{ ");
+            for (int i = 0; i < row.length; i++) {
+                System.out.printf("%.3f", row[i]);
+                if (i < row.length - 1) {
                     System.out.print(", ");
                 }
             }
-            System.out.println("]");
+            System.out.println(" }");
+        }
+    }
+
+    /**
+     * Helper method to multiply two matrices.
+     *
+     * @param a first matrix
+     * @param b second matrix
+     * @return the product matrix
+     */
+    private static double[][] multiplyMatrices(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;
+    }
+
+    /**
+     * Result class to hold L and U matrices from decomposition.
+     */
+    public static class Result {
+        private final double[][] lMatrix;
+        private final double[][] uMatrix;
+
+        /**
+         * Constructor for Result.
+         *
+         * @param l Lower triangular matrix
+         * @param u Upper triangular matrix
+         */
+        public Result(double[][] l, double[][] u) {
+            this.lMatrix = l;
+            this.uMatrix = u;
+        }
+
+        /**
+         * Gets the lower triangular matrix.
+         *
+         * @return L matrix
+         */
+        public double[][] getL() {
+            return lMatrix;
+        }
+
+        /**
+         * Gets the upper triangular matrix.
+         *
+         * @return U matrix
+         */
+        public double[][] getU() {
+            return uMatrix;
         }
     }
 }

src/test/java/com/thealgorithms/matrix/LUDecompositionTest.java

@@ -1,40 +1,107 @@
 package com.thealgorithms.matrix;
 
 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 {
+class LUDecompositionTest {
+
+    private static final double EPSILON = 1e-6;
 
     @Test
-    public void testLUDecomposition() {
-        double[][] a = {{4, 3}, {6, 3}};
+    void testBasicLUDecomposition() {
+        double[][] matrix = {
+            {2, -1, -2},
+            {-4, 6, 3},
+            {-4, -2, 8}
+        };
 
-        // Perform LU decomposition
-        LUDecomposition.LU lu = LUDecomposition.decompose(a);
-        double[][] l = lu.l;
-        double[][] u = lu.u;
+        LUDecomposition.Result result = LUDecomposition.decompose(matrix);
 
-        // Reconstruct a from l and u
-        double[][] reconstructed = multiplyMatrices(l, u);
+        double[][] expectedL = {
+            {1.0, 0.0, 0.0},
+            {-2.0, 1.0, 0.0},
+            {-2.0, -1.0, 1.0}
+        };
 
-        // Assert that reconstructed matrix matches original a
-        for (int i = 0; i < a.length; i++) {
-            assertArrayEquals(a[i], reconstructed[i], 1e-9);
-        }
+        double[][] expectedU = {
+            {2.0, -1.0, -2.0},
+            {0.0, 4.0, -1.0},
+            {0.0, 0.0, 3.0}
+        };
+
+        assertMatrixEquals(expectedL, result.getL());
+        assertMatrixEquals(expectedU, result.getU());
+    }
+
+    @Test
+    void testIdentityMatrix() {
+        double[][] identity = {
+            {1, 0, 0},
+            {0, 1, 0},
+            {0, 0, 1}
+        };
+
+        LUDecomposition.Result result = LUDecomposition.decompose(identity);
+
+        assertMatrixEquals(identity, result.getL());
+        assertMatrixEquals(identity, result.getU());
+    }
+
+    @Test
+    void testTwoByTwoMatrix() {
+        double[][] matrix = {
+            {4, 3},
+            {6, 3}
+        };
+
+        LUDecomposition.Result result = LUDecomposition.decompose(matrix);
+
+        double[][] expectedL = {
+            {1.0, 0.0},
+            {1.5, 1.0}
+        };
+
+        double[][] expectedU = {
+            {4.0, 3.0},
+            {0.0, -1.5}
+        };
+
+        assertMatrixEquals(expectedL, result.getL());
+        assertMatrixEquals(expectedU, result.getU());
+    }
+
+    @Test
+    void testNonSquareMatrix() {
+        double[][] nonSquare = {
+            {1, 2, 3},
+            {4, 5, 6}
+        };
+
+        assertThrows(IllegalArgumentException.class, () -> LUDecomposition.decompose(nonSquare));
+    }
+
+    @Test
+    void testEmptyMatrix() {
+        double[][] empty = {};
+
+        assertThrows(IllegalArgumentException.class, () -> LUDecomposition.decompose(empty));
+    }
+
+    @Test
+    void testSingularMatrix() {
+        double[][] singular = {
+            {1, 2, 3},
+            {2, 4, 6},
+            {3, 6, 9}
+        };
+
+        assertThrows(IllegalArgumentException.class, () -> LUDecomposition.decompose(singular));
     }
 
-    // Helper method to multiply two matrices
-    private double[][] multiplyMatrices(double[][] a, double[][] b) {
-        int n = a.length;
-        double[][] c = 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++) {
-                    c[i][j] += a[i][k] * b[k][j];
-                }
-            }
+    private void assertMatrixEquals(double[][] expected, double[][] actual) {
+        for (int i = 0; i < expected.length; i++) {
+            assertArrayEquals(expected[i], actual[i], EPSILON);
         }
-        return c;
     }
 }