001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import java.util.Arrays;
021    
022    import org.apache.commons.math.MathRuntimeException;
023    import org.apache.commons.math.exception.util.LocalizedFormats;
024    import org.apache.commons.math.util.FastMath;
025    
026    
027    /**
028     * Calculates the QR-decomposition of a matrix.
029     * <p>The QR-decomposition of a matrix A consists of two matrices Q and R
030     * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
031     * upper triangular. If A is m&times;n, Q is m&times;m and R m&times;n.</p>
032     * <p>This class compute the decomposition using Householder reflectors.</p>
033     * <p>For efficiency purposes, the decomposition in packed form is transposed.
034     * This allows inner loop to iterate inside rows, which is much more cache-efficient
035     * in Java.</p>
036     *
037     * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
038     * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
039     *
040     * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
041     * @since 1.2
042     */
043    public class QRDecompositionImpl implements QRDecomposition {
044    
045        /**
046         * A packed TRANSPOSED representation of the QR decomposition.
047         * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
048         * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
049         * from which an explicit form of Q can be recomputed if desired.</p>
050         */
051        private double[][] qrt;
052    
053        /** The diagonal elements of R. */
054        private double[] rDiag;
055    
056        /** Cached value of Q. */
057        private RealMatrix cachedQ;
058    
059        /** Cached value of QT. */
060        private RealMatrix cachedQT;
061    
062        /** Cached value of R. */
063        private RealMatrix cachedR;
064    
065        /** Cached value of H. */
066        private RealMatrix cachedH;
067    
068        /**
069         * Calculates the QR-decomposition of the given matrix.
070         * @param matrix The matrix to decompose.
071         */
072        public QRDecompositionImpl(RealMatrix matrix) {
073    
074            final int m = matrix.getRowDimension();
075            final int n = matrix.getColumnDimension();
076            qrt = matrix.transpose().getData();
077            rDiag = new double[FastMath.min(m, n)];
078            cachedQ  = null;
079            cachedQT = null;
080            cachedR  = null;
081            cachedH  = null;
082    
083            /*
084             * The QR decomposition of a matrix A is calculated using Householder
085             * reflectors by repeating the following operations to each minor
086             * A(minor,minor) of A:
087             */
088            for (int minor = 0; minor < FastMath.min(m, n); minor++) {
089    
090                final double[] qrtMinor = qrt[minor];
091    
092                /*
093                 * Let x be the first column of the minor, and a^2 = |x|^2.
094                 * x will be in the positions qr[minor][minor] through qr[m][minor].
095                 * The first column of the transformed minor will be (a,0,0,..)'
096                 * The sign of a is chosen to be opposite to the sign of the first
097                 * component of x. Let's find a:
098                 */
099                double xNormSqr = 0;
100                for (int row = minor; row < m; row++) {
101                    final double c = qrtMinor[row];
102                    xNormSqr += c * c;
103                }
104                final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
105                rDiag[minor] = a;
106    
107                if (a != 0.0) {
108    
109                    /*
110                     * Calculate the normalized reflection vector v and transform
111                     * the first column. We know the norm of v beforehand: v = x-ae
112                     * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
113                     * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
114                     * Here <x, e> is now qr[minor][minor].
115                     * v = x-ae is stored in the column at qr:
116                     */
117                    qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
118    
119                    /*
120                     * Transform the rest of the columns of the minor:
121                     * They will be transformed by the matrix H = I-2vv'/|v|^2.
122                     * If x is a column vector of the minor, then
123                     * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
124                     * Therefore the transformation is easily calculated by
125                     * subtracting the column vector (2<x,v>/|v|^2)v from x.
126                     *
127                     * Let 2<x,v>/|v|^2 = alpha. From above we have
128                     * |v|^2 = -2a*(qr[minor][minor]), so
129                     * alpha = -<x,v>/(a*qr[minor][minor])
130                     */
131                    for (int col = minor+1; col < n; col++) {
132                        final double[] qrtCol = qrt[col];
133                        double alpha = 0;
134                        for (int row = minor; row < m; row++) {
135                            alpha -= qrtCol[row] * qrtMinor[row];
136                        }
137                        alpha /= a * qrtMinor[minor];
138    
139                        // Subtract the column vector alpha*v from x.
140                        for (int row = minor; row < m; row++) {
141                            qrtCol[row] -= alpha * qrtMinor[row];
142                        }
143                    }
144                }
145            }
146        }
147    
148        /** {@inheritDoc} */
149        public RealMatrix getR() {
150    
151            if (cachedR == null) {
152    
153                // R is supposed to be m x n
154                final int n = qrt.length;
155                final int m = qrt[0].length;
156                cachedR = MatrixUtils.createRealMatrix(m, n);
157    
158                // copy the diagonal from rDiag and the upper triangle of qr
159                for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
160                    cachedR.setEntry(row, row, rDiag[row]);
161                    for (int col = row + 1; col < n; col++) {
162                        cachedR.setEntry(row, col, qrt[col][row]);
163                    }
164                }
165    
166            }
167    
168            // return the cached matrix
169            return cachedR;
170    
171        }
172    
173        /** {@inheritDoc} */
174        public RealMatrix getQ() {
175            if (cachedQ == null) {
176                cachedQ = getQT().transpose();
177            }
178            return cachedQ;
179        }
180    
181        /** {@inheritDoc} */
182        public RealMatrix getQT() {
183    
184            if (cachedQT == null) {
185    
186                // QT is supposed to be m x m
187                final int n = qrt.length;
188                final int m = qrt[0].length;
189                cachedQT = MatrixUtils.createRealMatrix(m, m);
190    
191                /*
192                 * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
193                 * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
194                 * succession to the result
195                 */
196                for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
197                    cachedQT.setEntry(minor, minor, 1.0);
198                }
199    
200                for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
201                    final double[] qrtMinor = qrt[minor];
202                    cachedQT.setEntry(minor, minor, 1.0);
203                    if (qrtMinor[minor] != 0.0) {
204                        for (int col = minor; col < m; col++) {
205                            double alpha = 0;
206                            for (int row = minor; row < m; row++) {
207                                alpha -= cachedQT.getEntry(col, row) * qrtMinor[row];
208                            }
209                            alpha /= rDiag[minor] * qrtMinor[minor];
210    
211                            for (int row = minor; row < m; row++) {
212                                cachedQT.addToEntry(col, row, -alpha * qrtMinor[row]);
213                            }
214                        }
215                    }
216                }
217    
218            }
219    
220            // return the cached matrix
221            return cachedQT;
222    
223        }
224    
225        /** {@inheritDoc} */
226        public RealMatrix getH() {
227    
228            if (cachedH == null) {
229    
230                final int n = qrt.length;
231                final int m = qrt[0].length;
232                cachedH = MatrixUtils.createRealMatrix(m, n);
233                for (int i = 0; i < m; ++i) {
234                    for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
235                        cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]);
236                    }
237                }
238    
239            }
240    
241            // return the cached matrix
242            return cachedH;
243    
244        }
245    
246        /** {@inheritDoc} */
247        public DecompositionSolver getSolver() {
248            return new Solver(qrt, rDiag);
249        }
250    
251        /** Specialized solver. */
252        private static class Solver implements DecompositionSolver {
253    
254            /**
255             * A packed TRANSPOSED representation of the QR decomposition.
256             * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
257             * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
258             * from which an explicit form of Q can be recomputed if desired.</p>
259             */
260            private final double[][] qrt;
261    
262            /** The diagonal elements of R. */
263            private final double[] rDiag;
264    
265            /**
266             * Build a solver from decomposed matrix.
267             * @param qrt packed TRANSPOSED representation of the QR decomposition
268             * @param rDiag diagonal elements of R
269             */
270            private Solver(final double[][] qrt, final double[] rDiag) {
271                this.qrt   = qrt;
272                this.rDiag = rDiag;
273            }
274    
275            /** {@inheritDoc} */
276            public boolean isNonSingular() {
277    
278                for (double diag : rDiag) {
279                    if (diag == 0) {
280                        return false;
281                    }
282                }
283                return true;
284    
285            }
286    
287            /** {@inheritDoc} */
288            public double[] solve(double[] b)
289            throws IllegalArgumentException, InvalidMatrixException {
290    
291                final int n = qrt.length;
292                final int m = qrt[0].length;
293                if (b.length != m) {
294                    throw MathRuntimeException.createIllegalArgumentException(
295                            LocalizedFormats.VECTOR_LENGTH_MISMATCH,
296                            b.length, m);
297                }
298                if (!isNonSingular()) {
299                    throw new SingularMatrixException();
300                }
301    
302                final double[] x = new double[n];
303                final double[] y = b.clone();
304    
305                // apply Householder transforms to solve Q.y = b
306                for (int minor = 0; minor < FastMath.min(m, n); minor++) {
307    
308                    final double[] qrtMinor = qrt[minor];
309                    double dotProduct = 0;
310                    for (int row = minor; row < m; row++) {
311                        dotProduct += y[row] * qrtMinor[row];
312                    }
313                    dotProduct /= rDiag[minor] * qrtMinor[minor];
314    
315                    for (int row = minor; row < m; row++) {
316                        y[row] += dotProduct * qrtMinor[row];
317                    }
318    
319                }
320    
321                // solve triangular system R.x = y
322                for (int row = rDiag.length - 1; row >= 0; --row) {
323                    y[row] /= rDiag[row];
324                    final double yRow   = y[row];
325                    final double[] qrtRow = qrt[row];
326                    x[row] = yRow;
327                    for (int i = 0; i < row; i++) {
328                        y[i] -= yRow * qrtRow[i];
329                    }
330                }
331    
332                return x;
333    
334            }
335    
336            /** {@inheritDoc} */
337            public RealVector solve(RealVector b)
338            throws IllegalArgumentException, InvalidMatrixException {
339                try {
340                    return solve((ArrayRealVector) b);
341                } catch (ClassCastException cce) {
342                    return new ArrayRealVector(solve(b.getData()), false);
343                }
344            }
345    
346            /** Solve the linear equation A &times; X = B.
347             * <p>The A matrix is implicit here. It is </p>
348             * @param b right-hand side of the equation A &times; X = B
349             * @return a vector X that minimizes the two norm of A &times; X - B
350             * @throws IllegalArgumentException if matrices dimensions don't match
351             * @throws InvalidMatrixException if decomposed matrix is singular
352             */
353            public ArrayRealVector solve(ArrayRealVector b)
354            throws IllegalArgumentException, InvalidMatrixException {
355                return new ArrayRealVector(solve(b.getDataRef()), false);
356            }
357    
358            /** {@inheritDoc} */
359            public RealMatrix solve(RealMatrix b)
360            throws IllegalArgumentException, InvalidMatrixException {
361    
362                final int n = qrt.length;
363                final int m = qrt[0].length;
364                if (b.getRowDimension() != m) {
365                    throw MathRuntimeException.createIllegalArgumentException(
366                            LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
367                            b.getRowDimension(), b.getColumnDimension(), m, "n");
368                }
369                if (!isNonSingular()) {
370                    throw new SingularMatrixException();
371                }
372    
373                final int columns        = b.getColumnDimension();
374                final int blockSize      = BlockRealMatrix.BLOCK_SIZE;
375                final int cBlocks        = (columns + blockSize - 1) / blockSize;
376                final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
377                final double[][] y       = new double[b.getRowDimension()][blockSize];
378                final double[]   alpha   = new double[blockSize];
379    
380                for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
381                    final int kStart = kBlock * blockSize;
382                    final int kEnd   = FastMath.min(kStart + blockSize, columns);
383                    final int kWidth = kEnd - kStart;
384    
385                    // get the right hand side vector
386                    b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
387    
388                    // apply Householder transforms to solve Q.y = b
389                    for (int minor = 0; minor < FastMath.min(m, n); minor++) {
390                        final double[] qrtMinor = qrt[minor];
391                        final double factor     = 1.0 / (rDiag[minor] * qrtMinor[minor]);
392    
393                        Arrays.fill(alpha, 0, kWidth, 0.0);
394                        for (int row = minor; row < m; ++row) {
395                            final double   d    = qrtMinor[row];
396                            final double[] yRow = y[row];
397                            for (int k = 0; k < kWidth; ++k) {
398                                alpha[k] += d * yRow[k];
399                            }
400                        }
401                        for (int k = 0; k < kWidth; ++k) {
402                            alpha[k] *= factor;
403                        }
404    
405                        for (int row = minor; row < m; ++row) {
406                            final double   d    = qrtMinor[row];
407                            final double[] yRow = y[row];
408                            for (int k = 0; k < kWidth; ++k) {
409                                yRow[k] += alpha[k] * d;
410                            }
411                        }
412    
413                    }
414    
415                    // solve triangular system R.x = y
416                    for (int j = rDiag.length - 1; j >= 0; --j) {
417                        final int      jBlock = j / blockSize;
418                        final int      jStart = jBlock * blockSize;
419                        final double   factor = 1.0 / rDiag[j];
420                        final double[] yJ     = y[j];
421                        final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
422                        int index = (j - jStart) * kWidth;
423                        for (int k = 0; k < kWidth; ++k) {
424                            yJ[k]          *= factor;
425                            xBlock[index++] = yJ[k];
426                        }
427    
428                        final double[] qrtJ = qrt[j];
429                        for (int i = 0; i < j; ++i) {
430                            final double rIJ  = qrtJ[i];
431                            final double[] yI = y[i];
432                            for (int k = 0; k < kWidth; ++k) {
433                                yI[k] -= yJ[k] * rIJ;
434                            }
435                        }
436    
437                    }
438    
439                }
440    
441                return new BlockRealMatrix(n, columns, xBlocks, false);
442    
443            }
444    
445            /** {@inheritDoc} */
446            public RealMatrix getInverse()
447            throws InvalidMatrixException {
448                return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
449            }
450    
451        }
452    
453    }