1 /* linalg/qr.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include <config.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <gsl/gsl_math.h>
26 #include <gsl/gsl_vector.h>
27 #include <gsl/gsl_matrix.h>
28 #include <gsl/gsl_blas.h>
29
30 #include <gsl/gsl_linalg.h>
31
32 #include "givens.c"
33 #include "apply_givens.c"
34
35 /* Factorise a general M x N matrix A into
36 *
37 * A = Q R
38 *
39 * where Q is orthogonal (M x M) and R is upper triangular (M x N).
40 *
41 * Q is stored as a packed set of Householder transformations in the
42 * strict lower triangular part of the input matrix.
43 *
44 * R is stored in the diagonal and upper triangle of the input matrix.
45 *
46 * The full matrix for Q can be obtained as the product
47 *
48 * Q = Q_k .. Q_2 Q_1
49 *
50 * where k = MIN(M,N) and
51 *
52 * Q_i = (I - tau_i * v_i * v_i')
53 *
54 * and where v_i is a Householder vector
55 *
56 * v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
57 *
58 * This storage scheme is the same as in LAPACK. */
59
60 int
gsl_linalg_QR_decomp(gsl_matrix * A,gsl_vector * tau)61 gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
62 {
63 const size_t M = A->size1;
64 const size_t N = A->size2;
65
66 if (tau->size != GSL_MIN (M, N))
67 {
68 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
69 }
70 else
71 {
72 size_t i;
73
74 for (i = 0; i < GSL_MIN (M, N); i++)
75 {
76 /* Compute the Householder transformation to reduce the j-th
77 column of the matrix to a multiple of the j-th unit vector */
78
79 gsl_vector_view c_full = gsl_matrix_column (A, i);
80 gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
81
82 double tau_i = gsl_linalg_householder_transform (&(c.vector));
83
84 gsl_vector_set (tau, i, tau_i);
85
86 /* Apply the transformation to the remaining columns and
87 update the norms */
88
89 if (i + 1 < N)
90 {
91 gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
92 gsl_linalg_householder_hm (tau_i, &(c.vector), &(m.matrix));
93 }
94 }
95
96 return GSL_SUCCESS;
97 }
98 }
99
100 /* Solves the system A x = b using the QR factorisation,
101
102 * R x = Q^T b
103 *
104 * to obtain x. Based on SLATEC code.
105 */
106
107 int
gsl_linalg_QR_solve(const gsl_matrix * QR,const gsl_vector * tau,const gsl_vector * b,gsl_vector * x)108 gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
109 {
110 if (QR->size1 != QR->size2)
111 {
112 GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
113 }
114 else if (QR->size1 != b->size)
115 {
116 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
117 }
118 else if (QR->size2 != x->size)
119 {
120 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
121 }
122 else
123 {
124 /* Copy x <- b */
125
126 gsl_vector_memcpy (x, b);
127
128 /* Solve for x */
129
130 gsl_linalg_QR_svx (QR, tau, x);
131
132 return GSL_SUCCESS;
133 }
134 }
135
136 /* Solves the system A x = b in place using the QR factorisation,
137
138 * R x = Q^T b
139 *
140 * to obtain x. Based on SLATEC code.
141 */
142
143 int
gsl_linalg_QR_svx(const gsl_matrix * QR,const gsl_vector * tau,gsl_vector * x)144 gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
145 {
146
147 if (QR->size1 != QR->size2)
148 {
149 GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
150 }
151 else if (QR->size1 != x->size)
152 {
153 GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
154 }
155 else
156 {
157 /* compute rhs = Q^T b */
158
159 gsl_linalg_QR_QTvec (QR, tau, x);
160
161 /* Solve R x = rhs, storing x in-place */
162
163 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
164
165 return GSL_SUCCESS;
166 }
167 }
168
169
170 /* Find the least squares solution to the overdetermined system
171 *
172 * A x = b
173 *
174 * for M >= N using the QR factorization A = Q R.
175 */
176
177 int
gsl_linalg_QR_lssolve(const gsl_matrix * QR,const gsl_vector * tau,const gsl_vector * b,gsl_vector * x,gsl_vector * residual)178 gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
179 {
180 const size_t M = QR->size1;
181 const size_t N = QR->size2;
182
183 if (M < N)
184 {
185 GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
186 }
187 else if (M != b->size)
188 {
189 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
190 }
191 else if (N != x->size)
192 {
193 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
194 }
195 else if (M != residual->size)
196 {
197 GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
198 }
199 else
200 {
201 gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
202 gsl_vector_view c = gsl_vector_subvector(residual, 0, N);
203
204 gsl_vector_memcpy(residual, b);
205
206 /* compute rhs = Q^T b */
207
208 gsl_linalg_QR_QTvec (QR, tau, residual);
209
210 /* Solve R x = rhs */
211
212 gsl_vector_memcpy(x, &(c.vector));
213
214 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);
215
216 /* Compute residual = b - A x = Q (Q^T b - R x) */
217
218 gsl_vector_set_zero(&(c.vector));
219
220 gsl_linalg_QR_Qvec(QR, tau, residual);
221
222 return GSL_SUCCESS;
223 }
224 }
225
226
227 int
gsl_linalg_QR_Rsolve(const gsl_matrix * QR,const gsl_vector * b,gsl_vector * x)228 gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
229 {
230 if (QR->size1 != QR->size2)
231 {
232 GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
233 }
234 else if (QR->size1 != b->size)
235 {
236 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
237 }
238 else if (QR->size2 != x->size)
239 {
240 GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
241 }
242 else
243 {
244 /* Copy x <- b */
245
246 gsl_vector_memcpy (x, b);
247
248 /* Solve R x = b, storing x in-place */
249
250 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
251
252 return GSL_SUCCESS;
253 }
254 }
255
256
257 int
gsl_linalg_QR_Rsvx(const gsl_matrix * QR,gsl_vector * x)258 gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
259 {
260 if (QR->size1 != QR->size2)
261 {
262 GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
263 }
264 else if (QR->size1 != x->size)
265 {
266 GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
267 }
268 else
269 {
270 /* Solve R x = b, storing x in-place */
271
272 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
273
274 return GSL_SUCCESS;
275 }
276 }
277
278 int
gsl_linalg_R_solve(const gsl_matrix * R,const gsl_vector * b,gsl_vector * x)279 gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
280 {
281 if (R->size1 != R->size2)
282 {
283 GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
284 }
285 else if (R->size1 != b->size)
286 {
287 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
288 }
289 else if (R->size2 != x->size)
290 {
291 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
292 }
293 else
294 {
295 /* Copy x <- b */
296
297 gsl_vector_memcpy (x, b);
298
299 /* Solve R x = b, storing x inplace in b */
300
301 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
302
303 return GSL_SUCCESS;
304 }
305 }
306
307 int
gsl_linalg_R_svx(const gsl_matrix * R,gsl_vector * x)308 gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x)
309 {
310 if (R->size1 != R->size2)
311 {
312 GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
313 }
314 else if (R->size2 != x->size)
315 {
316 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
317 }
318 else
319 {
320 /* Solve R x = b, storing x inplace in b */
321
322 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
323
324 return GSL_SUCCESS;
325 }
326 }
327
328
329
330 /* Form the product Q^T v from a QR factorized matrix
331 */
332
333 int
gsl_linalg_QR_QTvec(const gsl_matrix * QR,const gsl_vector * tau,gsl_vector * v)334 gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
335 {
336 const size_t M = QR->size1;
337 const size_t N = QR->size2;
338
339 if (tau->size != GSL_MIN (M, N))
340 {
341 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
342 }
343 else if (v->size != M)
344 {
345 GSL_ERROR ("vector size must be M", GSL_EBADLEN);
346 }
347 else
348 {
349 size_t i;
350
351 /* compute Q^T v */
352
353 for (i = 0; i < GSL_MIN (M, N); i++)
354 {
355 gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
356 gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
357 gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
358 double ti = gsl_vector_get (tau, i);
359 gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
360 }
361 return GSL_SUCCESS;
362 }
363 }
364
365
366 int
gsl_linalg_QR_Qvec(const gsl_matrix * QR,const gsl_vector * tau,gsl_vector * v)367 gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
368 {
369 const size_t M = QR->size1;
370 const size_t N = QR->size2;
371
372 if (tau->size != GSL_MIN (M, N))
373 {
374 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
375 }
376 else if (v->size != M)
377 {
378 GSL_ERROR ("vector size must be M", GSL_EBADLEN);
379 }
380 else
381 {
382 size_t i;
383
384 /* compute Q^T v */
385
386 for (i = GSL_MIN (M, N); i-- > 0;)
387 {
388 gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
389 gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
390 i, M - i);
391 gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
392 double ti = gsl_vector_get (tau, i);
393 gsl_linalg_householder_hv (ti, &h.vector, &w.vector);
394 }
395 return GSL_SUCCESS;
396 }
397 }
398
399 /* Form the product Q^T A from a QR factorized matrix */
400
401 int
gsl_linalg_QR_QTmat(const gsl_matrix * QR,const gsl_vector * tau,gsl_matrix * A)402 gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
403 {
404 const size_t M = QR->size1;
405 const size_t N = QR->size2;
406
407 if (tau->size != GSL_MIN (M, N))
408 {
409 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
410 }
411 else if (A->size1 != M)
412 {
413 GSL_ERROR ("matrix must have M rows", GSL_EBADLEN);
414 }
415 else
416 {
417 size_t i;
418
419 /* compute Q^T A */
420
421 for (i = 0; i < GSL_MIN (M, N); i++)
422 {
423 gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
424 gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
425 gsl_matrix_view m = gsl_matrix_submatrix(A, i, 0, M - i, A->size2);
426 double ti = gsl_vector_get (tau, i);
427 gsl_linalg_householder_hm (ti, &(h.vector), &(m.matrix));
428 }
429 return GSL_SUCCESS;
430 }
431 }
432
433
434 /* Form the orthogonal matrix Q from the packed QR matrix */
435
436 int
gsl_linalg_QR_unpack(const gsl_matrix * QR,const gsl_vector * tau,gsl_matrix * Q,gsl_matrix * R)437 gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
438 {
439 const size_t M = QR->size1;
440 const size_t N = QR->size2;
441
442 if (Q->size1 != M || Q->size2 != M)
443 {
444 GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
445 }
446 else if (R->size1 != M || R->size2 != N)
447 {
448 GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
449 }
450 else if (tau->size != GSL_MIN (M, N))
451 {
452 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
453 }
454 else
455 {
456 size_t i, j;
457
458 /* Initialize Q to the identity */
459
460 gsl_matrix_set_identity (Q);
461
462 for (i = GSL_MIN (M, N); i-- > 0;)
463 {
464 gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
465 gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
466 i, M - i);
467 gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
468 double ti = gsl_vector_get (tau, i);
469 gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
470 }
471
472 /* Form the right triangular matrix R from a packed QR matrix */
473
474 for (i = 0; i < M; i++)
475 {
476 for (j = 0; j < i && j < N; j++)
477 gsl_matrix_set (R, i, j, 0.0);
478
479 for (j = i; j < N; j++)
480 gsl_matrix_set (R, i, j, gsl_matrix_get (QR, i, j));
481 }
482
483 return GSL_SUCCESS;
484 }
485 }
486
487
488 /* Update a QR factorisation for A= Q R , A' = A + u v^T,
489
490 * Q' R' = QR + u v^T
491 * = Q (R + Q^T u v^T)
492 * = Q (R + w v^T)
493 *
494 * where w = Q^T u.
495 *
496 * Algorithm from Golub and Van Loan, "Matrix Computations", Section
497 * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
498 */
499
500 int
gsl_linalg_QR_update(gsl_matrix * Q,gsl_matrix * R,gsl_vector * w,const gsl_vector * v)501 gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
502 gsl_vector * w, const gsl_vector * v)
503 {
504 const size_t M = R->size1;
505 const size_t N = R->size2;
506
507 if (Q->size1 != M || Q->size2 != M)
508 {
509 GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
510 }
511 else if (w->size != M)
512 {
513 GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
514 }
515 else if (v->size != N)
516 {
517 GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
518 }
519 else
520 {
521 size_t j, k;
522 double w0;
523
524 /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
525
526 J_1^T .... J_(n-1)^T w = +/- |w| e_1
527
528 simultaneously applied to R, H = J_1^T ... J^T_(n-1) R
529 so that H is upper Hessenberg. (12.5.2) */
530
531 for (k = M - 1; k > 0; k--) /* loop from k = M-1 to 1 */
532 {
533 double c, s;
534 double wk = gsl_vector_get (w, k);
535 double wkm1 = gsl_vector_get (w, k - 1);
536
537 create_givens (wkm1, wk, &c, &s);
538 apply_givens_vec (w, k - 1, k, c, s);
539 apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
540 }
541
542 w0 = gsl_vector_get (w, 0);
543
544 /* Add in w v^T (Equation 12.5.3) */
545
546 for (j = 0; j < N; j++)
547 {
548 double r0j = gsl_matrix_get (R, 0, j);
549 double vj = gsl_vector_get (v, j);
550 gsl_matrix_set (R, 0, j, r0j + w0 * vj);
551 }
552
553 /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
554 Equation 12.5.4 */
555
556 for (k = 1; k < GSL_MIN(M,N+1); k++)
557 {
558 double c, s;
559 double diag = gsl_matrix_get (R, k - 1, k - 1);
560 double offdiag = gsl_matrix_get (R, k, k - 1);
561
562 create_givens (diag, offdiag, &c, &s);
563 apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
564
565 gsl_matrix_set (R, k, k - 1, 0.0); /* exact zero of G^T */
566 }
567
568 return GSL_SUCCESS;
569 }
570 }
571
572 int
gsl_linalg_QR_QRsolve(gsl_matrix * Q,gsl_matrix * R,const gsl_vector * b,gsl_vector * x)573 gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
574 {
575 const size_t M = R->size1;
576 const size_t N = R->size2;
577
578 if (M != N)
579 {
580 return GSL_ENOTSQR;
581 }
582 else if (Q->size1 != M || b->size != M || x->size != M)
583 {
584 return GSL_EBADLEN;
585 }
586 else
587 {
588 /* compute sol = Q^T b */
589
590 gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);
591
592 /* Solve R x = sol, storing x in-place */
593
594 gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
595
596 return GSL_SUCCESS;
597 }
598 }
599