diff --git a/java/com.sap.sailing.polars.test/src/com/sap/sailing/polars/regression/test/IncrementalLeastSquaresTest.java b/java/com.sap.sailing.polars.test/src/com/sap/sailing/polars/regression/test/IncrementalLeastSquaresTest.java index 0e44953b422..55e87027925 100644 --- a/java/com.sap.sailing.polars.test/src/com/sap/sailing/polars/regression/test/IncrementalLeastSquaresTest.java +++ b/java/com.sap.sailing.polars.test/src/com/sap/sailing/polars/regression/test/IncrementalLeastSquaresTest.java @@ -87,7 +87,7 @@ public class IncrementalLeastSquaresTest { } @Test - public void testIncrementalRegressionWithRealDataNoIntercept() throws NotEnoughDataHasBeenAddedException { + public void testIncrementalRegressionWithRealDataNoInterceptWithSymbollicInversion() throws NotEnoughDataHasBeenAddedException { IncrementalLeastSquares leastSquares = new IncrementalAnyOrderLeastSquaresImpl(3, false); fillRegression(leastSquares); @@ -97,6 +97,65 @@ public class IncrementalLeastSquaresTest { assertThat(coeffs[1], closeTo(0.9317759507, ERROR)); assertThat(coeffs[2], closeTo(-0.0364029315, ERROR)); assertThat(coeffs[3], closeTo(0.0005326266, ERROR)); + + leastSquares.addData(0, 0); + long startTime = System.nanoTime(); + leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + System.out.println("With symbollic inversion no intercept" + (System.nanoTime() - startTime) + "ns"); + } + + @Test + public void testIncrementalRegressionWithRealDataNoInterceptWithoutSymbollicInversion() throws NotEnoughDataHasBeenAddedException { + + IncrementalLeastSquares leastSquares = new IncrementalAnyOrderLeastSquaresImpl(3, false, false); + fillRegression(leastSquares); + double[] coeffs = leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + + assertThat(coeffs[0], closeTo(0, ERROR)); + assertThat(coeffs[1], closeTo(0.9317759507, ERROR)); + assertThat(coeffs[2], closeTo(-0.0364029315, ERROR)); + assertThat(coeffs[3], closeTo(0.0005326266, ERROR)); + + leastSquares.addData(0, 0); + long startTime = System.nanoTime(); + leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + System.out.println("With numeric inversion no intercept " + (System.nanoTime() - startTime) + "ns"); + } + + @Test + public void testIncrementalRegressionWithRealDataWithInterceptWithSymbollicInversion() throws NotEnoughDataHasBeenAddedException { + + IncrementalLeastSquares leastSquares = new IncrementalAnyOrderLeastSquaresImpl(3, true, true); + fillRegression(leastSquares); + double[] coeffs = leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + + assertThat(coeffs[0], closeTo(6.319899225, ERROR)); + assertThat(coeffs[1], closeTo(-0.427890096, ERROR)); + assertThat(coeffs[2], closeTo(0.057763706, ERROR)); + assertThat(coeffs[3], closeTo(-0.001574209, ERROR)); + + leastSquares.addData(0, 0); + long startTime = System.nanoTime(); + leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + System.out.println("With symbollic inversion and intercept " + (System.nanoTime() - startTime) + "ns"); + } + + @Test + public void testIncrementalRegressionWithRealDataWithInterceptWithoutSymbollicInversion() throws NotEnoughDataHasBeenAddedException { + + IncrementalLeastSquares leastSquares = new IncrementalAnyOrderLeastSquaresImpl(3, true, false); + fillRegression(leastSquares); + double[] coeffs = leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + + assertThat(coeffs[0], closeTo(6.319899225, ERROR)); + assertThat(coeffs[1], closeTo(-0.427890096, ERROR)); + assertThat(coeffs[2], closeTo(0.057763706, ERROR)); + assertThat(coeffs[3], closeTo(-0.001574209, ERROR)); + + leastSquares.addData(0, 0); + long startTime = System.nanoTime(); + leastSquares.getOrCreatePolynomialFunction().getCoefficients(); + System.out.println("With numeric inversion and intercept " + (System.nanoTime() - startTime) + "ns"); } private void fillRegression(IncrementalLeastSquares leastSquares) { diff --git a/java/com.sap.sailing.polars/src/com/sap/sailing/polars/regression/impl/IncrementalAnyOrderLeastSquaresImpl.java b/java/com.sap.sailing.polars/src/com/sap/sailing/polars/regression/impl/IncrementalAnyOrderLeastSquaresImpl.java index 7a824bfc331..32f8dd8a7b4 100644 --- a/java/com.sap.sailing.polars/src/com/sap/sailing/polars/regression/impl/IncrementalAnyOrderLeastSquaresImpl.java +++ b/java/com.sap.sailing.polars/src/com/sap/sailing/polars/regression/impl/IncrementalAnyOrderLeastSquaresImpl.java @@ -25,7 +25,7 @@ import com.sap.sse.concurrent.NamedReentrantReadWriteLock; * */ public class IncrementalAnyOrderLeastSquaresImpl implements IncrementalLeastSquares { - + private double[][] matrixOfXSums; private double[] vectorOfXYMultSums; @@ -45,10 +45,12 @@ public class IncrementalAnyOrderLeastSquaresImpl implements IncrementalLeastSqua private PolynomialFunction cachedFunction; private final boolean hasIntercept; + private final boolean useSymbollicInversionIfPossible; - public IncrementalAnyOrderLeastSquaresImpl(int polynomialOrder, boolean hasIntercept) { + public IncrementalAnyOrderLeastSquaresImpl(int polynomialOrder, boolean hasIntercept, boolean useSymbollicInversionIfPossible) { this.hasIntercept = hasIntercept; this.polynomialOrder = polynomialOrder; + this.useSymbollicInversionIfPossible = useSymbollicInversionIfPossible; if (hasIntercept) { matrixOfXSums = new double[polynomialOrder + 1][polynomialOrder + 1]; vectorOfXYMultSums = new double[polynomialOrder + 1]; @@ -59,7 +61,11 @@ public class IncrementalAnyOrderLeastSquaresImpl implements IncrementalLeastSqua } public IncrementalAnyOrderLeastSquaresImpl(int polynomialOrder) { - this(polynomialOrder, true); + this(polynomialOrder, true, true); + } + + public IncrementalAnyOrderLeastSquaresImpl(int polynomialOrder, boolean hasIntercept) { + this(polynomialOrder, hasIntercept, true); } @Override @@ -127,7 +133,14 @@ public class IncrementalAnyOrderLeastSquaresImpl implements IncrementalLeastSqua } finally { LockUtil.unlockAfterRead(lock); } - RealMatrix inversedMatrix = new LUDecompositionImpl(matrixOfXSumsCopy).getSolver().getInverse(); + LeastSquaresInversionHelperType type = LeastSquaresInversionHelperType.getType(polynomialOrder, hasIntercept); + RealMatrix inversedMatrix = null; + if (useSymbollicInversionIfPossible && type != null) { + inversedMatrix = type.inverseMatrix(matrixOfXSumsCopy); + } + if (inversedMatrix == null) { + inversedMatrix = new LUDecompositionImpl(matrixOfXSumsCopy).getSolver().getInverse(); + } RealVector coeffs = inversedMatrix.operate(vectorOfXYMultSumsCopy); if (!hasIntercept) { resultFunction = createFunctionWithZeroIntercept(coeffs); @@ -159,5 +172,90 @@ public class IncrementalAnyOrderLeastSquaresImpl implements IncrementalLeastSqua } return new PolynomialFunction(coeff); } + + private enum LeastSquaresInversionHelperType { + + CUBIC, + CUBIC_NO_INTERCEPT; + + + /** + * + * @param realMatrix matrix to inverse + * @return Inversed matrix if symbolic inversion was found. Else null + */ + private RealMatrix inverseMatrix(RealMatrix realMatrix) { + RealMatrix result = null; + switch (this) { + case CUBIC_NO_INTERCEPT: + double[][] originDataCNI = realMatrix.getData(); + double a = originDataCNI[0][0]; + double b = originDataCNI[0][1]; + double c = originDataCNI[0][2]; + double d = originDataCNI[1][2]; + double e = originDataCNI[2][2]; + double[][] matrixDataCNI = { { d * d - c * e, b * e - c * d, c * c - b * d }, + { b * e - c * d, c * c - a * e, a * d - b * c }, + { c * c - b * d, a * d - b * c, b * b - a * c } }; + RealMatrix intermediateMatrixCNI = MatrixUtils.createRealMatrix(matrixDataCNI); + double factorCNI = (1.0 / (-a*c*e + a*d*d + b*b*e-2*b*c*d+c*c*c)); + result = intermediateMatrixCNI.scalarMultiply(factorCNI); + break; + case CUBIC: + double[][] originDataC = realMatrix.getData(); + double n = originDataC[0][0]; + a = originDataC[0][1]; + b = originDataC[0][2]; + c = originDataC[0][3]; + d = originDataC[1][3]; + e = originDataC[2][3]; + double f = originDataC[3][3]; + double[][] matrixDataC = { + { -d * d * d + 2 * c * e * d + b * f * d - b * e * e - c * c * f, + -e * c * c + d * d * c + b * f * c + a * e * e - b * d * e - a * d * f, + -f * b * b + d * d * b + c * e * b - c * c * d - a * d * e + a * c * f, + c * c * c - 2 * b * d * c - a * e * c + a * d * d + b * b * e }, + { -e * c * c + d * d * c + b * f * c + a * e * e - b * d * e - a * d * f, + -f * b * b + 2 * c * e * b - c * c * d - e * e * n + d * f * n, + c * c * c - b * d * c - a * e * c - f * n * c + a * b * f + d * e * n, + d * b * b - c * c * b - a * e * b + a * c * d - d * d * n + c * e * n }, + { -f * b * b + d * d * b + c * e * b - c * c * d - a * d * e + a * c * f, + c * c * c - b * d * c - a * e * c - f * n * c + a * b * f + d * e * n, + -f * a * a + 2 * c * d * a - b * c * c - d * d * n + b * f * n, + e * a * a - c * c * a - b * d * a + b * b * c + c * d * n - b * e * n }, + { c * c * c - 2 * b * d * c - a * e * c + a * d * d + b * b * e, + d * b * b - c * c * b - a * e * b + a * c * d - d * d * n + c * e * n, + e * a * a - c * c * a - b * d * a + b * b * c + c * d * n - b * e * n, + -b * b * b + 2 * a * c * b + d * n * b - a * a * d - c * c * n } }; + RealMatrix intermediateMatrixC = MatrixUtils.createRealMatrix(matrixDataC); + double factorC = 1.0 / (-a * a * d * f + a * a * e * e + 2 * c + * (a * (b * f + d * d) + e * (b * b + d * n)) - c * c * (2 * a * e + 3 * b * d + f * n) + b + * (-2 * a * d * e + d * f * n + e * e * -n) - b * b * b * f + b * b * d * d + c * c * c * c - d + * d * d * n); + result = intermediateMatrixC.scalarMultiply(factorC); + break; + default: + break; + } + return result; + } + + /** + * + * @param order Order of the polynomial + * @param hasIntercept If the least squares is set to no intercept + * @return type of the polynomial or null if it doesn't exist in this implementation + */ + private static LeastSquaresInversionHelperType getType(int order, boolean hasIntercept) { + LeastSquaresInversionHelperType result = null; + switch(order) { + case 3: + result = hasIntercept ? CUBIC : CUBIC_NO_INTERCEPT; + break; + } + return result; + } + + } }