Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added a decompose method to the base DecompositionSolver interface 

to allow a solver to decompose a matrix after construction.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@690803 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-08-31 22:28:55 +00:00
parent c1f2e60e81
commit 481a091105
7 changed files with 518 additions and 333 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
src/java/org/apache/commons/math/linear/DecompositionSolver.java
View File

@ -22,48 +22,68 @@ import java.io.Serializable;
/**
* 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;
}

49
src/java/org/apache/commons/math/linear/LUDecomposition.java
View File

@ -24,10 +24,14 @@ package org.apache.commons.math.linear;
* 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>
@ -36,19 +40,38 @@ package org.apache.commons.math.linear;
*/
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 @@ public interface LUDecomposition extends DecompositionSolver {
* <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;
}

517
src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
View File

@ -32,19 +32,13 @@ package org.apache.commons.math.linear;
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;
@ -64,265 +58,65 @@ public class LUDecompositionImpl implements LUDecomposition {
/** Default bound to determine effective singularity in LU decomposition */
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
* @exception InvalidMatrixException if matrix is not square
*/
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();
}
/** {@inheritDoc} */
public RealMatrix getL() {
if ((cachedL == null) && !singular) {
final double[][] lData = new double[m][m];
for (int i = 0; i < m; ++i) {
System.arraycopy(lu[i], 0, lData[i], 0, i);
lData[i][i] = 1.0;
}
cachedL = new RealMatrixImpl(lData, false);
}
return cachedL;
}
/** {@inheritDoc} */
public RealMatrix getU() {
if ((cachedU == null) && !singular) {
final double[][] uData = new double[m][m];
for (int i = 0; i < m; ++i) {
System.arraycopy(lu[i], i, uData[i], i, m - i);
}
cachedU = new RealMatrixImpl(uData, false);
}
return cachedU;
}
/** {@inheritDoc} */
public RealMatrix getP() {
if ((cachedP == null) && !singular) {
final double[][] pData = new double[m][m];
for (int i = 0; i < m; ++i) {
pData[i][pivot[i]] = 1.0;
}
cachedP = new RealMatrixImpl(pData, false);
}
return cachedP;
}
/** {@inheritDoc} */
public int[] getPivot() {
return pivot.clone();
}
/** {@inheritDoc} */
public boolean isNonSingular() {
return !singular;
}
/** {@inheritDoc} */
public double getDeterminant() {
if (singular) {
return 0;
} else {
double determinant = parity;
for (int i = 0; i < m; i++) {
determinant *= lu[i][i];
}
return determinant;
}
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalArgumentException, InvalidMatrixException {
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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 IllegalArgumentException, InvalidMatrixException {
try {
return solve((RealVectorImpl) b);
} catch (ClassCastException cce) {
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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 &times; X = B.
* <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 such that A &times; X = B
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealVectorImpl solve(RealVectorImpl b)
throws IllegalArgumentException, InvalidMatrixException {
return new RealVectorImpl(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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);
}
/**
* 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>
*/
private void luDecompose() {
// Initialize permutation array and parity
for (int row = 0; row < m; row++) {
pivot[row] = row;
@ -389,6 +183,247 @@ public class LUDecompositionImpl implements LUDecomposition {
lu[row][col] /= luDiag;
}
}
}
/** {@inheritDoc} */
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);
lData[i][i] = 1.0;
}
cachedL = new RealMatrixImpl(lData, false);
}
return cachedL;
}
/** {@inheritDoc} */
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);
}
cachedU = new RealMatrixImpl(uData, false);
}
return cachedU;
}
/** {@inheritDoc} */
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;
}
cachedP = new RealMatrixImpl(pData, false);
}
return cachedP;
}
/** {@inheritDoc} */
public int[] getPivot()
throws IllegalStateException {
checkDecomposed();
return pivot.clone();
}
/** {@inheritDoc} */
public boolean isNonSingular()
throws IllegalStateException {
checkDecomposed();
return !singular;
}
/** {@inheritDoc} */
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];
}
return determinant;
}
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalStateException, IllegalArgumentException, InvalidMatrixException {
checkDecomposed();
final int m = pivot.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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) {
checkDecomposed();
final int m = pivot.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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 &times; X = B.
* <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 such that A &times; X = B
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws 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 {
checkDecomposed();
final int m = pivot.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (singular) {
throw new InvalidMatrixException("Matrix is singular.");
}
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);
}
/**
* 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 checkDecomposed()
throws IllegalStateException {
if (lu == null) {
throw new IllegalStateException("no matrix have been decomposed yet");
}
}
}

25
src/java/org/apache/commons/math/linear/QRDecomposition.java
View File

@ -20,8 +20,13 @@ package org.apache.commons.math.linear;
/**
* 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 @@ public interface QRDecomposition extends DecompositionSolver {
* 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 @@ public interface QRDecomposition extends DecompositionSolver {
* 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;
}

165
src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
View File

@ -34,7 +34,7 @@ package org.apache.commons.math.linear;
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 @@ public class QRDecompositionImpl implements QRDecomposition {
* 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 @@ public class QRDecompositionImpl implements QRDecomposition {
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;
public QRDecompositionImpl() {
}
/**
* The column dimension of the given matrix. The size of R will be m x n.
*/
private final int n;
/**
* 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 @@ public class QRDecompositionImpl implements QRDecomposition {
}
/** {@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 @@ public class QRDecompositionImpl implements QRDecomposition {
}
/** {@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 @@ public class QRDecompositionImpl implements QRDecomposition {
}
/** {@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 @@ public class QRDecompositionImpl implements QRDecomposition {
}
/** {@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");
checkDecomposed();
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");
}
final double[] x = new double[n];
final double[] y = b.clone();
// apply Householder transforms to solve Q.y = b
for (int minor = 0; minor < Math.min(m, n); minor++) {
double dotProduct = 0;
for (int row = minor; row < m; row++) {
dotProduct += y[row] * qr[row][minor];
}
if (!isFullRank()) {
throw new InvalidMatrixException("Matrix is rank-deficient");
dotProduct /= rDiag[minor] * qr[minor][minor];
for (int row = minor; row < m; row++) {
y[row] += dotProduct * qr[row][minor];
}
final double[] x = new double[n];
final double[] y = b.clone();
// apply Householder transforms to solve Q.y = b
for (int minor = 0; minor < Math.min(m, n); 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 @@ public class QRDecompositionImpl implements QRDecomposition {
* <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 @@ public class QRDecompositionImpl implements QRDecomposition {
}
/**
* 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");
}
}
}

23
src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
View File

@ -88,6 +88,29 @@ public class LUDecompositionImplTest extends TestCase {
}
}
/** 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);

12
src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
View File

@ -349,4 +349,16 @@ public class QRDecompositionImplTest extends TestCase {
}
/** 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");
}
}
}