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 getQT() method to QRDecomposition interface 

compute Q by transposing QT rather than the other way round for efficiency
JIRA: MATH-223

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@699850 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-09-28 16:17:41 +00:00
parent 19a6674a57
commit 8260edaaf6
3 changed files with 49 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -57,6 +57,15 @@ public interface QRDecomposition extends DecompositionSolver {
*/
RealMatrix getQ() throws IllegalStateException;
/**
* Returns the transpose of 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 getQT() throws IllegalStateException;
/**
* Returns the Householder reflector vectors.
* <p>H is a lower trapezoidal matrix whose columns represent

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

@ -52,6 +52,9 @@ public class QRDecompositionImpl implements QRDecomposition {
/** Cached value of Q. */
private RealMatrix cachedQ;
/** Cached value of QT. */
private RealMatrix cachedQT;
/** Cached value of R. */
private RealMatrix cachedR;
@ -87,9 +90,10 @@ public class QRDecompositionImpl implements QRDecomposition {
final int n = matrix.getColumnDimension();
qrt = matrix.transpose().getData();
rDiag = new double[n];
cachedQ = null;
cachedR = null;
cachedH = null;
cachedQ = null;
cachedQT = null;
cachedR = null;
cachedH = null;
/*
* The QR decomposition of a matrix A is calculated using Householder
@ -191,15 +195,24 @@ public class QRDecompositionImpl implements QRDecomposition {
/** {@inheritDoc} */
public RealMatrix getQ()
throws IllegalStateException {
if (cachedQ == null) {
cachedQ = getQT().transpose();
}
return cachedQ;
}
/** {@inheritDoc} */
public RealMatrix getQT()
throws IllegalStateException {
if (cachedQ == null) {
checkDecomposed();
// Q is supposed to be m x m
// QT is supposed to be m x m
final int n = qrt.length;
final int m = qrt[0].length;
double[][] q = new double[m][m];
double[][] qT = new double[m][m];
/*
* Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
@ -207,34 +220,35 @@ public class QRDecompositionImpl implements QRDecomposition {
* succession to the result
*/
for (int minor = m - 1; minor >= Math.min(m, n); minor--) {
q[minor][minor]=1;
qT[minor][minor]=1;
}
for (int minor = Math.min(m,n)-1; minor >= 0; minor--){
final double[] qrtMinor = qrt[minor];
q[minor][minor] = 1;
qT[minor][minor] = 1;
if (qrtMinor[minor] != 0.0) {
for (int col = minor; col < m; col++) {
final double[] qTCol = qT[col];
double alpha = 0;
for (int row = minor; row < m; row++) {
alpha -= q[row][col] * qrtMinor[row];
alpha -= qTCol[row] * qrtMinor[row];
}
alpha /= rDiag[minor] * qrtMinor[minor];
for (int row = minor; row < m; row++) {
q[row][col] -= alpha * qrtMinor[row];
qTCol[row] -= alpha * qrtMinor[row];
}
}
}
}
// cache the matrix for subsequent calls
cachedQ = new RealMatrixImpl(q, false);
cachedQT = new RealMatrixImpl(qT, false);
}
// return the cached matrix
return cachedQ;
return cachedQT;
}

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

@ -109,31 +109,31 @@ public class QRDecompositionImplTest extends TestCase {
/** test the orthogonality of Q */
public void testQOrthogonal() {
RealMatrix matrix = new RealMatrixImpl(testData3x3NonSingular, false);
matrix = new QRDecompositionImpl(matrix).getQ();
RealMatrix q = new QRDecompositionImpl(matrix).getQ();
RealMatrix qT = new QRDecompositionImpl(matrix).getQT();
RealMatrix eye = MatrixUtils.createRealIdentityMatrix(3);
double norm = matrix.transpose().multiply(matrix).subtract(eye)
.getNorm();
double norm = qT.multiply(q).subtract(eye).getNorm();
assertEquals("3x3 nonsingular Q'Q = I", 0, norm, normTolerance);
matrix = new RealMatrixImpl(testData3x3Singular, false);
matrix = new QRDecompositionImpl(matrix).getQ();
q = new QRDecompositionImpl(matrix).getQ();
qT = new QRDecompositionImpl(matrix).getQT();
eye = MatrixUtils.createRealIdentityMatrix(3);
norm = matrix.transpose().multiply(matrix).subtract(eye)
.getNorm();
norm = qT.multiply(q).subtract(eye).getNorm();
assertEquals("3x3 singular Q'Q = I", 0, norm, normTolerance);
matrix = new RealMatrixImpl(testData3x4, false);
matrix = new QRDecompositionImpl(matrix).getQ();
q = new QRDecompositionImpl(matrix).getQ();
qT = new QRDecompositionImpl(matrix).getQT();
eye = MatrixUtils.createRealIdentityMatrix(3);
norm = matrix.transpose().multiply(matrix).subtract(eye)
.getNorm();
norm = qT.multiply(q).subtract(eye).getNorm();
assertEquals("3x4 Q'Q = I", 0, norm, normTolerance);
matrix = new RealMatrixImpl(testData4x3, false);
matrix = new QRDecompositionImpl(matrix).getQ();
q = new QRDecompositionImpl(matrix).getQ();
qT = new QRDecompositionImpl(matrix).getQT();
eye = MatrixUtils.createRealIdentityMatrix(4);
norm = matrix.transpose().multiply(matrix).subtract(eye)
.getNorm();
norm = qT.multiply(q).subtract(eye).getNorm();
assertEquals("4x3 Q'Q = I", 0, norm, normTolerance);
}
@ -337,6 +337,8 @@ public class QRDecompositionImplTest extends TestCase {
// check values against known references
RealMatrix q = qr.getQ();
assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
RealMatrix qT = qr.getQT();
assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13);
RealMatrix r = qr.getR();
assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
RealMatrix h = qr.getH();