diff --git a/src/java/org/apache/commons/math/linear/QRDecomposition.java b/src/java/org/apache/commons/math/linear/QRDecomposition.java index b2adb68b5..f6daa18e9 100644 --- a/src/java/org/apache/commons/math/linear/QRDecomposition.java +++ b/src/java/org/apache/commons/math/linear/QRDecomposition.java @@ -57,6 +57,15 @@ public interface QRDecomposition extends DecompositionSolver { */ RealMatrix getQ() throws IllegalStateException; + /** + * Returns the transpose of the matrix Q of the decomposition. + *
Q is an orthogonal matrix
+ * @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. *H is a lower trapezoidal matrix whose columns represent diff --git a/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java b/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java index fa85e900b..af47babaf 100644 --- a/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java +++ b/src/java/org/apache/commons/math/linear/QRDecompositionImpl.java @@ -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; } diff --git a/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java b/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java index 2349bb1a7..0c05c423d 100644 --- a/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java +++ b/src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java @@ -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();