From 1047f93f6886aaa421c4db73274ce384de1f82e1 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Fri, 29 Apr 2011 08:48:38 +0000 Subject: [PATCH] Added solve methods using double[][] to linear algebra decomposition solvers. JIRA: MATH-564 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1097730 13f79535-47bb-0310-9956-ffa450edef68 --- .../linear/CholeskyDecompositionImpl.java | 43 ++++++++++++--- .../math/linear/DecompositionSolver.java | 12 ++++ .../math/linear/EigenDecompositionImpl.java | 55 +++++++++++++------ .../math/linear/LUDecompositionImpl.java | 18 ++++-- .../math/linear/QRDecompositionImpl.java | 5 ++ .../SingularValueDecompositionImpl.java | 16 ++++++ src/site/xdoc/changes.xml | 3 + .../math/linear/CholeskySolverTest.java | 3 + .../commons/math/linear/EigenSolverTest.java | 4 ++ .../commons/math/linear/LUSolverTest.java | 3 + .../commons/math/linear/QRSolverTest.java | 3 + .../math/linear/SingularValueSolverTest.java | 3 + 12 files changed, 138 insertions(+), 30 deletions(-) diff --git a/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java index 5e6c88f6c..736a97463 100644 --- a/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java @@ -273,15 +273,34 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition { return new ArrayRealVector(solve(b.getDataRef()), false); } - /** {@inheritDoc} */ - public RealMatrix solve(RealMatrix b) { + /** Solve the linear equation A × X = B for matrices A. + *

The A matrix is implicit, it is provided by the underlying + * decomposition algorithm.

+ * @param b right-hand side of the equation A × X = B + * @param reUseB if true, the b array will be reused and returned, + * instead of being copied + * @return a matrix X that minimizes the two norm of A × X - B + * @throws org.apache.commons.math.exception.DimensionMismatchException + * if the matrices dimensions do not match. + * @throws SingularMatrixException + * if the decomposed matrix is singular. + */ + private double[][] solve(double[][] b, boolean reUseB) { final int m = lTData.length; - if (b.getRowDimension() != m) { - throw new DimensionMismatchException(b.getRowDimension(), m); + if (b.length != m) { + throw new DimensionMismatchException(b.length, m); } - final int nColB = b.getColumnDimension(); - double[][] x = b.getData(); + final int nColB = b[0].length; + final double[][] x; + if (reUseB) { + x = b; + } else { + x = new double[b.length][nColB]; + for (int i = 0; i < b.length; ++i) { + System.arraycopy(b[i], 0, x[i], 0, nColB); + } + } // Solve LY = b for (int j = 0; j < m; j++) { @@ -316,10 +335,20 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition { } } - return new Array2DRowRealMatrix(x, false); + return x; } + /** {@inheritDoc} */ + public double[][] solve(double[][] b) { + return solve(b, false); + } + + /** {@inheritDoc} */ + public RealMatrix solve(RealMatrix b) { + return new Array2DRowRealMatrix(solve(b.getData(), true), false); + } + /** {@inheritDoc} */ public RealMatrix getInverse() { return solve(MatrixUtils.createRealIdentityMatrix(lTData.length)); diff --git a/src/main/java/org/apache/commons/math/linear/DecompositionSolver.java b/src/main/java/org/apache/commons/math/linear/DecompositionSolver.java index aab8fb91b..919399291 100644 --- a/src/main/java/org/apache/commons/math/linear/DecompositionSolver.java +++ b/src/main/java/org/apache/commons/math/linear/DecompositionSolver.java @@ -59,6 +59,18 @@ public interface DecompositionSolver { */ RealVector solve(final RealVector b); + /** Solve the linear equation A × X = B for matrices A. + *

The A matrix is implicit, it is provided by the underlying + * decomposition algorithm.

+ * @param b right-hand side of the equation A × X = B + * @return a matrix X that minimizes the two norm of A × X - B + * @throws org.apache.commons.math.exception.DimensionMismatchException + * if the matrices dimensions do not match. + * @throws SingularMatrixException + * if the decomposed matrix is singular. + */ + double[][] solve(final double[][] b); + /** Solve the linear equation A × X = B for matrices A. *

The A matrix is implicit, it is provided by the underlying * decomposition algorithm.

diff --git a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java index 3f89d68a9..4fa192d2c 100644 --- a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java @@ -338,37 +338,48 @@ public class EigenDecompositionImpl implements EigenDecomposition { return new ArrayRealVector(bp, false); } - /** - * Solve the linear equation A × X = B for symmetric matrices A. - *

- * This method only find exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0. - *

- * @param b Right-hand side of the equation A × X = B - * @return a Matrix X that minimizes the two norm of A × X - B - * @throws DimensionMismatchException if the matrices dimensions do not match. - * @throws SingularMatrixException if the decomposed matrix is singular. + /** Solve the linear equation A × X = B for matrices A. + *

The A matrix is implicit, it is provided by the underlying + * decomposition algorithm.

+ * @param b right-hand side of the equation A × X = B + * @param reUseB if true, the b array will be reused and returned, + * instead of being copied + * @return a matrix X that minimizes the two norm of A × X - B + * @throws org.apache.commons.math.exception.DimensionMismatchException + * if the matrices dimensions do not match. + * @throws SingularMatrixException + * if the decomposed matrix is singular. */ - public RealMatrix solve(final RealMatrix b) { + private double[][] solve(double[][] b, boolean reUseB) { if (!isNonSingular()) { throw new SingularMatrixException(); } final int m = realEigenvalues.length; - if (b.getRowDimension() != m) { - throw new DimensionMismatchException(b.getRowDimension(), m); + if (b.length != m) { + throw new DimensionMismatchException(b.length, m); } - final int nColB = b.getColumnDimension(); - final double[][] bp = new double[m][nColB]; + final int nColB = b[0].length; + final double[][] bp; + if (reUseB) { + bp = b; + } else { + bp = new double[m][nColB]; + } + final double[] tmpCol = new double[m]; for (int k = 0; k < nColB; ++k) { + for (int i = 0; i < m; ++i) { + tmpCol[i] = b[i][k]; + bp[i][k] = 0; + } for (int i = 0; i < m; ++i) { final ArrayRealVector v = eigenvectors[i]; final double[] vData = v.getDataRef(); double s = 0; for (int j = 0; j < m; ++j) { - s += v.getEntry(j) * b.getEntry(j, k); + s += v.getEntry(j) * tmpCol[j]; } s /= realEigenvalues[i]; for (int j = 0; j < m; ++j) { @@ -377,10 +388,20 @@ public class EigenDecompositionImpl implements EigenDecomposition { } } - return MatrixUtils.createRealMatrix(bp); + return bp; } + /** {@inheritDoc} */ + public double[][] solve(double[][] b) { + return solve(b, false); + } + + /** {@inheritDoc} */ + public RealMatrix solve(RealMatrix b) { + return new Array2DRowRealMatrix(solve(b.getData(), true), false); + } + /** * Check if the decomposed matrix is non-singular. * @return true if the decomposed matrix is non-singular diff --git a/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java index 01fa921df..d3113a007 100644 --- a/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java @@ -337,17 +337,17 @@ public class LUDecompositionImpl implements LUDecomposition { } /** {@inheritDoc} */ - public RealMatrix solve(RealMatrix b) { + public double[][] solve(double[][] b) { final int m = pivot.length; - if (b.getRowDimension() != m) { - throw new DimensionMismatchException(b.getRowDimension(), m); + if (b.length != m) { + throw new DimensionMismatchException(b.length, m); } if (singular) { throw new SingularMatrixException(); } - final int nColB = b.getColumnDimension(); + final int nColB = b[0].length; // Apply permutations to b final double[][] bp = new double[m][nColB]; @@ -355,7 +355,7 @@ public class LUDecompositionImpl implements LUDecomposition { final double[] bpRow = bp[row]; final int pRow = pivot[row]; for (int col = 0; col < nColB; col++) { - bpRow[col] = b.getEntry(pRow, col); + bpRow[col] = b[pRow][col]; } } @@ -387,7 +387,13 @@ public class LUDecompositionImpl implements LUDecomposition { } } - return new Array2DRowRealMatrix(bp, false); + return bp; + + } + + /** {@inheritDoc} */ + public RealMatrix solve(RealMatrix b) { + return new Array2DRowRealMatrix(solve(b.getData()), false); } /** {@inheritDoc} */ diff --git a/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java index ba56b3cf7..5d2740242 100644 --- a/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java @@ -337,6 +337,11 @@ public class QRDecompositionImpl implements QRDecomposition { return new ArrayRealVector(solve(b.getDataRef()), false); } + /** {@inheritDoc} */ + public double[][] solve(double[][] b) { + return solve(new BlockRealMatrix(b)).getData(); + } + /** {@inheritDoc} */ public RealMatrix solve(RealMatrix b) { final int n = qrt.length; diff --git a/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java index e0d50ce92..811eaa71a 100644 --- a/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java @@ -301,6 +301,22 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio return pseudoInverse.operate(b); } + /** + * Solve the linear equation A × X = B in least square sense. + *

+ * The m×n matrix A may not be square, the solution X is such that + * ||A × X - B|| is minimal. + *

+ * + * @param b Right-hand side of the equation A × X = B + * @return a matrix X that minimizes the two norm of A × X - B + * @throws org.apache.commons.math.exception.DimensionMismatchException + * if the matrices dimensions do not match. + */ + public double[][] solve(final double[][] b) { + return pseudoInverse.multiply(MatrixUtils.createRealMatrix(b)).getData(); + } + /** * Solve the linear equation A × X = B in least square sense. *

diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 6e17ad20a..6b12a9b15 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,9 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added solve methods using double[][] to linear algebra decomposition solvers. + Added an interpolator adapter for data with known period. diff --git a/src/test/java/org/apache/commons/math/linear/CholeskySolverTest.java b/src/test/java/org/apache/commons/math/linear/CholeskySolverTest.java index a4cca7938..1f1cdb25d 100644 --- a/src/test/java/org/apache/commons/math/linear/CholeskySolverTest.java +++ b/src/test/java/org/apache/commons/math/linear/CholeskySolverTest.java @@ -81,6 +81,9 @@ public class CholeskySolverTest { // using RealMatrix Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13); + // using double[][] + Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 1.0e-13); + // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { Assert.assertEquals(0, diff --git a/src/test/java/org/apache/commons/math/linear/EigenSolverTest.java b/src/test/java/org/apache/commons/math/linear/EigenSolverTest.java index 70a279603..a8829c59d 100644 --- a/src/test/java/org/apache/commons/math/linear/EigenSolverTest.java +++ b/src/test/java/org/apache/commons/math/linear/EigenSolverTest.java @@ -120,6 +120,10 @@ public class EigenSolverTest { RealMatrix solution=es.solve(b); Assert.assertEquals(0, solution.subtract(xRef).getNorm(), 2.5e-12); + // using double[][] + solution = MatrixUtils.createRealMatrix(es.solve(b.getData())); + Assert.assertEquals(0, solution.subtract(xRef).getNorm(), 2.5e-12); + // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { Assert.assertEquals(0, diff --git a/src/test/java/org/apache/commons/math/linear/LUSolverTest.java b/src/test/java/org/apache/commons/math/linear/LUSolverTest.java index 61beccadd..c8d5fc681 100644 --- a/src/test/java/org/apache/commons/math/linear/LUSolverTest.java +++ b/src/test/java/org/apache/commons/math/linear/LUSolverTest.java @@ -143,6 +143,9 @@ public class LUSolverTest { // using RealMatrix Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13); + // using double[][] + Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 1.0e-13); + // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { Assert.assertEquals(0, diff --git a/src/test/java/org/apache/commons/math/linear/QRSolverTest.java b/src/test/java/org/apache/commons/math/linear/QRSolverTest.java index 0d2d4a42d..70fc1f336 100644 --- a/src/test/java/org/apache/commons/math/linear/QRSolverTest.java +++ b/src/test/java/org/apache/commons/math/linear/QRSolverTest.java @@ -137,6 +137,9 @@ public class QRSolverTest { // using RealMatrix Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm()); + // using double[][] + Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm()); + // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { final double[] x = solver.solve(b.getColumn(i)); diff --git a/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java b/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java index 4198cc324..4d301be70 100644 --- a/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java +++ b/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java @@ -101,6 +101,9 @@ public class SingularValueSolverTest { // using RealMatrix Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance); + // using double[][] + Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), normTolerance); + // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { Assert.assertEquals(0,