Essentia  2.1-beta5-dev
jama_qr.h
Go to the documentation of this file.
1 #ifndef JAMA_QR_H
2 #define JAMA_QR_H
3 
4 #include "tnt_array1d.h"
5 #include "tnt_array2d.h"
6 #include "tnt_math_utils.h"
7 
8 namespace JAMA
9 {
10 
34 template <class Real>
35 class QR {
36 
37 
43 
48  int m, n;
49 
54 
55 
56 public:
57 
63  QR(const TNT::Array2D<Real> &A) /* constructor */
64  {
65  QR_ = A.copy();
66  m = A.dim1();
67  n = A.dim2();
68  Rdiag = TNT::Array1D<Real>(n);
69  int i=0, j=0, k=0;
70 
71  // Main loop.
72  for (k = 0; k < n; k++) {
73  // Compute 2-norm of k-th column without under/overflow.
74  Real nrm = 0;
75  for (i = k; i < m; i++) {
76  nrm = TNT::hypot(nrm,QR_[i][k]);
77  }
78 
79  if (nrm != 0.0) {
80  // Form k-th Householder vector.
81  if (QR_[k][k] < 0) {
82  nrm = -nrm;
83  }
84  for (i = k; i < m; i++) {
85  QR_[i][k] /= nrm;
86  }
87  QR_[k][k] += 1.0;
88 
89  // Apply transformation to remaining columns.
90  for (j = k+1; j < n; j++) {
91  Real s = 0.0;
92  for (i = k; i < m; i++) {
93  s += QR_[i][k]*QR_[i][j];
94  }
95  s = -s/QR_[k][k];
96  for (i = k; i < m; i++) {
97  QR_[i][j] += s*QR_[i][k];
98  }
99  }
100  }
101  Rdiag[k] = -nrm;
102  }
103  }
104 
105 
111  int isFullRank() const
112  {
113  for (int j = 0; j < n; j++)
114  {
115  if (Rdiag[j] == 0)
116  return 0;
117  }
118  return 1;
119  }
120 
121 
122 
123 
131  {
132  TNT::Array2D<Real> H(m,n);
133 
134  /* note: H is completely filled in by algorithm, so
135  initializaiton of H is not necessary.
136  */
137  for (int i = 0; i < m; i++)
138  {
139  for (int j = 0; j < n; j++)
140  {
141  if (i >= j) {
142  H[i][j] = QR_[i][j];
143  } else {
144  H[i][j] = 0.0;
145  }
146  }
147  }
148  return H;
149  }
150 
151 
152 
158  {
159  TNT::Array2D<Real> R(n,n);
160  for (int i = 0; i < n; i++) {
161  for (int j = 0; j < n; j++) {
162  if (i < j) {
163  R[i][j] = QR_[i][j];
164  } else if (i == j) {
165  R[i][j] = Rdiag[i];
166  } else {
167  R[i][j] = 0.0;
168  }
169  }
170  }
171  return R;
172  }
173 
174 
175 
176 
177 
184  {
185  int i=0, j=0, k=0;
186 
187  TNT::Array2D<Real> Q(m,n);
188  for (k = n-1; k >= 0; k--) {
189  for (i = 0; i < m; i++) {
190  Q[i][k] = 0.0;
191  }
192  Q[k][k] = 1.0;
193  for (j = k; j < n; j++) {
194  if (QR_[k][k] != 0) {
195  Real s = 0.0;
196  for (i = k; i < m; i++) {
197  s += QR_[i][k]*Q[i][j];
198  }
199  s = -s/QR_[k][k];
200  for (i = k; i < m; i++) {
201  Q[i][j] += s*QR_[i][k];
202  }
203  }
204  }
205  }
206  return Q;
207  }
208 
209 
218  {
219  if (b.dim1() != m) /* arrays must be conformant */
220  return TNT::Array1D<Real>();
221 
222  if ( !isFullRank() ) /* matrix is rank deficient */
223  {
224  return TNT::Array1D<Real>();
225  }
226 
227  TNT::Array1D<Real> x = b.copy();
228 
229  // Compute Y = transpose(Q)*b
230  for (int k = 0; k < n; k++)
231  {
232  Real s = 0.0;
233  for (int i = k; i < m; i++)
234  {
235  s += QR_[i][k]*x[i];
236  }
237  s = -s/QR_[k][k];
238  for (int i = k; i < m; i++)
239  {
240  x[i] += s*QR_[i][k];
241  }
242  }
243  // Solve R*X = Y;
244  for (int k = n-1; k >= 0; k--)
245  {
246  x[k] /= Rdiag[k];
247  for (int i = 0; i < k; i++) {
248  x[i] -= x[k]*QR_[i][k];
249  }
250  }
251 
252 
253  /* return n x nx portion of X */
254  TNT::Array1D<Real> x_(n);
255  for (int i=0; i<n; i++)
256  x_[i] = x[i];
257 
258  return x_;
259  }
260 
269  {
270  if (B.dim1() != m) /* arrays must be conformant */
271  return TNT::Array2D<Real>(0,0);
272 
273  if ( !isFullRank() ) /* matrix is rank deficient */
274  {
275  return TNT::Array2D<Real>(0,0);
276  }
277 
278  int nx = B.dim2();
279  TNT::Array2D<Real> X = B.copy();
280  int i=0, j=0, k=0;
281 
282  // Compute Y = transpose(Q)*B
283  for (k = 0; k < n; k++) {
284  for (j = 0; j < nx; j++) {
285  Real s = 0.0;
286  for (i = k; i < m; i++) {
287  s += QR_[i][k]*X[i][j];
288  }
289  s = -s/QR_[k][k];
290  for (i = k; i < m; i++) {
291  X[i][j] += s*QR_[i][k];
292  }
293  }
294  }
295  // Solve R*X = Y;
296  for (k = n-1; k >= 0; k--) {
297  for (j = 0; j < nx; j++) {
298  X[k][j] /= Rdiag[k];
299  }
300  for (i = 0; i < k; i++) {
301  for (j = 0; j < nx; j++) {
302  X[i][j] -= X[k][j]*QR_[i][k];
303  }
304  }
305  }
306 
307 
308  /* return n x nx portion of X */
309  TNT::Array2D<Real> X_(n,nx);
310  for (i=0; i<n; i++)
311  for (j=0; j<nx; j++)
312  X_[i][j] = X[i][j];
313 
314  return X_;
315  }
316 
317 
318 };
319 
320 
321 }
322 // namespace JAMA
323 
324 #endif
325 // JAMA_QR__H
326 
TNT::Array2D< Real > getHouseholder(void) const
Definition: jama_qr.h:130
int dim1() const
Definition: tnt_array2d.h:231
int n
Definition: jama_qr.h:48
Array2D copy() const
Definition: tnt_array2d.h:180
Definition: jama_cholesky.h:8
Definition: jama_qr.h:35
Real hypot(const Real &a, const Real &b)
Definition: tnt_math_utils.h:18
int dim2() const
Definition: tnt_array2d.h:234
int dim1() const
Definition: tnt_array1d.h:218
float Real
Definition: types.h:68
int m
Definition: jama_qr.h:48
TNT::Array2D< Real > solve(const TNT::Array2D< Real > &B) const
Definition: jama_qr.h:268
TNT::Array2D< Real > QR_
Definition: jama_qr.h:42
TNT::Array1D< Real > solve(const TNT::Array1D< Real > &b) const
Definition: jama_qr.h:217
TNT::Array2D< Real > getR() const
Definition: jama_qr.h:157
int isFullRank() const
Definition: jama_qr.h:111
Array1D copy() const
Definition: tnt_array1d.h:176
QR(const TNT::Array2D< Real > &A)
Definition: jama_qr.h:63
TNT::Array2D< Real > getQ() const
Definition: jama_qr.h:183
TNT::Array1D< Real > Rdiag
Definition: jama_qr.h:53