commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r690803 - in /commons/proper/math/branches/MATH_2_0/src: java/org/apache/commons/math/linear/ test/org/apache/commons/math/linear/
Date Sun, 31 Aug 2008 22:28:56 GMT
Author: luc
Date: Sun Aug 31 15:28:55 2008
New Revision: 690803

URL: http://svn.apache.org/viewvc?rev=690803&view=rev
Log:
added a decompose method to the base DecompositionSolver interface
to allow a solver to decompose a matrix after construction.

Modified:
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/DecompositionSolver.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecomposition.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/DecompositionSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/DecompositionSolver.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/DecompositionSolver.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/DecompositionSolver.java
Sun Aug 31 15:28:55 2008
@@ -22,48 +22,68 @@
 /**
  * A base interface to decomposition algorithms that can solve A × X = B.
  * <p>This interface is the common base of decomposition algorithms like
- * {@link QRDecomposition} or {@link LUDecomposition}. All these algorithms
- * decompose an A matrix has a product of several specific matrices from
- * which they can solve A &times; X = B.</p>
- * <p>Depending on the solver, the solution is either an exact linear solution
- * or a least squares solution. When an exact linear solution exist, both the
- * linear and the least squares solution are equal. When no exact linear solution
- * exist, a least square solution gives an X which such that A &times; X is the
- * closest possible to B.</p>
+ * {@link QRDecomposition}, {@link LUDecomposition} or {@link
+ * SingularValueDecomposition}. All these algorithms decompose an A matrix has a
+ * product of several specific matrices from which they can solve A &times; X = B
+ * in least squares sense: they find X such that ||A &times; X - B|| is minimal.</p>
+ * <p>Some solvers like {@link LUDecomposition} can only find the solution for
+ * square matrices and when the solution is an exact linear solution, i.e. when
+ * ||A &times; X - B|| is exactly 0. Other solvers can also find solutions
+ * with non-square matrix A and with non-null minimal norm. If an exact linear
+ * solution exists it is also the minimal norm solution.</p>
  *   
  * @version $Revision$ $Date$
  * @since 2.0
  */
 public interface DecompositionSolver extends Serializable {
 
+    /**
+     * Decompose a matrix.
+     * @param matrix
+     * @exception InvalidMatrixException if matrix does not fulfill
+     * the decomposition requirements (for example non-square matrix
+     * for {@link LUDecomposition})
+     */
+    void decompose(RealMatrix matrix)
+      throws InvalidMatrixException;
+
     /** Solve the linear equation A &times; X = B.
-     * <p>The A matrix is implicit here. It is </p>
+     * <p>The A matrix is implicit here. It <strong>must</strong> have
+     * already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
      * @param b right-hand side of the equation A &times; X = B
      * @return a vector X that minimizes the two norm of A &times; X - B
-     * @throws IllegalArgumentException if matrices dimensions don't match
-     * @throws InvalidMatrixException if decomposed matrix is singular
+     * @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
      */
     double[] solve(double[] b)
-      throws IllegalArgumentException, InvalidMatrixException;
+      throws IllegalStateException, IllegalArgumentException, InvalidMatrixException;
 
     /** Solve the linear equation A &times; X = B.
-     * <p>The A matrix is implicit here. It is </p>
+     * <p>The A matrix is implicit here. It <strong>must</strong> have
+     * already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
      * @param b right-hand side of the equation A &times; X = B
      * @return a vector X that minimizes the two norm of A &times; X - B
-     * @throws IllegalArgumentException if matrices dimensions don't match
-     * @throws InvalidMatrixException if decomposed matrix is singular
+     * @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
      */
     RealVector solve(RealVector b)
-      throws IllegalArgumentException, InvalidMatrixException;
+      throws IllegalStateException, IllegalArgumentException, InvalidMatrixException;
 
     /** Solve the linear equation A &times; X = B.
-     * <p>The A matrix is implicit here. It is </p>
+     * <p>The A matrix is implicit here. It <strong>must</strong> have
+     * already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
      * @param b right-hand side of the equation A &times; X = B
      * @return a matrix X that minimizes the two norm of A &times; X - B
-     * @throws IllegalArgumentException if matrices dimensions don't match
-     * @throws InvalidMatrixException if decomposed matrix is singular
+     * @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
      */
     RealMatrix solve(RealMatrix b)
-      throws IllegalArgumentException, InvalidMatrixException;
+      throws IllegalStateException, IllegalArgumentException, InvalidMatrixException;
 
 }

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
Sun Aug 31 15:28:55 2008
@@ -24,10 +24,14 @@
  * such that P&times;A = L&times;U. P is a rows permutation matrix that is used
  * to rearrange the rows of A before so that it can be decomposed. L is a lower
  * triangular matrix with unit diagonal terms and U is an upper triangular matrix.</p>
- * <p>This interface is similar to the class with similar name from the now defunct
+ * <p>This interface is based on the class with similar name from the now defunct
  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
- * exception of the <code>det</code> method which has been renamed as {@link
- * #getDeterminant() getDeterminant}.</p>
+ * following changes:</p>
+ * <ul>
+ *   <li>several signatures have been added for the <code>solve</code>
methods (in the superinterface),</code>
+ *   <li>a <code>decompose</code> method has been added (in the superinterface),</code>
+ *   <li>the <code>det</code> method has been renamed as {@link #getDeterminant()
getDeterminant}.</li>
+ * </ul>
  *   
  * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
  * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
@@ -37,18 +41,37 @@
 public interface LUDecomposition extends DecompositionSolver {
 
     /**
+     * Computes a new
+     * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
+     * LU decomposition</a> for this matrix, storing the result for use by other methods.
+     * <p>
+     * <strong>Implementation Note</strong>:<br>
+     * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
+     * Crout's algorithm</a>, with partial pivoting.</p>
+     * @param matrix The matrix to decompose.
+     * @param singularityThreshold threshold (based on partial row norm)
+     * under which a matrix is considered singular
+     * @exception InvalidMatrixException if matrix is not square
+     */
+    void decompose(RealMatrix matrix, double singularityThreshold);
+
+    /**
      * Returns the matrix L of the decomposition. 
      * <p>L is an lower-triangular matrix</p>
      * @return the L matrix (or null if decomposed matrix is singular)
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    RealMatrix getL();
+    RealMatrix getL() throws IllegalStateException;
 
     /**
      * Returns the matrix U of the decomposition. 
      * <p>U is an upper-triangular matrix</p>
      * @return the U matrix (or null if decomposed matrix is singular)
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    RealMatrix getU();
+    RealMatrix getU() throws IllegalStateException;
 
     /**
      * Returns the P rows permutation matrix.
@@ -57,29 +80,37 @@
      * <p>The positions of the 1 elements are given by the {@link #getPivot()
      * pivot permutation vector}.</p>
      * @return the P rows permutation matrix (or null if decomposed matrix is singular)
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      * @see #getPivot()
      */
-    RealMatrix getP();
+    RealMatrix getP() throws IllegalStateException;
 
     /**
      * Returns the pivot permutation vector.
      * @return the pivot permutation vector
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      * @see #getPermutation()
      */
-    int[] getPivot();
+    int[] getPivot() throws IllegalStateException;
 
     /**
      * Check if the decomposed matrix is non-singular.
      * @return true if the decomposed matrix is non-singular
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      * @see #getDeterminant()
      */
-    boolean isNonSingular();
+    boolean isNonSingular() throws IllegalStateException;
 
     /**
      * Return the determinant of the matrix
      * @return determinant of the matrix
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      * @see #isNonSingular()
      */
-    double getDeterminant();
+    double getDeterminant() throws IllegalStateException;
 
 }

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
Sun Aug 31 15:28:55 2008
@@ -32,19 +32,13 @@
 public class LUDecompositionImpl implements LUDecomposition {
 
     /** Serializable version identifier. */
-    private static final long serialVersionUID = -1606789599960880183L;
-
-    /** Bound to determine effective singularity in LU decomposition */
-    private final double singularityThreshold;
-
-    /** Size of the matrix. */
-    private final int m;
+    private static final long serialVersionUID = -9052751605297201067L;
 
     /** Entries of LU decomposition. */
-    private final double lu[][];
+    private double lu[][];
 
     /** Pivot permutation associated with LU decomposition */
-    private final int[] pivot;
+    private int[] pivot;
 
     /** Parity of the permutation associated with the LU decomposition */
     private int parity;
@@ -65,19 +59,35 @@
     private static final double DEFAULT_TOO_SMALL = 10E-12;
 
     /**
+     * Build a new instance.
+     * <p>Note that either {@link #decompose(RealMatrix)} or
+     * {@link #decompose(RealMatrix, double)} <strong>must</strong> be called
+     * before any of the {@link #getP()}, {@link #getPivot()}, {@link #getL()},
+     * {@link #getU()}, {@link #getDeterminant()}, {@link #isNonSingular()},
+     * {@link #solve(double[])}, {@link #solve(RealMatrix)}, {@link #solve(RealVector)}
+     * or {@link #solve(RealVectorImpl)} methods can be called.</p>
+     * @see #decompose(RealMatrix)
+     * @see #decompose(RealMatrix, double)
+     */
+    public LUDecompositionImpl() {
+    }
+
+    /**
      * 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 {
-        this(matrix, DEFAULT_TOO_SMALL);
+        decompose(matrix);
     }
 
     /**
      * 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
@@ -85,25 +95,103 @@
      */
     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 InvalidMatrixException("LU decomposition requires that the matrix be
square.");
         }
-        this.singularityThreshold = singularityThreshold;
-        m = matrix.getColumnDimension();
+        final int m = matrix.getColumnDimension();
         lu = matrix.getData();
         pivot = new int[m];
         cachedL = null;
         cachedU = null;
         cachedP = null;
 
-        // perform decomposition
-        luDecompose();
+        // Initialize permutation array and parity
+        for (int row = 0; row < m; row++) {
+            pivot[row] = row;
+        }
+        parity = 1;
+        singular = false;
+
+        // Loop over columns
+        for (int col = 0; col < m; col++) {
+
+            double sum = 0;
+
+            // upper
+            for (int row = 0; row < col; row++) {
+                final double[] luRow = lu[row];
+                sum = luRow[col];
+                for (int i = 0; i < row; i++) {
+                    sum -= luRow[i] * lu[i][col];
+                }
+                luRow[col] = sum;
+            }
+
+            // lower
+            int max = col; // permutation row
+            double largest = Double.NEGATIVE_INFINITY;
+            for (int row = col; row < m; row++) {
+                final double[] luRow = lu[row];
+                sum = luRow[col];
+                for (int i = 0; i < col; i++) {
+                    sum -= luRow[i] * lu[i][col];
+                }
+                luRow[col] = sum;
+
+                // maintain best permutation choice
+                if (Math.abs(sum) > largest) {
+                    largest = Math.abs(sum);
+                    max = row;
+                }
+            }
+
+            // Singularity check
+            if (Math.abs(lu[max][col]) < singularityThreshold) {
+                singular = true;
+                return;
+            }
+
+            // Pivot if necessary
+            if (max != col) {
+                double tmp = 0;
+                for (int i = 0; i < m; i++) {
+                    tmp = lu[max][i];
+                    lu[max][i] = lu[col][i];
+                    lu[col][i] = tmp;
+                }
+                int temp = pivot[max];
+                pivot[max] = pivot[col];
+                pivot[col] = temp;
+                parity = -parity;
+            }
+
+            // Divide the lower elements by the "winning" diagonal elt.
+            final double luDiag = lu[col][col];
+            for (int row = col + 1; row < m; row++) {
+                lu[row][col] /= luDiag;
+            }
+        }
 
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getL() {
+    public RealMatrix getL()
+        throws IllegalStateException {
+        checkDecomposed();
         if ((cachedL == null) && !singular) {
+            final int m = pivot.length;
             final double[][] lData = new double[m][m];
             for (int i = 0; i < m; ++i) {
                 System.arraycopy(lu[i], 0, lData[i], 0, i);
@@ -115,8 +203,11 @@
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getU() {
+    public RealMatrix getU()
+        throws IllegalStateException {
+        checkDecomposed();
         if ((cachedU == null) && !singular) {
+            final int m = pivot.length;
             final double[][] uData = new double[m][m];
             for (int i = 0; i < m; ++i) {
                 System.arraycopy(lu[i], i, uData[i], i, m - i);
@@ -127,8 +218,11 @@
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getP() {
+    public RealMatrix getP()
+        throws IllegalStateException {
+        checkDecomposed();
         if ((cachedP == null) && !singular) {
+            final int m = pivot.length;
             final double[][] pData = new double[m][m];
             for (int i = 0; i < m; ++i) {
                 pData[i][pivot[i]] = 1.0;
@@ -139,20 +233,27 @@
     }
 
     /** {@inheritDoc} */
-    public int[] getPivot() {
+    public int[] getPivot()
+        throws IllegalStateException {
+        checkDecomposed();
         return pivot.clone();
     }
 
     /** {@inheritDoc} */
-    public boolean isNonSingular() {
+    public boolean isNonSingular()
+        throws IllegalStateException {
+        checkDecomposed();
         return !singular;
     }
 
     /** {@inheritDoc} */
-    public double getDeterminant() {
+    public double getDeterminant()
+        throws IllegalStateException {
+        checkDecomposed();
         if (singular) {
             return 0;
         } else {
+            final int m = pivot.length;
             double determinant = parity;
             for (int i = 0; i < m; i++) {
                 determinant *= lu[i][i];
@@ -163,8 +264,10 @@
 
     /** {@inheritDoc} */
     public double[] solve(double[] b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
 
+        checkDecomposed();
+        final int m = pivot.length;
         if (b.length != m) {
             throw new IllegalArgumentException("constant vector has wrong length");
         }
@@ -200,11 +303,13 @@
 
     /** {@inheritDoc} */
     public RealVector solve(RealVector b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
         try {
             return solve((RealVectorImpl) b);
         } catch (ClassCastException cce) {
 
+            checkDecomposed();
+            final int m = pivot.length;
             if (b.getDimension() != m) {
                 throw new IllegalArgumentException("constant vector has wrong length");
             }
@@ -247,13 +352,16 @@
      * @throws InvalidMatrixException if decomposed matrix is singular
      */
     public RealVectorImpl solve(RealVectorImpl b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
         return new RealVectorImpl(solve(b.getDataRef()), false);
     }
 
     /** {@inheritDoc} */
     public RealMatrix solve(RealMatrix b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
+
+        checkDecomposed();
+        final int m = pivot.length;
         if (b.getRowDimension() != m) {
             throw new IllegalArgumentException("Incorrect row dimension");
         }
@@ -306,88 +414,15 @@
     }
 
     /**
-     * Computes a new
-     * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
-     * LU decomposition</a> for this matrix, storing the result for use by other methods.
-     * <p>
-     * <strong>Implementation Note</strong>:<br>
-     * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
-     * Crout's algorithm</a>, with partial pivoting.</p>
-     * <p>
-     * <strong>Usage Note</strong>:<br>
-     * This method should rarely be invoked directly. Its only use is
-     * to force recomputation of the LU decomposition when changes have been
-     * made to the underlying data using direct array references. Changes
-     * made using setXxx methods will trigger recomputation when needed
-     * automatically.</p>
+     * Check if either {@link #decompose(RealMatrix)} or {@link
+     * #decompose(RealMatrix, double) has been called.
+     * @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
+     * has not been called
      */
-    private void luDecompose() {
-
-        // Initialize permutation array and parity
-        for (int row = 0; row < m; row++) {
-            pivot[row] = row;
-        }
-        parity = 1;
-        singular = false;
-
-        // Loop over columns
-        for (int col = 0; col < m; col++) {
-
-            double sum = 0;
-
-            // upper
-            for (int row = 0; row < col; row++) {
-                final double[] luRow = lu[row];
-                sum = luRow[col];
-                for (int i = 0; i < row; i++) {
-                    sum -= luRow[i] * lu[i][col];
-                }
-                luRow[col] = sum;
-            }
-
-            // lower
-            int max = col; // permutation row
-            double largest = Double.NEGATIVE_INFINITY;
-            for (int row = col; row < m; row++) {
-                final double[] luRow = lu[row];
-                sum = luRow[col];
-                for (int i = 0; i < col; i++) {
-                    sum -= luRow[i] * lu[i][col];
-                }
-                luRow[col] = sum;
-
-                // maintain best permutation choice
-                if (Math.abs(sum) > largest) {
-                    largest = Math.abs(sum);
-                    max = row;
-                }
-            }
-
-            // Singularity check
-            if (Math.abs(lu[max][col]) < singularityThreshold) {
-                singular = true;
-                return;
-            }
-
-            // Pivot if necessary
-            if (max != col) {
-                double tmp = 0;
-                for (int i = 0; i < m; i++) {
-                    tmp = lu[max][i];
-                    lu[max][i] = lu[col][i];
-                    lu[col][i] = tmp;
-                }
-                int temp = pivot[max];
-                pivot[max] = pivot[col];
-                pivot[col] = temp;
-                parity = -parity;
-            }
-
-            // Divide the lower elements by the "winning" diagonal elt.
-            final double luDiag = lu[col][col];
-            for (int row = col + 1; row < m; row++) {
-                lu[row][col] /= luDiag;
-            }
+    private void checkDecomposed()
+        throws IllegalStateException {
+        if (lu == null) {
+            throw new IllegalStateException("no matrix have been decomposed yet");
         }
     }
 

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecomposition.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecomposition.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecomposition.java
Sun Aug 31 15:28:55 2008
@@ -20,8 +20,13 @@
 /**
  * An interface to classes that implement a algorithm to calculate the 
  * QR-decomposition of a real matrix.
- * <p>This interface is similar to the class with similar name from the now defunct
- * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
+ * <p>This interface is based on the class with similar name from the now defunct
+ * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
+ * following changes:</p>
+ * <ul>
+ *   <li>several signatures have been added for the <code>solve</code>
methods (in the superinterface),</code>
+ *   <li>a <code>decompose</code> method has been added (in the superinterface),</code>
+ * </ul>
  *   
  * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
  * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
@@ -34,15 +39,19 @@
      * Returns the matrix R of the decomposition. 
      * <p>R is an upper-triangular matrix</p>
      * @return the R matrix
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    RealMatrix getR();
+    RealMatrix getR() throws IllegalStateException;
 
     /**
      * Returns the matrix Q of the decomposition.
      * <p>Q is an orthogonal matrix</p>
      * @return the Q matrix
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    RealMatrix getQ();
+    RealMatrix getQ() throws IllegalStateException;
 
     /**
      * Returns the Householder reflector vectors.
@@ -50,13 +59,17 @@
      * each successive Householder reflector vector. This matrix is used
      * to compute Q.</p>
      * @return a matrix containing the Householder reflector vectors
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    RealMatrix getH();
+    RealMatrix getH() throws IllegalStateException;
 
     /**
      * Check if the decomposed matrix is full rank.
      * @return true if the decomposed matrix is full rank
+     * @exception IllegalStateException if {@link
+     * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
      */
-    boolean isFullRank();
+    boolean isFullRank() throws IllegalStateException;
 
 }

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
Sun Aug 31 15:28:55 2008
@@ -34,7 +34,7 @@
 public class QRDecompositionImpl implements QRDecomposition {
 
     /** Serializable version identifier. */
-    private static final long serialVersionUID = 3965943878043764074L;
+    private static final long serialVersionUID = 7125583145349720380L;
 
     /**
      * A packed representation of the QR decomposition. The elements above the 
@@ -42,12 +42,10 @@
      * are the Householder reflector vectors of which an explicit form of Q can
      * be calculated. 
      */
-    private final double[][] qr;
+    private double[][] qr;
 
-    /**
-     * The diagonal elements of R.
-     */
-    private final double[] rDiag;
+    /** The diagonal elements of R. */
+    private double[] rDiag;
 
     /** Cached value of Q. */
     private RealMatrix cachedQ;
@@ -59,24 +57,33 @@
     private RealMatrix cachedH;
 
     /**
-     * The row dimension of the given matrix. The size of Q will be m x m, the 
-     * size of R will be m x n. 
+     * Build a new instance.
+     * <p>Note that {@link #decompose(RealMatrix)} <strong>must</strong>
be called
+     * before any of the {@link #getQ()}, {@link #getR()}, {@link #getH()},
+     * {@link #isFullRank()}, {@link #solve(double[])}, {@link #solve(RealMatrix)},
+     * {@link #solve(RealVector)} or {@link #solve(RealVectorImpl)} methods can be
+     * called.</p>
+     * @see #decompose(RealMatrix)
      */
-    private final int m;
-
-    /**
-     * The column dimension of the given matrix. The size of R will be m x n. 
-     */
-    private final int n;
+    public QRDecompositionImpl() {
+    }
 
     /**
-     * Calculates the QR decomposition of the given matrix. 
-     * 
+     * Calculates the QR-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 QRDecompositionImpl(RealMatrix matrix) {
-        m = matrix.getRowDimension();
-        n = matrix.getColumnDimension();
+    public QRDecompositionImpl(RealMatrix matrix)
+        throws InvalidMatrixException {
+        decompose(matrix);
+    }
+
+    /** {@inheritDoc} */
+    public void decompose(RealMatrix matrix) {
+        final int m = matrix.getRowDimension();
+        final int n = matrix.getColumnDimension();
         qr = matrix.getData();
         rDiag = new double[n];
         cachedQ = null;
@@ -147,11 +154,16 @@
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getR() {
+    public RealMatrix getR()
+        throws IllegalStateException {
 
         if (cachedR == null) {
 
+            checkDecomposed();
+
             // R is supposed to be m x n
+            final int m = qr.length;
+            final int n = qr[0].length;
             double[][] r = new double[m][n];
 
             // copy the diagonal from rDiag and the upper triangle of qr
@@ -172,11 +184,16 @@
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getQ() {
+    public RealMatrix getQ()
+        throws IllegalStateException {
 
         if (cachedQ == null) {
 
+            checkDecomposed();
+
             // Q is supposed to be m x m
+            final int m = qr.length;
+            final int n = qr[0].length;
             double[][] Q = new double[m][m];
 
             /* 
@@ -216,9 +233,15 @@
     }
 
     /** {@inheritDoc} */
-    public RealMatrix getH() {
+    public RealMatrix getH()
+        throws IllegalStateException {
+
         if (cachedH == null) {
 
+            checkDecomposed();
+
+            final int m = qr.length;
+            final int n = qr[0].length;
             double[][] hData = new double[m][n];
             for (int i = 0; i < m; ++i) {
                 for (int j = 0; j < Math.min(i + 1, n); ++j) {
@@ -237,64 +260,74 @@
     }
 
     /** {@inheritDoc} */
-    public boolean isFullRank() {
+    public boolean isFullRank()
+        throws IllegalStateException {
+
+        checkDecomposed();
+
         for (double diag : rDiag) {
             if (diag == 0) {
                 return false;
             }
         }
         return true;
+
     }
 
     /** {@inheritDoc} */
     public double[] solve(double[] b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
  
-            if (b.length != m) {
-                throw new IllegalArgumentException("Incorrect row dimension");
-            }
-            if (!isFullRank()) {
-                throw new InvalidMatrixException("Matrix is rank-deficient");
-            }
+        checkDecomposed();
 
-            final double[] x = new double[n];
-            final double[] y = b.clone();
+        final int m = qr.length;
+        final int n = qr[0].length;
+        if (b.length != m) {
+            throw new IllegalArgumentException("Incorrect row dimension");
+        }
+        if (!isFullRank()) {
+            throw new InvalidMatrixException("Matrix is rank-deficient");
+        }
 
-            // apply Householder transforms to solve Q.y = b
-            for (int minor = 0; minor < Math.min(m, n); minor++) {
+        final double[] x = new double[n];
+        final double[] y = b.clone();
 
-                double dotProduct = 0;
-                for (int row = minor; row < m; row++) {
-                    dotProduct += y[row] * qr[row][minor];
-                }
-                dotProduct /= rDiag[minor] * qr[minor][minor];
+        // apply Householder transforms to solve Q.y = b
+        for (int minor = 0; minor < Math.min(m, n); minor++) {
 
-                for (int row = minor; row < m; row++) {
-                    y[row] += dotProduct * qr[row][minor];
-                }
+            double dotProduct = 0;
+            for (int row = minor; row < m; row++) {
+                dotProduct += y[row] * qr[row][minor];
+            }
+            dotProduct /= rDiag[minor] * qr[minor][minor];
 
+            for (int row = minor; row < m; row++) {
+                y[row] += dotProduct * qr[row][minor];
             }
 
-            // solve triangular system R.x = y
-            for (int row = n - 1; row >= 0; --row) {
-                y[row] /= rDiag[row];
-                final double yRow = y[row];
-                x[row] = yRow;
-                for (int i = 0; i < row; i++) {
-                    y[i] -= yRow * qr[i][row];
-                }
+        }
+
+        // solve triangular system R.x = y
+        for (int row = n - 1; row >= 0; --row) {
+            y[row] /= rDiag[row];
+            final double yRow = y[row];
+            x[row] = yRow;
+            for (int i = 0; i < row; i++) {
+                y[i] -= yRow * qr[i][row];
             }
+        }
 
-            return x;
+        return x;
 
     }
 
     /** {@inheritDoc} */
     public RealVector solve(RealVector b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
         try {
             return solve((RealVectorImpl) b);
         } catch (ClassCastException cce) {
+            checkDecomposed();
             return new RealVectorImpl(solve(b.getData()), false);
         }
     }
@@ -303,18 +336,24 @@
      * <p>The A matrix is implicit here. It is </p>
      * @param b right-hand side of the equation A &times; X = B
      * @return a vector X that minimizes the two norm of A &times; X - B
+     * @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
+     * has not been called
      * @throws IllegalArgumentException if matrices dimensions don't match
      * @throws InvalidMatrixException if decomposed matrix is singular
      */
     public RealVectorImpl solve(RealVectorImpl b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
         return new RealVectorImpl(solve(b.getDataRef()), false);
     }
 
     /** {@inheritDoc} */
     public RealMatrix solve(RealMatrix b)
-        throws IllegalArgumentException, InvalidMatrixException {
+        throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
 
+        checkDecomposed();
+
+        final int m = qr.length;
+        final int n = qr[0].length;
         if (b.getRowDimension() != m) {
             throw new IllegalArgumentException("Incorrect row dimension");
         }
@@ -364,4 +403,16 @@
 
     }
 
+    /**
+     * Check if {@link #decompose(RealMatrix)} has been called.
+     * @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
+     * has not been called
+     */
+    private void checkDecomposed()
+        throws IllegalStateException {
+        if (qr == null) {
+            throw new IllegalStateException("no matrix have been decomposed yet");
+        }
+    }
+
 }

Modified: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
Sun Aug 31 15:28:55 2008
@@ -88,6 +88,29 @@
         }
     }
 
+    /** test no call to decompose */
+    public void testNoDecompose() {
+        try {
+            new LUDecompositionImpl().getPivot();
+            fail("an exception should have been caught");
+        } catch (IllegalStateException ise) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    /** test threshold impact */
+    public void testThreshold() {
+        final RealMatrix matrix = new RealMatrixImpl(new double[][] {
+                                                       { 1.0, 2.0, 3.0},
+                                                       { 2.0, 5.0, 3.0},
+                                                       { 4.000001, 9.0, 9.0}
+                                                     }, false);
+        assertFalse(new LUDecompositionImpl(matrix, 1.0e-5).isNonSingular());
+        assertTrue(new LUDecompositionImpl(matrix, 1.0e-10).isNonSingular());
+    }
+
     /** test PA = LU */
     public void testPAEqualLU() {
         RealMatrix matrix = new RealMatrixImpl(testData, false);

Modified: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java?rev=690803&r1=690802&r2=690803&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
Sun Aug 31 15:28:55 2008
@@ -349,4 +349,16 @@
         
     }
 
+    /** test no call to decompose */
+    public void testNoDecompose() {
+        try {
+            new QRDecompositionImpl().isFullRank();
+            fail("an exception should have been caught");
+        } catch (IllegalStateException ise) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
 }



Mime
View raw message