Author: luc
Date: Sun Sep 28 09:17:41 2008
New Revision: 699850
URL: http://svn.apache.org/viewvc?rev=699850&view=rev
Log:
added a getQT() method to QRDecomposition interface
compute Q by transposing QT rather than the other way round for efficiency
JIRA: MATH-223
Modified:
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/QRDecompositionImplTest.java
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=699850&r1=699849&r2=699850&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 Sep 28 09:17:41 2008
@@ -58,6 +58,15 @@
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 [EMAIL PROTECTED]
+ * 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
* each successive Householder reflector vector. This matrix is used
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=699850&r1=699849&r2=699850&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 Sep 28 09:17:41 2008
@@ -52,6 +52,9 @@
/** 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 @@
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 @@
/** [EMAIL PROTECTED] */
public RealMatrix getQ()
throws IllegalStateException {
-
if (cachedQ == null) {
+ cachedQ = getQT().transpose();
+ }
+ return cachedQ;
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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 @@
* 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;
}
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=699850&r1=699849&r2=699850&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 Sep 28 09:17:41 2008
@@ -109,31 +109,31 @@
/** 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 @@
// 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();
