Fixed MatrixDense solution to underdetermined systems. Added Gaussian…

… Elimination

src/main/java/com/wildbitsfoundry/etk4j/math/interpolation/QuadraticSpline.java

@@ -1,27 +1,22 @@
 package com.wildbitsfoundry.etk4j.math.interpolation;
 
-import com.wildbitsfoundry.etk4j.math.linearalgebra.LUDecomposition;
+import com.wildbitsfoundry.etk4j.math.linearalgebra.GaussianElimination;
 import com.wildbitsfoundry.etk4j.math.linearalgebra.LUDecompositionDense;
 import com.wildbitsfoundry.etk4j.math.linearalgebra.MatrixDense;
 import com.wildbitsfoundry.etk4j.util.DoubleArrays;
 
-import java.lang.reflect.Array;
 import java.util.Arrays;
 
+import static com.wildbitsfoundry.etk4j.math.linearalgebra.Matrices.solveLDLtTridiagonalSystem;
+import static com.wildbitsfoundry.etk4j.math.linearalgebra.Matrices.solveLDUTridiagonalSystem;
+
 public class QuadraticSpline extends Spline {
 
     private static final double P5 = 0.5, P33 = 1.0 / 3.0;
 
-    protected QuadraticSpline(double[] x, double[] y, double[] coefficients) {
+    protected QuadraticSpline(double[] x, double[] coefficients) {
         super(x, 3);
-        final int n = this.x.length;
-        // compute coefficients
-        this.coefficients = new double[(n - 1) * 3];
-        for (int i = 0, j = 0; i < n - 1; ++i, ++j) {
-            this.coefficients[j] = coefficients[i + n - 1];
-            this.coefficients[++j] = coefficients[i];
-            this.coefficients[++j] = y[i];
-        }
+        this.coefficients = coefficients;
     }
 
     public static QuadraticSpline newNaturalSpline(double[] x, double[] y) {
@@ -48,17 +43,18 @@ public static QuadraticSpline newNaturalSplineInPlace(double[] x, double[] y) {
             A.set(j + n - 1, j + n - 1, 2 * hx);
         }
 
-        A.set(2 * n - 4, n - 1, 1); // c_0 = 0
-        A.set(2 * n - 3, 2 * n - 3, 1); // c_n-1 = 0
+        A.set(2 * n - 4, n - 1, 1); // c[0] = 0
+        A.set(2 * n - 3, 2 * n - 3, 1); // c[n - 1] = 0
 
-        LUDecompositionDense LU = new LUDecompositionDense(A);
-        if(LU.isSingular()) {
-            // Apply a little regularization factor if the x vectors is evenly spaced to help with singularity
-            A.addEquals(MatrixDense.identity(A.getRowCount(), A.getColumnCount()).multiply(1e-10));
-            LU = new LUDecompositionDense(A);
+        double[] solution = GaussianElimination.solve(A, b);
+        // compute coefficients
+        double[] coefficients = new double[(n - 1) * 3];
+        for (int i = 0, j = 0; i < n - 1; ++i, ++j) {
+            coefficients[j] = solution[i + n - 1];
+            coefficients[++j] = solution[i];
+            coefficients[++j] = y[i];
         }
-        double[] coefficients = LU.solve(b).getArray();
-        return new QuadraticSpline(x, y, coefficients);
+        return new QuadraticSpline(x, coefficients);
     }
 
     public static QuadraticSpline newClampedSpline(double[] x, double[] y, double m0, double mn) {
@@ -92,14 +88,15 @@ public static QuadraticSpline newClampedSplineInPlace(double[] x, double[] y, do
         A.set(2 * n - 3, 2 * n - 3, 2 * (x[n - 1] - x[n - 2]));
         b[2 * n - 3] = mn; // derivative at the end
 
-        LUDecompositionDense LU = new LUDecompositionDense(A);
-        if(LU.isSingular()) {
-            // Apply a little regularization factor if the x vectors is evenly spaced to help with singularity
-            A.addEquals(MatrixDense.identity(A.getRowCount(), A.getColumnCount()).multiply(1e-10));
-            LU = new LUDecompositionDense(A);
+        double[] solution = GaussianElimination.solve(A, b);
+        // compute coefficients
+        double[] coefficients = new double[(n - 1) * 3];
+        for (int i = 0, j = 0; i < n - 1; ++i, ++j) {
+            coefficients[j] = solution[i + n - 1];
+            coefficients[++j] = solution[i];
+            coefficients[++j] = y[i];
         }
-        double[] coefficients = LU.solve(b).getArray();
-        return new QuadraticSpline(x, y, coefficients);
+        return new QuadraticSpline(x, coefficients);
     }
 
     @Override
@@ -141,4 +138,45 @@ public String toString() {
         return sb.toString().replace("+ -", "- ").replace("- -", "+ ")
                 .replace("=  + ", "= ").replace("=  - ", "= -");
     }
+
+    public static void main(String[] args) {
+        double[] x = {0.0, 10.0, 15.0, 20.0, 22.5, 30.0};
+        double[] y = {0.0, 227.04, 362.78, 517.35, 602.97, 901.67};
+//        double[] x = {0, 1, 2, 3};
+//        double[] y = {1, Math.exp(1), Math.exp(2), Math.exp(3)};
+//        double[] x = {0, 1, 2, 3, 4};
+//        double[] y = {1, 2, 0, 2, 1};
+
+//        QuadraticSpline qs2 = QuadraticSpline.newNaturalSpline(x, y);
+//
+//        System.out.println(qs2.evaluateAt(3));
+//
+//        CubicSpline cs = CubicSpline.newNaturalSpline(x, y);
+//
+//        System.out.println(cs.evaluateAt(3));
+
+        QuadraticSpline qs3 = QuadraticSpline.newClampedSpline(x, y, 1, Math.exp(3));
+
+        System.out.println(qs3.evaluateAt(3));
+
+        System.out.println(qs3);
+
+
+//        CubicSpline cs3 = CubicSpline.newClampedSpline(x, y, 1, Math.exp(3));
+//
+//        System.out.println(cs3.evaluateAt(3));
+
+//        QuadraticSpline qs3 = QuadraticSpline.newNaturalSpline(x, y);
+//
+//        System.out.println(qs3.evaluateAt(3));
+//
+//        System.out.println(qs3);
+//
+//        QuadraticSpline qs4 = QuadraticSpline.newMeowMix2(x, y);
+//
+//        System.out.println(qs4.evaluateAt(3));
+//
+//        System.out.println(qs4);
+
+    }
 }

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/GaussianElimination.java

@@ -0,0 +1,53 @@
+package com.wildbitsfoundry.etk4j.math.linearalgebra;
+
+public final class GaussianElimination {
+
+    private GaussianElimination() {}
+
+    public static double[] solve(MatrixDense A, double[] b) {
+        return solve(A, b, 1e-10);
+    }
+
+    public static double[] solve(MatrixDense A, double[] b, double EPSILON) {
+        int n = b.length;
+
+        for (int p = 0; p < n; p++) {
+            // Find pivot row and swap
+            int max = p;
+            for (int i = p + 1; i < n; i++) {
+                if (Math.abs(A.unsafeGet(i, p)) > Math.abs(A.unsafeGet(max, p))) {
+                    max = i;
+                }
+            }
+            double[] temp = A.getRow(p);
+            A.setRow(p, A.getRow(max));
+            A.setRow(max, temp);
+            double t = b[p]; b[p] = b[max]; b[max] = t;
+
+            // Check for singular or nearly singular matrix
+            if (Math.abs(A.unsafeGet(p, p)) <= EPSILON) {
+                throw new ArithmeticException("Matrix is singular or nearly singular");
+            }
+
+            // Pivot within A and b
+            for (int i = p + 1; i < n; i++) {
+                double alpha = A.unsafeGet(i, p) / A.unsafeGet(p, p);
+                b[i] -= alpha * b[p];
+                for (int j = p; j < n; j++) {
+                    A.unsafeSet(i, j, A.unsafeGet(i, j) - alpha * A.unsafeGet(p, j));
+                }
+            }
+        }
+
+        // Back substitution
+        double[] x = new double[n];
+        for (int i = n - 1; i >= 0; i--) {
+            double sum = 0.0;
+            for (int j = i + 1; j < n; j++) {
+                sum += A.unsafeGet(i, j) * x[j];
+            }
+            x[i] = (b[i] - sum) / A.unsafeGet(i, i);
+        }
+        return x;
+    }
+}

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDense.java

@@ -923,19 +923,20 @@ public ComplexMatrixDense multiply(ComplexMatrixDense matrix) {
      * @return The solution to {@code Ax = b}
      */
     public MatrixDense solve(MatrixDense b) {
-
         if (rows == cols) { // Matrix is Squared
             return new LUDecompositionDense(this).solve(b);
         } else if (rows > cols) { // Matrix is thin (Overdetermined system)
             return new QRDecompositionDense(this).solve(b);
         } else { // Matrix is short and wide (Under-determined system)
-            QRDecompositionDense qr = this.transpose().QR();
-            MatrixDense R1 = forwardSubstitutionSolve(qr.getRT(), b);
-            R1.appendRows(cols - R1.rows);
-            return qr.QmultiplyX(R1);
+            return this.pinv().multiply(b);
         }
     }
 
+    // TODO document and write test
+    public MatrixDense solve(double[] b) {
+        return solve(new MatrixDense(b, b.length));
+    }
+
     /**
      * Solve X*A = B, which is also A'*X' = B'
      *