1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
|
package tim.prune.function.estimate.jama;
/**
* QR Decomposition.
*
* For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
* orthogonal matrix Q and an n-by-n upper triangular matrix R so that A = Q*R.
*
* The QR decomposition always exists, even if the matrix does not have full
* rank, so the constructor will never fail. The primary use of the QR
* decomposition is in the least squares solution of nonsquare systems of
* simultaneous linear equations. This will fail if isFullRank() returns false.
*
* Original authors The MathWorks, Inc. and the National Institute of Standards and Technology
* The original public domain code has now been modified and reduced,
* and is placed under GPL2 with the rest of the GpsPrune code.
*/
public class QRDecomposition
{
/** Array for internal storage of decomposition */
private double[][] _QR;
/** Row and column dimensions */
private int _m, _n;
/** Array for internal storage of diagonal of R */
private double[] _Rdiag;
/**
* QR Decomposition, computed by Householder reflections.
*
* @param inA Rectangular matrix
* @return Structure to access R and the Householder vectors and compute Q.
*/
public QRDecomposition(Matrix inA)
{
// Initialize.
_QR = inA.getArrayCopy();
_m = inA.getNumRows();
_n = inA.getNumColumns();
_Rdiag = new double[_n];
// Main loop.
for (int k = 0; k < _n; k++)
{
// Compute 2-norm of k-th column without under/overflow.
double nrm = 0;
for (int i = k; i < _m; i++) {
nrm = Maths.pythag(nrm, _QR[i][k]);
}
if (nrm != 0.0)
{
// Form k-th Householder vector.
if (_QR[k][k] < 0) {
nrm = -nrm;
}
for (int i = k; i < _m; i++) {
_QR[i][k] /= nrm;
}
_QR[k][k] += 1.0;
// Apply transformation to remaining columns.
for (int j = k + 1; j < _n; j++)
{
double s = 0.0;
for (int i = k; i < _m; i++) {
s += _QR[i][k] * _QR[i][j];
}
s = -s / _QR[k][k];
for (int i = k; i < _m; i++) {
_QR[i][j] += s * _QR[i][k];
}
}
}
_Rdiag[k] = -nrm;
}
}
/*
* ------------------------ Public Methods ------------------------
*/
/**
* Is the matrix full rank?
* @return true if R, and hence A, has full rank.
*/
public boolean isFullRank()
{
for (int j = 0; j < _n; j++) {
if (_Rdiag[j] == 0)
return false;
}
return true;
}
/**
* Return the Householder vectors
* @deprecated
* @return Lower trapezoidal matrix whose columns define the reflections
*/
private Matrix getH()
{
Matrix X = new Matrix(_m, _n);
double[][] H = X.getArray();
for (int i = 0; i < _m; i++) {
for (int j = 0; j < _n; j++) {
if (i >= j) {
H[i][j] = _QR[i][j];
} else {
H[i][j] = 0.0;
}
}
}
return X;
}
/**
* Return the upper triangular factor
* @deprecated
* @return R
*/
private Matrix getR()
{
Matrix X = new Matrix(_n, _n);
double[][] R = X.getArray();
for (int i = 0; i < _n; i++) {
for (int j = 0; j < _n; j++) {
if (i < j) {
R[i][j] = _QR[i][j];
} else if (i == j) {
R[i][j] = _Rdiag[i];
} else {
R[i][j] = 0.0;
}
}
}
return X;
}
/**
* Generate and return the (economy-sized) orthogonal factor
* @deprecated
* @return Q
*/
private Matrix getQ()
{
Matrix X = new Matrix(_m, _n);
double[][] Q = X.getArray();
for (int k = _n - 1; k >= 0; k--) {
for (int i = 0; i < _m; i++) {
Q[i][k] = 0.0;
}
Q[k][k] = 1.0;
for (int j = k; j < _n; j++) {
if (_QR[k][k] != 0) {
double s = 0.0;
for (int i = k; i < _m; i++) {
s += _QR[i][k] * Q[i][j];
}
s = -s / _QR[k][k];
for (int i = k; i < _m; i++) {
Q[i][j] += s * _QR[i][k];
}
}
}
}
return X;
}
/**
* Least squares solution of A*X = B
* @param B A Matrix with as many rows as A and any number of columns
* @return X that minimizes the two norm of Q*R*X-B
* @exception IllegalArgumentException if matrix dimensions don't agree
* @exception RuntimeException if Matrix is rank deficient.
*/
public Matrix solve(Matrix B)
{
if (B.getNumRows() != _m) {
throw new IllegalArgumentException("Matrix row dimensions must agree.");
}
if (!isFullRank()) {
throw new RuntimeException("Matrix is rank deficient.");
}
// Copy right hand side
int nx = B.getNumColumns();
double[][] X = B.getArrayCopy();
// Compute Y = transpose(Q)*B
for (int k = 0; k < _n; k++) {
for (int j = 0; j < nx; j++) {
double s = 0.0;
for (int i = k; i < _m; i++) {
s += _QR[i][k] * X[i][j];
}
s = -s / _QR[k][k];
for (int i = k; i < _m; i++) {
X[i][j] += s * _QR[i][k];
}
}
}
// Solve R*X = Y;
for (int k = _n - 1; k >= 0; k--) {
for (int j = 0; j < nx; j++) {
X[k][j] /= _Rdiag[k];
}
for (int i = 0; i < k; i++) {
for (int j = 0; j < nx; j++) {
X[i][j] -= X[k][j] * _QR[i][k];
}
}
}
return (new Matrix(X, _n, nx).getMatrix(0, _n - 1, 0, nx - 1));
}
}
|
