You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2008/12/21 23:05:06 UTC
svn commit: r728528 - in /commons/proper/math/trunk/src:
java/org/apache/commons/math/linear/LUDecompositionImpl.java
java/org/apache/commons/math/linear/LUSolver.java
site/xdoc/userguide/linear.xml
Author: luc
Date: Sun Dec 21 14:05:06 2008
New Revision: 728528
URL: http://svn.apache.org/viewvc?rev=728528&view=rev
Log:
added getSolver() into LUDecomposition
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUSolver.java
commons/proper/math/trunk/src/site/xdoc/userguide/linear.xml
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java?rev=728528&r1=728527&r2=728528&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java Sun Dec 21 14:05:06 2008
@@ -17,7 +17,6 @@
package org.apache.commons.math.linear;
-
/**
* Calculates the LUP-decomposition of a square matrix.
* <p>The LUP-decomposition of a matrix A consists of three matrices
@@ -33,7 +32,7 @@
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 @@
/**
* Calculates the LU-decomposition of the given matrix.
- * <p>Calling this constructor is equivalent to first call the no-arguments
- * constructor and then call {@link #decompose(RealMatrix)}.</p>
* @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.
- * <p>Calling this constructor is equivalent to first call the no-arguments
- * constructor and then call {@link #decompose(RealMatrix, double)}.</p>
* @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 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 @@
/** {@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 @@
}
/** {@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.
+ * <p>The A matrix is implicit here. It is </p>
+ * @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));
+ }
+
}
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUSolver.java?rev=728528&r1=728527&r2=728528&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUSolver.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/linear/LUSolver.java Sun Dec 21 14:05:06 2008
@@ -29,17 +29,21 @@
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 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 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 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 @@
* @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 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;
}
}
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/linear.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/linear.xml?rev=728528&r1=728527&r2=728528&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/linear.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/linear.xml Sun Dec 21 14:05:06 2008
@@ -67,7 +67,7 @@
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();
</source>
</p>
<p>
@@ -115,7 +115,7 @@
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();
</source>
Next create a <code>RealVector</code> array to represent the constant
vector B and use <code>solve(RealVector)</code> to solve the system
@@ -132,8 +132,9 @@
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
- <code>LUSolver</code> can only find the solution for square matrices and when
- the solution is an exact linear solution. Other solvers like <code>QRDecomposition</code>
+ the one obtained from <code>LUDecomposition</code> can only find the solution
+ for square matrices and when the solution is an exact linear solution. Other
+ solvers like the one obtained from <code>QRDecomposition</code>
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).
</p>