htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/lq.c

   1 /* linalg/lq.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
   4  * Copyright (C) 2004 Joerg Wensch, modifications for LQ.
   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 #include <config.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_vector.h>
  26 #include <gsl/gsl_matrix.h>
  27 #include <gsl/gsl_blas.h>
  28
  29 #include <gsl/gsl_linalg.h>
  30
  31 #define REAL double
  32
  33 #include "givens.c"
  34 #include "apply_givens.c"
  35
  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
  38  * column vectors.
  39  *
  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
  45  * of A^T
  46  *
  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
  52  * vectors.
  53  *
  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
  57  * time.  */
  58
  59 /* Factorise a general N x M matrix A into
  60  *
  61  *   A = L Q
  62  *
  63  * where Q is orthogonal (M x M) and L is lower triangular (N x M).
  64  *
  65  * Q is stored as a packed set of Householder transformations in the
  66  * strict upper triangular part of the input matrix.
  67  *
  68  * R is stored in the diagonal and lower triangle of the input matrix.
  69  *
  70  * The full matrix for Q can be obtained as the product
  71  *
  72  *       Q = Q_k .. Q_2 Q_1
  73  *
  74  * where k = MIN(M,N) and
  75  *
  76  *       Q_i = (I - tau_i * v_i * v_i')
  77  *
  78  * and where v_i is a Householder vector
  79  *
  80  *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
  81  *
  82  * This storage scheme is the same as in LAPACK.  */
  83
  84 int
  85 gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau)
  86 {
  87   const size_t N = A->size1;
  88   const size_t M = A->size2;
  89
  90   if (tau->size != GSL_MIN (M, N))
  91     {
  92       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  93     }
  94   else
  95     {
  96       size_t i;
  97
  98       for (i = 0; i < GSL_MIN (M, N); i++)
  99         {
 100           /* Compute the Householder transformation to reduce the j-th
 101              column of the matrix to a multiple of the j-th unit vector */
 102
 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);
 105
 106           double tau_i = gsl_linalg_householder_transform (&(c.vector));
 107
 108           gsl_vector_set (tau, i, tau_i);
 109
 110           /* Apply the transformation to the remaining columns and
 111              update the norms */
 112
 113           if (i + 1 < N)
 114             {
 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));
 117             }
 118         }
 119
 120       return GSL_SUCCESS;
 121     }
 122 }
 123
 124 /* Solves the system x^T A = b^T using the LQ factorisation,
 125
 126  *  x^T L = b^T Q^T
 127  *
 128  * to obtain x. Based on SLATEC code.
 129  */
 130
 131
 132 int
 133 gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
 134 {
 135   if (LQ->size1 != LQ->size2)
 136     {
 137       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 138     }
 139   else if (LQ->size2 != b->size)
 140     {
 141       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 142     }
 143   else if (LQ->size1 != x->size)
 144     {
 145       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 146     }
 147   else
 148     {
 149       /* Copy x <- b */
 150
 151       gsl_vector_memcpy (x, b);
 152
 153       /* Solve for x */
 154
 155       gsl_linalg_LQ_svx_T (LQ, tau, x);
 156
 157       return GSL_SUCCESS;
 158     }
 159 }
 160
 161 /* Solves the system x^T A = b^T in place using the LQ factorisation,
 162  *
 163  *  x^T L = b^T Q^T
 164  *
 165  * to obtain x. Based on SLATEC code.
 166  */
 167
 168 int
 169 gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * x)
 170 {
 171
 172   if (LQ->size1 != LQ->size2)
 173     {
 174       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 175     }
 176   else if (LQ->size1 != x->size)
 177     {
 178       GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
 179     }
 180   else
 181     {
 182       /* compute rhs = Q^T b */
 183
 184       gsl_linalg_LQ_vecQT (LQ, tau, x);
 185
 186       /* Solve R x = rhs, storing x in-place */
 187
 188       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 189
 190       return GSL_SUCCESS;
 191     }
 192 }
 193
 194
 195 /* Find the least squares solution to the overdetermined system
 196  *
 197  *   x^T A = b^T
 198  *
 199  * for M >= N using the LQ factorization A = L Q.
 200  */
 201
 202 int
 203 gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
 204 {
 205   const size_t N = LQ->size1;
 206   const size_t M = LQ->size2;
 207
 208   if (M < N)
 209     {
 210       GSL_ERROR ("LQ matrix must have M>=N", GSL_EBADLEN);
 211     }
 212   else if (M != b->size)
 213     {
 214       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 215     }
 216   else if (N != x->size)
 217     {
 218       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 219     }
 220   else if (M != residual->size)
 221     {
 222       GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
 223     }
 224   else
 225     {
 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);
 228
 229       gsl_vector_memcpy(residual, b);
 230
 231       /* compute rhs = b^T Q^T */
 232
 233       gsl_linalg_LQ_vecQT (LQ, tau, residual);
 234
 235       /* Solve x^T L = rhs */
 236
 237       gsl_vector_memcpy(x, &(c.vector));
 238
 239       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, &(L.matrix), x);
 240
 241       /* Compute residual = b^T - x^T A = (b^T Q^T - x^T L) Q */
 242
 243       gsl_vector_set_zero(&(c.vector));
 244
 245       gsl_linalg_LQ_vecQ(LQ, tau, residual);
 246
 247       return GSL_SUCCESS;
 248     }
 249 }
 250
 251
 252 int
 253 gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b, gsl_vector * x)
 254 {
 255   if (LQ->size1 != LQ->size2)
 256     {
 257       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 258     }
 259   else if (LQ->size1 != b->size)
 260     {
 261       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 262     }
 263   else if (LQ->size1 != x->size)
 264     {
 265       GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
 266     }
 267   else
 268     {
 269       /* Copy x <- b */
 270
 271       gsl_vector_memcpy (x, b);
 272
 273       /* Solve R x = b, storing x in-place */
 274
 275       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 276
 277       return GSL_SUCCESS;
 278     }
 279 }
 280
 281
 282 int
 283 gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x)
 284 {
 285   if (LQ->size1 != LQ->size2)
 286     {
 287       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 288     }
 289   else if (LQ->size2 != x->size)
 290     {
 291       GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
 292     }
 293   else
 294     {
 295       /* Solve x^T L = b^T, storing x in-place */
 296
 297       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 298
 299       return GSL_SUCCESS;
 300     }
 301 }
 302
 303 int
 304 gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b, gsl_vector * x)
 305 {
 306   if (L->size1 != L->size2)
 307     {
 308       GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
 309     }
 310   else if (L->size2 != b->size)
 311     {
 312       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 313     }
 314   else if (L->size1 != x->size)
 315     {
 316       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 317     }
 318   else
 319     {
 320       /* Copy x <- b */
 321
 322       gsl_vector_memcpy (x, b);
 323
 324       /* Solve R x = b, storing x inplace in b */
 325
 326       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);
 327
 328       return GSL_SUCCESS;
 329     }
 330 }
 331
 332
 333
 334
 335 int
 336 gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v)
 337 {
 338   const size_t N = LQ->size1;
 339   const size_t M = LQ->size2;
 340
 341   if (tau->size != GSL_MIN (M, N))
 342     {
 343       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 344     }
 345   else if (v->size != M)
 346     {
 347       GSL_ERROR ("vector size must be M", GSL_EBADLEN);
 348     }
 349   else
 350     {
 351       size_t i;
 352
 353       /* compute v Q^T  */
 354
 355       for (i = 0; i < GSL_MIN (M, N); i++)
 356         {
 357           gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
 358           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
 359                                                                 i, M - i);
 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));
 363         }
 364       return GSL_SUCCESS;
 365     }
 366 }
 367
 368 int
 369 gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v)
 370 {
 371   const size_t N = LQ->size1;
 372   const size_t M = LQ->size2;
 373
 374   if (tau->size != GSL_MIN (M, N))
 375     {
 376       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 377     }
 378   else if (v->size != M)
 379     {
 380       GSL_ERROR ("vector size must be M", GSL_EBADLEN);
 381     }
 382   else
 383     {
 384       size_t i;
 385
 386       /* compute v Q^T  */
 387
 388       for (i =  GSL_MIN (M, N); i > 0 && i--;)
 389         {
 390           gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
 391           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
 392                                                                 i, M - i);
 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));
 396         }
 397       return GSL_SUCCESS;
 398     }
 399 }
 400
 401
 402 /*  Form the orthogonal matrix Q from the packed LQ matrix */
 403
 404 int
 405 gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L)
 406 {
 407   const size_t N = LQ->size1;
 408   const size_t M = LQ->size2;
 409
 410   if (Q->size1 != M || Q->size2 != M)
 411     {
 412       GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
 413     }
 414   else if (L->size1 != N || L->size2 != M)
 415     {
 416       GSL_ERROR ("R matrix must be N x M", GSL_ENOTSQR);
 417     }
 418   else if (tau->size != GSL_MIN (M, N))
 419     {
 420       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 421     }
 422   else
 423     {
 424       size_t i, j, l_border;
 425
 426       /* Initialize Q to the identity */
 427
 428       gsl_matrix_set_identity (Q);
 429
 430       for (i = GSL_MIN (M, N); i > 0 && i--;)
 431         {
 432           gsl_vector_const_view c = gsl_matrix_const_row (LQ, i);
 433           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
 434                                                                 i, M - i);
 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);
 438         }
 439
 440       /*  Form the lower triangular matrix L from a packed LQ matrix */
 441
 442       for (i = 0; i < N; i++)
 443         {
 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));
 447
 448             for (j = l_border+1; j < M; j++)
 449                 gsl_matrix_set (L, i, j, 0.0);
 450         }
 451
 452       return GSL_SUCCESS;
 453     }
 454 }
 455
 456
 457 /* Update a LQ factorisation for A= L Q ,  A' = A + v u^T,
 458
 459  * L' Q' = LQ + v u^T
 460  *       = (L + v u^T Q^T) Q
 461  *       = (L + v w^T) Q
 462  *
 463  * where w = Q u.
 464  *
 465  * Algorithm from Golub and Van Loan, "Matrix Computations", Section
 466  * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
 467  */
 468
 469 int
 470 gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * L,
 471                       const gsl_vector * v, gsl_vector * w)
 472 {
 473   const size_t N = L->size1;
 474   const size_t M = L->size2;
 475
 476   if (Q->size1 != M || Q->size2 != M)
 477     {
 478       GSL_ERROR ("Q matrix must be N x N if L is M x N", GSL_ENOTSQR);
 479     }
 480   else if (w->size != M)
 481     {
 482       GSL_ERROR ("w must be length N if L is M x N", GSL_EBADLEN);
 483     }
 484   else if (v->size != N)
 485     {
 486       GSL_ERROR ("v must be length M if L is M x N", GSL_EBADLEN);
 487     }
 488   else
 489     {
 490       size_t j, k;
 491       double w0;
 492
 493       /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
 494
 495          J_1^T .... J_(n-1)^T w = +/- |w| e_1
 496
 497          simultaneously applied to L,  H = J_1^T ... J^T_(n-1) L
 498          so that H is upper Hessenberg.  (12.5.2) */
 499
 500       for (k = M - 1; k > 0; k--)
 501         {
 502           double c, s;
 503           double wk = gsl_vector_get (w, k);
 504           double wkm1 = gsl_vector_get (w, k - 1);
 505
 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);
 509        }
 510
 511       w0 = gsl_vector_get (w, 0);
 512
 513       /* Add in v w^T  (Equation 12.5.3) */
 514
 515       for (j = 0; j < N; j++)
 516         {
 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);
 520         }
 521
 522       /* Apply Givens transformations L' = G_(n-1)^T ... G_1^T H
 523          Equation 12.5.4 */
 524
 525       for (k = 1; k < GSL_MIN(M,N+1); k++)
 526         {
 527           double c, s;
 528           double diag = gsl_matrix_get (L, k - 1, k - 1);
 529           double offdiag = gsl_matrix_get (L, k - 1 , k);
 530
 531           create_givens (diag, offdiag, &c, &s);
 532           apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
 533
 534           gsl_matrix_set (L, k - 1, k, 0.0);    /* exact zero of G^T */
 535         }
 536
 537       return GSL_SUCCESS;
 538     }
 539 }
 540
 541 int
 542 gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L, const gsl_vector * b, gsl_vector * x)
 543 {
 544   const size_t N = L->size1;
 545   const size_t M = L->size2;
 546
 547   if (M != N)
 548     {
 549       return GSL_ENOTSQR;
 550     }
 551   else if (Q->size1 != M || b->size != M || x->size != M)
 552     {
 553       return GSL_EBADLEN;
 554     }
 555   else
 556     {
 557       /* compute sol = b^T Q^T */
 558
 559       gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
 560
 561       /* Solve x^T L = sol, storing x in-place */
 562
 563       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);
 564
 565       return GSL_SUCCESS;
 566     }
 567 }