3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
4 * Copyright (C) 2004 Joerg Wensch, modifications for LQ.
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.
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.
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.
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_vector.h>
26 #include <gsl/gsl_matrix.h>
27 #include <gsl/gsl_blas.h>
29 #include <gsl/gsl_linalg.h>
34 #include "apply_givens.c"
36 /* Note: The standard in numerical linear algebra is to solve A x = b
37 * resp. ||A x - b||_2 -> min by QR-decompositions where x, b are
40 * When the matrix A has a large number of rows it is much more
41 * efficient to work with the transposed matrix A^T and to solve the
42 * system x^T A = b^T resp. ||x^T A - b^T||_2 -> min. This is caused
43 * by the row-oriented format in which GSL stores matrices. Therefore
44 * the QR-decomposition of A has to be replaced by a LQ decomposition
47 * The purpose of this package is to provide the algorithms to compute
48 * the LQ-decomposition and to solve the linear equations resp. least
49 * squares problems. The dimensions N, M of the matrix are switched
50 * because here A will probably be a transposed matrix. We write x^T,
51 * b^T,... for vectors the comments to emphasize that they are row
54 * It may even be useful to transpose your matrix explicitly (assumed
55 * that there are no memory restrictions) because this takes O(M x N)
56 * computing time where the decompostion takes O(M x N^2) computing
59 /* Factorise a general N x M matrix A into
63 * where Q is orthogonal (M x M) and L is lower triangular (N x M).
65 * Q is stored as a packed set of Householder transformations in the
66 * strict upper triangular part of the input matrix.
68 * R is stored in the diagonal and lower triangle of the input matrix.
70 * The full matrix for Q can be obtained as the product
74 * where k = MIN(M,N) and
76 * Q_i = (I - tau_i * v_i * v_i')
78 * and where v_i is a Householder vector
80 * v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
82 * This storage scheme is the same as in LAPACK. */
85 gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau)
87 const size_t N = A->size1;
88 const size_t M = A->size2;
90 if (tau->size != GSL_MIN (M, N))
92 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
98 for (i = 0; i < GSL_MIN (M, N); i++)
100 /* Compute the Householder transformation to reduce the j-th
101 column of the matrix to a multiple of the j-th unit vector */
103 gsl_vector_view c_full = gsl_matrix_row (A, i);
104 gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
106 double tau_i = gsl_linalg_householder_transform (&(c.vector));
108 gsl_vector_set (tau, i, tau_i);
110 /* Apply the transformation to the remaining columns and
115 gsl_matrix_view m = gsl_matrix_submatrix (A, i + 1, i, N - (i + 1), M - i );
116 gsl_linalg_householder_mh (tau_i, &(c.vector), &(m.matrix));
124 /* Solves the system x^T A = b^T using the LQ factorisation,
128 * to obtain x. Based on SLATEC code.
133 gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
135 if (LQ->size1 != LQ->size2)
137 GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
139 else if (LQ->size2 != b->size)
141 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
143 else if (LQ->size1 != x->size)
145 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
151 gsl_vector_memcpy (x, b);
155 gsl_linalg_LQ_svx_T (LQ, tau, x);
161 /* Solves the system x^T A = b^T in place using the LQ factorisation,
165 * to obtain x. Based on SLATEC code.
169 gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * x)
172 if (LQ->size1 != LQ->size2)
174 GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
176 else if (LQ->size1 != x->size)
178 GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
182 /* compute rhs = Q^T b */
184 gsl_linalg_LQ_vecQT (LQ, tau, x);
186 /* Solve R x = rhs, storing x in-place */
188 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
195 /* Find the least squares solution to the overdetermined system
199 * for M >= N using the LQ factorization A = L Q.
203 gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
205 const size_t N = LQ->size1;
206 const size_t M = LQ->size2;
210 GSL_ERROR ("LQ matrix must have M>=N", GSL_EBADLEN);
212 else if (M != b->size)
214 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
216 else if (N != x->size)
218 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
220 else if (M != residual->size)
222 GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
226 gsl_matrix_const_view L = gsl_matrix_const_submatrix (LQ, 0, 0, N, N);
227 gsl_vector_view c = gsl_vector_subvector(residual, 0, N);
229 gsl_vector_memcpy(residual, b);
231 /* compute rhs = b^T Q^T */
233 gsl_linalg_LQ_vecQT (LQ, tau, residual);
235 /* Solve x^T L = rhs */
237 gsl_vector_memcpy(x, &(c.vector));
239 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, &(L.matrix), x);
241 /* Compute residual = b^T - x^T A = (b^T Q^T - x^T L) Q */
243 gsl_vector_set_zero(&(c.vector));
245 gsl_linalg_LQ_vecQ(LQ, tau, residual);
253 gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b, gsl_vector * x)
255 if (LQ->size1 != LQ->size2)
257 GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
259 else if (LQ->size1 != b->size)
261 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
263 else if (LQ->size1 != x->size)
265 GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
271 gsl_vector_memcpy (x, b);
273 /* Solve R x = b, storing x in-place */
275 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
283 gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x)
285 if (LQ->size1 != LQ->size2)
287 GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
289 else if (LQ->size2 != x->size)
291 GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
295 /* Solve x^T L = b^T, storing x in-place */
297 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
304 gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b, gsl_vector * x)
306 if (L->size1 != L->size2)
308 GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
310 else if (L->size2 != b->size)
312 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
314 else if (L->size1 != x->size)
316 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
322 gsl_vector_memcpy (x, b);
324 /* Solve R x = b, storing x inplace in b */
326 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);
336 gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v)
338 const size_t N = LQ->size1;
339 const size_t M = LQ->size2;
341 if (tau->size != GSL_MIN (M, N))
343 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
345 else if (v->size != M)
347 GSL_ERROR ("vector size must be M", GSL_EBADLEN);
355 for (i = 0; i < GSL_MIN (M, N); i++)
357 gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
358 gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
360 gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
361 double ti = gsl_vector_get (tau, i);
362 gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
369 gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v)
371 const size_t N = LQ->size1;
372 const size_t M = LQ->size2;
374 if (tau->size != GSL_MIN (M, N))
376 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
378 else if (v->size != M)
380 GSL_ERROR ("vector size must be M", GSL_EBADLEN);
388 for (i = GSL_MIN (M, N); i > 0 && i--;)
390 gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
391 gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
393 gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
394 double ti = gsl_vector_get (tau, i);
395 gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
402 /* Form the orthogonal matrix Q from the packed LQ matrix */
405 gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L)
407 const size_t N = LQ->size1;
408 const size_t M = LQ->size2;
410 if (Q->size1 != M || Q->size2 != M)
412 GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
414 else if (L->size1 != N || L->size2 != M)
416 GSL_ERROR ("R matrix must be N x M", GSL_ENOTSQR);
418 else if (tau->size != GSL_MIN (M, N))
420 GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
424 size_t i, j, l_border;
426 /* Initialize Q to the identity */
428 gsl_matrix_set_identity (Q);
430 for (i = GSL_MIN (M, N); i > 0 && i--;)
432 gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
433 gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
435 gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
436 double ti = gsl_vector_get (tau, i);
437 gsl_linalg_householder_mh (ti, &h.vector, &m.matrix);
440 /* Form the lower triangular matrix L from a packed LQ matrix */
442 for (i = 0; i < N; i++)
444 l_border=GSL_MIN(i,M-1);
445 for (j = 0; j <= l_border ; j++)
446 gsl_matrix_set (L, i, j, gsl_matrix_get (LQ, i, j));
448 for (j = l_border+1; j < M; j++)
449 gsl_matrix_set (L, i, j, 0.0);
457 /* Update a LQ factorisation for A= L Q , A' = A + v u^T,
460 * = (L + v u^T Q^T) Q
465 * Algorithm from Golub and Van Loan, "Matrix Computations", Section
466 * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
470 gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * L,
471 const gsl_vector * v, gsl_vector * w)
473 const size_t N = L->size1;
474 const size_t M = L->size2;
476 if (Q->size1 != M || Q->size2 != M)
478 GSL_ERROR ("Q matrix must be N x N if L is M x N", GSL_ENOTSQR);
480 else if (w->size != M)
482 GSL_ERROR ("w must be length N if L is M x N", GSL_EBADLEN);
484 else if (v->size != N)
486 GSL_ERROR ("v must be length M if L is M x N", GSL_EBADLEN);
493 /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
495 J_1^T .... J_(n-1)^T w = +/- |w| e_1
497 simultaneously applied to L, H = J_1^T ... J^T_(n-1) L
498 so that H is upper Hessenberg. (12.5.2) */
500 for (k = M - 1; k > 0; k--)
503 double wk = gsl_vector_get (w, k);
504 double wkm1 = gsl_vector_get (w, k - 1);
506 create_givens (wkm1, wk, &c, &s);
507 apply_givens_vec (w, k - 1, k, c, s);
508 apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
511 w0 = gsl_vector_get (w, 0);
513 /* Add in v w^T (Equation 12.5.3) */
515 for (j = 0; j < N; j++)
517 double lj0 = gsl_matrix_get (L, j, 0);
518 double vj = gsl_vector_get (v, j);
519 gsl_matrix_set (L, j, 0, lj0 + w0 * vj);
522 /* Apply Givens transformations L' = G_(n-1)^T ... G_1^T H
525 for (k = 1; k < GSL_MIN(M,N+1); k++)
528 double diag = gsl_matrix_get (L, k - 1, k - 1);
529 double offdiag = gsl_matrix_get (L, k - 1 , k);
531 create_givens (diag, offdiag, &c, &s);
532 apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
534 gsl_matrix_set (L, k - 1, k, 0.0); /* exact zero of G^T */
542 gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L, const gsl_vector * b, gsl_vector * x)
544 const size_t N = L->size1;
545 const size_t M = L->size2;
551 else if (Q->size1 != M || b->size != M || x->size != M)
557 /* compute sol = b^T Q^T */
559 gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
561 /* Solve x^T L = sol, storing x in-place */
563 gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);