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