From 008807938f612951c7757ab5f6af6edc9f8bca92 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sun, 21 Dec 2008 22:05:06 +0000 Subject: [PATCH] added getSolver() into LUDecomposition git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@728528 13f79535-47bb-0310-9956-ffa450edef68 --- .../math/linear/LUDecompositionImpl.java | 236 ++++++++++++++++-- .../apache/commons/math/linear/LUSolver.java | 173 ++----------- src/site/xdoc/userguide/linear.xml | 9 +- 3 files changed, 240 insertions(+), 178 deletions(-) diff --git a/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java b/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java index 631764a6f..dd3247d2a 100644 --- a/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java +++ b/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java @@ -17,7 +17,6 @@ package org.apache.commons.math.linear; - /** * Calculates the LUP-decomposition of a square matrix. *

The LUP-decomposition of a matrix A consists of three matrices @@ -33,7 +32,7 @@ package org.apache.commons.math.linear; public class LUDecompositionImpl implements LUDecomposition { /** Serializable version identifier. */ - private static final long serialVersionUID = 3446121671437672843L; + private static final long serialVersionUID = 1954692554563387537L; /** Entries of LU decomposition. */ private double lu[][]; @@ -61,20 +60,16 @@ public class LUDecompositionImpl implements LUDecomposition { /** * Calculates the LU-decomposition of the given matrix. - *

Calling this constructor is equivalent to first call the no-arguments - * constructor and then call {@link #decompose(RealMatrix)}.

* @param matrix The matrix to decompose. * @exception InvalidMatrixException if matrix is not square */ public LUDecompositionImpl(RealMatrix matrix) throws InvalidMatrixException { - decompose(matrix); + this(matrix, DEFAULT_TOO_SMALL); } /** * Calculates the LU-decomposition of the given matrix. - *

Calling this constructor is equivalent to first call the no-arguments - * constructor and then call {@link #decompose(RealMatrix, double)}.

* @param matrix The matrix to decompose. * @param singularityThreshold threshold (based on partial row norm) * under which a matrix is considered singular @@ -82,21 +77,11 @@ public class LUDecompositionImpl implements LUDecomposition { */ public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold) throws InvalidMatrixException { - decompose(matrix, singularityThreshold); - } - /** {@inheritDoc} */ - public void decompose(RealMatrix matrix) - throws InvalidMatrixException { - decompose(matrix, DEFAULT_TOO_SMALL); - } - - /** {@inheritDoc} */ - public void decompose(RealMatrix matrix, double singularityThreshold) - throws InvalidMatrixException { if (!matrix.isSquare()) { throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); } + final int m = matrix.getColumnDimension(); lu = matrix.getData(); pivot = new int[m]; @@ -222,7 +207,26 @@ public class LUDecompositionImpl implements LUDecomposition { /** {@inheritDoc} */ public int[] getPivot() throws IllegalStateException { - return pivot; + return pivot.clone(); + } + + /** {@inheritDoc} */ + public double getDeterminant() { + if (singular) { + return 0; + } else { + final int m = pivot.length; + double determinant = even ? 1 : -1; + for (int i = 0; i < m; i++) { + determinant *= lu[i][i]; + } + return determinant; + } + } + + /** {@inheritDoc} */ + public boolean isSingular() { + return singular; } /** {@inheritDoc} */ @@ -231,9 +235,197 @@ public class LUDecompositionImpl implements LUDecomposition { } /** {@inheritDoc} */ - public boolean isSingular() - throws IllegalStateException { - return singular; + public DecompositionSolver getSolver() { + return new Solver(lu, pivot, singular); + } + + private static class Solver implements DecompositionSolver { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -6353105415121373022L; + + /** Entries of LU decomposition. */ + private final double lu[][]; + + /** Pivot permutation associated with LU decomposition. */ + private final int[] pivot; + + /** Singularity indicator. */ + private final boolean singular; + + /** + * Build a solver from decomposed matrix. + * @param lu entries of LU decomposition + * @param pivot pivot permutation associated with LU decomposition + * @param singular singularity indicator + */ + private Solver(final double[][] lu, final int[] pivot, final boolean singular) { + this.lu = lu; + this.pivot = pivot; + this.singular = singular; + } + + /** {@inheritDoc} */ + public boolean isNonSingular() { + return !singular; + } + + /** {@inheritDoc} */ + public double[] solve(double[] b) + throws IllegalStateException, IllegalArgumentException, InvalidMatrixException { + + final int m = pivot.length; + if (b.length != m) { + throw new IllegalArgumentException("constant vector has wrong length"); + } + if (singular) { + throw new SingularMatrixException(); + } + + final double[] bp = new double[m]; + + // Apply permutations to b + for (int row = 0; row < m; row++) { + bp[row] = b[pivot[row]]; + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + for (int i = col + 1; i < m; i++) { + bp[i] -= bp[col] * lu[i][col]; + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + bp[col] /= lu[col][col]; + for (int i = 0; i < col; i++) { + bp[i] -= bp[col] * lu[i][col]; + } + } + + return bp; + + } + + /** {@inheritDoc} */ + public RealVector solve(RealVector b) + throws IllegalStateException, IllegalArgumentException, InvalidMatrixException { + try { + return solve((RealVectorImpl) b); + } catch (ClassCastException cce) { + + final int m = pivot.length; + if (b.getDimension() != m) { + throw new IllegalArgumentException("constant vector has wrong length"); + } + if (singular) { + throw new SingularMatrixException(); + } + + final double[] bp = new double[m]; + + // Apply permutations to b + for (int row = 0; row < m; row++) { + bp[row] = b.getEntry(pivot[row]); + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + for (int i = col + 1; i < m; i++) { + bp[i] -= bp[col] * lu[i][col]; + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + bp[col] /= lu[col][col]; + for (int i = 0; i < col; i++) { + bp[i] -= bp[col] * lu[i][col]; + } + } + + return new RealVectorImpl(bp, false); + + } + } + + /** Solve the linear equation A × X = B. + *

The A matrix is implicit here. It is

+ * @param b right-hand side of the equation A × X = B + * @return a vector X such that A × X = B + * @exception IllegalStateException if {@link #decompose(RealMatrix) decompose} + * has not been called + * @exception IllegalArgumentException if matrices dimensions don't match + * @exception InvalidMatrixException if decomposed matrix is singular + */ + public RealVectorImpl solve(RealVectorImpl b) + throws IllegalStateException, IllegalArgumentException, InvalidMatrixException { + return new RealVectorImpl(solve(b.getDataRef()), false); + } + + /** {@inheritDoc} */ + public RealMatrix solve(RealMatrix b) + throws IllegalStateException, IllegalArgumentException, InvalidMatrixException { + + final int m = pivot.length; + if (b.getRowDimension() != m) { + throw new IllegalArgumentException("Incorrect row dimension"); + } + if (singular) { + throw new SingularMatrixException(); + } + + final int nColB = b.getColumnDimension(); + + // Apply permutations to b + final double[][] bp = new double[m][nColB]; + for (int row = 0; row < m; row++) { + final double[] bpRow = bp[row]; + final int pRow = pivot[row]; + for (int col = 0; col < nColB; col++) { + bpRow[col] = b.getEntry(pRow, col); + } + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + final double[] bpCol = bp[col]; + for (int i = col + 1; i < m; i++) { + final double[] bpI = bp[i]; + final double luICol = lu[i][col]; + for (int j = 0; j < nColB; j++) { + bpI[j] -= bpCol[j] * luICol; + } + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + final double[] bpCol = bp[col]; + final double luDiag = lu[col][col]; + for (int j = 0; j < nColB; j++) { + bpCol[j] /= luDiag; + } + for (int i = 0; i < col; i++) { + final double[] bpI = bp[i]; + final double luICol = lu[i][col]; + for (int j = 0; j < nColB; j++) { + bpI[j] -= bpCol[j] * luICol; + } + } + } + + return new RealMatrixImpl(bp, false); + + } + + /** {@inheritDoc} */ + public RealMatrix getInverse() + throws IllegalStateException, InvalidMatrixException { + return solve(MatrixUtils.createRealIdentityMatrix(pivot.length)); + } + } } diff --git a/src/java/org/apache/commons/math/linear/LUSolver.java b/src/java/org/apache/commons/math/linear/LUSolver.java index 24a09b5d2..6bc6bb34e 100644 --- a/src/java/org/apache/commons/math/linear/LUSolver.java +++ b/src/java/org/apache/commons/math/linear/LUSolver.java @@ -29,17 +29,21 @@ package org.apache.commons.math.linear; public class LUSolver implements DecompositionSolver { /** Serializable version identifier. */ - private static final long serialVersionUID = -8775006035077527661L; + private static final long serialVersionUID = -369589527412301256L; - /** Underlying decomposition. */ - private final LUDecomposition decomposition; + /** Underlying solver. */ + private final DecompositionSolver solver; + + /** Determinant. */ + private final double determinant; /** * Simple constructor. * @param decomposition decomposition to use */ public LUSolver(final LUDecomposition decomposition) { - this.decomposition = decomposition; + this.solver = decomposition.getSolver(); + this.determinant = decomposition.getDeterminant(); } /** Solve the linear equation A × X = B for square matrices A. @@ -52,42 +56,7 @@ public class LUSolver implements DecompositionSolver { */ public double[] solve(final double[] b) throws IllegalArgumentException, InvalidMatrixException { - - final int[] pivot = decomposition.getPivot(); - final int m = pivot.length; - if (b.length != m) { - throw new IllegalArgumentException("constant vector has wrong length"); - } - if (decomposition.isSingular()) { - throw new SingularMatrixException(); - } - - final double[] bp = new double[m]; - - // Apply permutations to b - for (int row = 0; row < m; row++) { - bp[row] = b[pivot[row]]; - } - - // Solve LY = b - final RealMatrix l = decomposition.getL(); - for (int col = 0; col < m; col++) { - for (int i = col + 1; i < m; i++) { - bp[i] -= bp[col] * l.getEntry(i, col); - } - } - - // Solve UX = Y - final RealMatrix u = decomposition.getU(); - for (int col = m - 1; col >= 0; col--) { - bp[col] /= u.getEntry(col, col); - for (int i = 0; i < col; i++) { - bp[i] -= bp[col] * u.getEntry(i, col); - } - } - - return bp; - + return solver.solve(b); } @@ -101,42 +70,7 @@ public class LUSolver implements DecompositionSolver { */ public RealVector solve(final RealVector b) throws IllegalArgumentException, InvalidMatrixException { - - final int[] pivot = decomposition.getPivot(); - final int m = pivot.length; - if (b.getDimension() != m) { - throw new IllegalArgumentException("constant vector has wrong length"); - } - if (decomposition.isSingular()) { - throw new SingularMatrixException(); - } - - final double[] bp = new double[m]; - - // Apply permutations to b - for (int row = 0; row < m; row++) { - bp[row] = b.getEntry(pivot[row]); - } - - // Solve LY = b - final RealMatrix l = decomposition.getL(); - for (int col = 0; col < m; col++) { - for (int i = col + 1; i < m; i++) { - bp[i] -= bp[col] * l.getEntry(i, col); - } - } - - // Solve UX = Y - final RealMatrix u = decomposition.getU(); - for (int col = m - 1; col >= 0; col--) { - bp[col] /= u.getEntry(col, col); - for (int i = 0; i < col; i++) { - bp[i] -= bp[col] * u.getEntry(i, col); - } - } - - return new RealVectorImpl(bp, false); - + return solver.solve(b); } /** Solve the linear equation A × X = B for square matrices A. @@ -149,79 +83,7 @@ public class LUSolver implements DecompositionSolver { */ public RealMatrix solve(final RealMatrix b) throws IllegalArgumentException, InvalidMatrixException { - - final int[] pivot = decomposition.getPivot(); - final int m = pivot.length; - if (b.getRowDimension() != m) { - throw new IllegalArgumentException("Incorrect row dimension"); - } - if (decomposition.isSingular()) { - throw new SingularMatrixException(); - } - - final int nColB = b.getColumnDimension(); - - // Apply permutations to b - final double[][] bp = new double[m][nColB]; - for (int row = 0; row < m; row++) { - final double[] bpRow = bp[row]; - final int pRow = pivot[row]; - for (int col = 0; col < nColB; col++) { - bpRow[col] = b.getEntry(pRow, col); - } - } - - // Solve LY = b - final RealMatrix l = decomposition.getL(); - for (int col = 0; col < m; col++) { - final double[] bpCol = bp[col]; - for (int i = col + 1; i < m; i++) { - final double[] bpI = bp[i]; - final double luICol = l.getEntry(i, col); - for (int j = 0; j < nColB; j++) { - bpI[j] -= bpCol[j] * luICol; - } - } - } - - // Solve UX = Y - final RealMatrix u = decomposition.getU(); - for (int col = m - 1; col >= 0; col--) { - final double[] bpCol = bp[col]; - final double luDiag = u.getEntry(col, col); - for (int j = 0; j < nColB; j++) { - bpCol[j] /= luDiag; - } - for (int i = 0; i < col; i++) { - final double[] bpI = bp[i]; - final double luICol = u.getEntry(i, col); - for (int j = 0; j < nColB; j++) { - bpI[j] -= bpCol[j] * luICol; - } - } - } - - return MatrixUtils.createRealMatrix(bp); - - } - - /** - * Return the determinant of the matrix - * @return determinant of the matrix - * @see #isNonSingular() - */ - public double getDeterminant() { - if (decomposition.isSingular()) { - return 0; - } else { - final int m = decomposition.getPivot().length; - final RealMatrix u = decomposition.getU(); - double determinant = decomposition.evenPermutation() ? 1 : -1; - for (int i = 0; i < m; i++) { - determinant *= u.getEntry(i, i); - } - return determinant; - } + return solver.solve(b); } /** @@ -229,7 +91,7 @@ public class LUSolver implements DecompositionSolver { * @return true if the decomposed matrix is non-singular */ public boolean isNonSingular() { - return !decomposition.isSingular(); + return solver.isNonSingular(); } /** Get the inverse of the decomposed matrix. @@ -238,8 +100,15 @@ public class LUSolver implements DecompositionSolver { */ public RealMatrix getInverse() throws InvalidMatrixException { - final int m = decomposition.getPivot().length; - return solve(MatrixUtils.createRealIdentityMatrix(m)); + return solver.getInverse(); + } + + /** + * Return the determinant of the matrix + * @return determinant of the matrix + */ + public double getDeterminant() { + return determinant; } } diff --git a/src/site/xdoc/userguide/linear.xml b/src/site/xdoc/userguide/linear.xml index f8c60aa24..441cd6744 100644 --- a/src/site/xdoc/userguide/linear.xml +++ b/src/site/xdoc/userguide/linear.xml @@ -67,7 +67,7 @@ System.out.println(p.getRowDimension()); // 2 System.out.println(p.getColumnDimension()); // 2 // Invert p, using LU decomposition -RealMatrix pInverse = new LUSolver(new LUDecompositionImpl(p))).getInverse(); +RealMatrix pInverse = new LUDecompositionImpl(p).getSolver().getInverse();

@@ -115,7 +115,7 @@ RealMatrix pInverse = new LUSolver(new LUDecompositionImpl(p))).getInverse(); RealMatrix coefficients = new RealMatrixImpl(new double[][] { { 2, 3, -2 }, { -1, 7, 6 }, { 4, -3, -5 } }, false); -LUSolver solver = new LUSolver(new LUDecompositionImpl(coefficients)); +DecompositionSolver solver = new LUDecompositionImpl(coefficients).getSolver(); Next create a RealVector array to represent the constant vector B and use solve(RealVector) to solve the system @@ -132,8 +132,9 @@ RealVector solution = solver.solve(constants); for X is such that the residual AX-B has minimal norm. If an exact solution exist (i.e. if for some X the residual AX-B is exactly 0), then this exact solution is also the solution in least square sense. Some solvers like - LUSolver can only find the solution for square matrices and when - the solution is an exact linear solution. Other solvers like QRDecomposition + the one obtained from LUDecomposition can only find the solution + for square matrices and when the solution is an exact linear solution. Other + solvers like the one obtained from QRDecomposition are more versatile and can also find solutions with non-square matrix A or when no exact solution exist (i.e. when the minimal value for AX-B norm is non-null).