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
- * [EMAIL PROTECTED] QRDecomposition} or [EMAIL PROTECTED] LUDecomposition}.
All these algorithms
- * decompose an A matrix has a product of several specific matrices from
- * which they can solve A × 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 × X is the
- * closest possible to B.</p>
+ * [EMAIL PROTECTED] QRDecomposition}, [EMAIL PROTECTED] LUDecomposition} or
[EMAIL PROTECTED]
+ * SingularValueDecomposition}. All these algorithms decompose an A matrix has
a
+ * product of several specific matrices from which they can solve A × X
= B
+ * in least squares sense: they find X such that ||A × X - B|| is
minimal.</p>
+ * <p>Some solvers like [EMAIL PROTECTED] LUDecomposition} can only find the
solution for
+ * square matrices and when the solution is an exact linear solution, i.e. when
+ * ||A × 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 [EMAIL PROTECTED] LUDecomposition})
+ */
+ void decompose(RealMatrix matrix)
+ throws InvalidMatrixException;
+
/** Solve the linear equation A × 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 [EMAIL PROTECTED]
#decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
- * @throws IllegalArgumentException if matrices dimensions don't match
- * @throws InvalidMatrixException if decomposed matrix is singular
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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 × 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 [EMAIL PROTECTED]
#decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
- * @throws IllegalArgumentException if matrices dimensions don't match
- * @throws InvalidMatrixException if decomposed matrix is singular
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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 × 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 [EMAIL PROTECTED]
#decompose(RealMatrix)}.</p>
* @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 IllegalArgumentException if matrices dimensions don't match
- * @throws InvalidMatrixException if decomposed matrix is singular
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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×A = L×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 [EMAIL
PROTECTED]
- * #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 [EMAIL PROTECTED]
#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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
#getPivot()
* pivot permutation vector}.</p>
* @return the P rows permutation matrix (or null if decomposed matrix is
singular)
+ * @exception IllegalStateException if [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED] #decompose(RealMatrix)} or
+ * [EMAIL PROTECTED] #decompose(RealMatrix, double)} <strong>must</strong>
be called
+ * before any of the [EMAIL PROTECTED] #getP()}, [EMAIL PROTECTED]
#getPivot()}, [EMAIL PROTECTED] #getL()},
+ * [EMAIL PROTECTED] #getU()}, [EMAIL PROTECTED] #getDeterminant()},
[EMAIL PROTECTED] #isNonSingular()},
+ * [EMAIL PROTECTED] #solve(double[])}, [EMAIL PROTECTED]
#solve(RealMatrix)}, [EMAIL PROTECTED] #solve(RealVector)}
+ * or [EMAIL PROTECTED] #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 [EMAIL PROTECTED] #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 [EMAIL PROTECTED] #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);
+ }
+
+ /** [EMAIL PROTECTED] */
+ public void decompose(RealMatrix matrix)
+ throws InvalidMatrixException {
+ decompose(matrix, DEFAULT_TOO_SMALL);
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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;
+ }
+ }
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- public int[] getPivot() {
+ public int[] getPivot()
+ throws IllegalStateException {
+ checkDecomposed();
return pivot.clone();
}
/** [EMAIL PROTECTED] */
- public boolean isNonSingular() {
+ public boolean isNonSingular()
+ throws IllegalStateException {
+ checkDecomposed();
return !singular;
}
/** [EMAIL PROTECTED] */
- 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 @@
/** [EMAIL PROTECTED] */
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 @@
/** [EMAIL PROTECTED] */
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);
}
/** [EMAIL PROTECTED] */
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 [EMAIL PROTECTED] #decompose(RealMatrix)} or [EMAIL
PROTECTED]
+ * #decompose(RealMatrix, double) has been called.
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED]
+ * 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 [EMAIL PROTECTED] #decompose(RealMatrix)}
<strong>must</strong> be called
+ * before any of the [EMAIL PROTECTED] #getQ()}, [EMAIL PROTECTED]
#getR()}, [EMAIL PROTECTED] #getH()},
+ * [EMAIL PROTECTED] #isFullRank()}, [EMAIL PROTECTED] #solve(double[])},
[EMAIL PROTECTED] #solve(RealMatrix)},
+ * [EMAIL PROTECTED] #solve(RealVector)} or [EMAIL PROTECTED]
#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 [EMAIL PROTECTED] #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);
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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 @@
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- 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 @@
}
/** [EMAIL PROTECTED] */
- public boolean isFullRank() {
+ public boolean isFullRank()
+ throws IllegalStateException {
+
+ checkDecomposed();
+
for (double diag : rDiag) {
if (diag == 0) {
return false;
}
}
return true;
+
}
/** [EMAIL PROTECTED] */
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;
}
/** [EMAIL PROTECTED] */
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 × X = B
* @return a vector X that minimizes the two norm of A × X - B
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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);
}
/** [EMAIL PROTECTED] */
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 [EMAIL PROTECTED] #decompose(RealMatrix)} has been called.
+ * @exception IllegalStateException if [EMAIL PROTECTED]
#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");
+ }
+ }
+
}
