2
* Licensed to the Apache Software Foundation (ASF) under one or more
3
* contributor license agreements. See the NOTICE file distributed with
4
* this work for additional information regarding copyright ownership.
5
* The ASF licenses this file to You under the Apache License, Version 2.0
6
* (the "License"); you may not use this file except in compliance with
7
* the License. You may obtain a copy of the License at
9
* http://www.apache.org/licenses/LICENSE-2.0
11
* Unless required by applicable law or agreed to in writing, software
12
* distributed under the License is distributed on an "AS IS" BASIS,
13
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14
* See the License for the specific language governing permissions and
15
* limitations under the License.
18
package org.apache.commons.math.linear;
21
* Calculates the QR-decomposition of a matrix. In the QR-decomposition of
22
* a matrix A consists of two matrices Q and R that satisfy: A = QR, Q is
23
* orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular. If A is
24
* m×n, Q is m×m and R m×n.
26
* Implemented using Householder reflectors.</p>
28
* @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
29
* @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
31
* @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
34
public class QRDecompositionImpl implements QRDecomposition {
37
* A packed representation of the QR decomposition. The elements above the
38
* diagonal are the elements of R, and the columns of the lower triangle
39
* are the Householder reflector vectors of which an explicit form of Q can
42
private double[][] qr;
45
* The diagonal elements of R.
47
private double[] rDiag;
50
* The row dimension of the given matrix. The size of Q will be m x m, the
51
* size of R will be m x n.
56
* The column dimension of the given matrix. The size of R will be m x n.
61
* Calculates the QR decomposition of the given matrix.
63
* @param matrix The matrix to decompose.
65
public QRDecompositionImpl(RealMatrix matrix) {
66
m = matrix.getRowDimension();
67
n = matrix.getColumnDimension();
68
qr = matrix.getData();
69
rDiag = new double[n];
72
* The QR decomposition of a matrix A is calculated using Householder
73
* reflectors by repeating the following operations to each minor
74
* A(minor,minor) of A:
76
for (int minor = 0; minor < Math.min(m, n); minor++) {
78
* Let x be the first column of the minor, and a^2 = |x|^2.
79
* x will be in the positions qr[minor][minor] through qr[m][minor].
80
* The first column of the transformed minor will be (a,0,0,..)'
81
* The sign of a is chosen to be opposite to the sign of the first
82
* component of x. Let's find a:
85
for (int row = minor; row < m; row++) {
86
xNormSqr += qr[row][minor]*qr[row][minor];
88
double a = Math.sqrt(xNormSqr);
89
if (qr[minor][minor] > 0) a = -a;
95
* Calculate the normalized reflection vector v and transform
96
* the first column. We know the norm of v beforehand: v = x-ae
97
* so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
98
* a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
99
* Here <x, e> is now qr[minor][minor].
100
* v = x-ae is stored in the column at qr:
102
qr[minor][minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
105
* Transform the rest of the columns of the minor:
106
* They will be transformed by the matrix H = I-2vv'/|v|^2.
107
* If x is a column vector of the minor, then
108
* Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
109
* Therefore the transformation is easily calculated by
110
* subtracting the column vector (2<x,v>/|v|^2)v from x.
112
* Let 2<x,v>/|v|^2 = alpha. From above we have
113
* |v|^2 = -2a*(qr[minor][minor]), so
114
* alpha = -<x,v>/(a*qr[minor][minor])
116
for (int col = minor+1; col < n; col++) {
118
for (int row = minor; row < m; row++) {
119
alpha -= qr[row][col]*qr[row][minor];
121
alpha /= a*qr[minor][minor];
123
// Subtract the column vector alpha*v from x.
124
for (int row = minor; row < m; row++) {
125
qr[row][col] -= alpha*qr[row][minor];
133
* Returns the matrix R of the QR-decomposition.
135
* @return the R matrix
137
public RealMatrix getR()
139
// R is supposed to be m x n
140
RealMatrixImpl ret = new RealMatrixImpl(m,n);
141
double[][] r = ret.getDataRef();
143
// copy the diagonal from rDiag and the upper triangle of qr
144
for (int row = Math.min(m,n)-1; row >= 0; row--) {
145
r[row][row] = rDiag[row];
146
for (int col = row+1; col < n; col++) {
147
r[row][col] = qr[row][col];
154
* Returns the matrix Q of the QR-decomposition.
156
* @return the Q matrix
158
public RealMatrix getQ()
160
// Q is supposed to be m x m
161
RealMatrixImpl ret = new RealMatrixImpl(m,m);
162
double[][] Q = ret.getDataRef();
165
* Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
166
* applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
167
* succession to the result
169
for (int minor = m-1; minor >= Math.min(m,n); minor--) {
173
for (int minor = Math.min(m,n)-1; minor >= 0; minor--){
175
if (qr[minor][minor] != 0.0) {
176
for (int col = minor; col < m; col++) {
178
for (int row = minor; row < m; row++) {
179
alpha -= Q[row][col] * qr[row][minor];
181
alpha /= rDiag[minor]*qr[minor][minor];
183
for (int row = minor; row < m; row++) {
184
Q[row][col] -= alpha*qr[row][minor];
added
removed
src/java/org/apache/commons/math/linear/QRDecompositionImpl.java
