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

   1 /* linalg/svd.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_vector.h>
  25 #include <gsl/gsl_matrix.h>
  26 #include <gsl/gsl_blas.h>
  27
  28 #include <gsl/gsl_linalg.h>
  29
  30 #include "givens.c"
  31 #include "svdstep.c"
  32
  33 /* Factorise a general M x N matrix A into,
  34  *
  35  *   A = U D V^T
  36  *
  37  * where U is a column-orthogonal M x N matrix (U^T U = I),
  38  * D is a diagonal N x N matrix,
  39  * and V is an N x N orthogonal matrix (V^T V = V V^T = I)
  40  *
  41  * U is stored in the original matrix A, which has the same size
  42  *
  43  * V is stored as a separate matrix (not V^T). You must take the
  44  * transpose to form the product above.
  45  *
  46  * The diagonal matrix D is stored in the vector S,  D_ii = S_i
  47  */
  48
  49 int
  50 gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S,
  51                       gsl_vector * work)
  52 {
  53   size_t a, b, i, j, iter;
  54
  55   const size_t M = A->size1;
  56   const size_t N = A->size2;
  57   const size_t K = GSL_MIN (M, N);
  58
  59   if (M < N)
  60     {
  61       GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
  62     }
  63   else if (V->size1 != N)
  64     {
  65       GSL_ERROR ("square matrix V must match second dimension of matrix A",
  66                  GSL_EBADLEN);
  67     }
  68   else if (V->size1 != V->size2)
  69     {
  70       GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
  71     }
  72   else if (S->size != N)
  73     {
  74       GSL_ERROR ("length of vector S must match second dimension of matrix A",
  75                  GSL_EBADLEN);
  76     }
  77   else if (work->size != N)
  78     {
  79       GSL_ERROR ("length of workspace must match second dimension of matrix A",
  80                  GSL_EBADLEN);
  81     }
  82
  83   /* Handle the case of N = 1 (SVD of a column vector) */
  84
  85   if (N == 1)
  86     {
  87       gsl_vector_view column = gsl_matrix_column (A, 0);
  88       double norm = gsl_blas_dnrm2 (&column.vector);
  89
  90       gsl_vector_set (S, 0, norm);
  91       gsl_matrix_set (V, 0, 0, 1.0);
  92
  93       if (norm != 0.0)
  94         {
  95           gsl_blas_dscal (1.0/norm, &column.vector);
  96         }
  97
  98       return GSL_SUCCESS;
  99     }
 100
 101   {
 102     gsl_vector_view f = gsl_vector_subvector (work, 0, K - 1);
 103
 104     /* bidiagonalize matrix A, unpack A into U S V */
 105
 106     gsl_linalg_bidiag_decomp (A, S, &f.vector);
 107     gsl_linalg_bidiag_unpack2 (A, S, &f.vector, V);
 108
 109     /* apply reduction steps to B=(S,Sd) */
 110
 111     chop_small_elements (S, &f.vector);
 112
 113     /* Progressively reduce the matrix until it is diagonal */
 114
 115     b = N - 1;
 116     iter = 0;
 117
 118     while (b > 0)
 119       {
 120         double fbm1 = gsl_vector_get (&f.vector, b - 1);
 121
 122         if (fbm1 == 0.0 || gsl_isnan (fbm1))
 123           {
 124             b--;
 125             continue;
 126           }
 127
 128         /* Find the largest unreduced block (a,b) starting from b
 129            and working backwards */
 130
 131         a = b - 1;
 132
 133         while (a > 0)
 134           {
 135             double fam1 = gsl_vector_get (&f.vector, a - 1);
 136
 137             if (fam1 == 0.0 || gsl_isnan (fam1))
 138               {
 139                 break;
 140               }
 141
 142             a--;
 143           }
 144
 145         iter++;
 146
 147         if (iter > 100 * N)
 148           {
 149             GSL_ERROR("SVD decomposition failed to converge", GSL_EMAXITER);
 150           }
 151
 152
 153         {
 154           const size_t n_block = b - a + 1;
 155           gsl_vector_view S_block = gsl_vector_subvector (S, a, n_block);
 156           gsl_vector_view f_block = gsl_vector_subvector (&f.vector, a, n_block - 1);
 157
 158           gsl_matrix_view U_block =
 159             gsl_matrix_submatrix (A, 0, a, A->size1, n_block);
 160           gsl_matrix_view V_block =
 161             gsl_matrix_submatrix (V, 0, a, V->size1, n_block);
 162
 163           qrstep (&S_block.vector, &f_block.vector, &U_block.matrix, &V_block.matrix);
 164
 165           /* remove any small off-diagonal elements */
 166
 167           chop_small_elements (&S_block.vector, &f_block.vector);
 168         }
 169       }
 170   }
 171   /* Make singular values positive by reflections if necessary */
 172
 173   for (j = 0; j < K; j++)
 174     {
 175       double Sj = gsl_vector_get (S, j);
 176
 177       if (Sj < 0.0)
 178         {
 179           for (i = 0; i < N; i++)
 180             {
 181               double Vij = gsl_matrix_get (V, i, j);
 182               gsl_matrix_set (V, i, j, -Vij);
 183             }
 184
 185           gsl_vector_set (S, j, -Sj);
 186         }
 187     }
 188
 189   /* Sort singular values into decreasing order */
 190
 191   for (i = 0; i < K; i++)
 192     {
 193       double S_max = gsl_vector_get (S, i);
 194       size_t i_max = i;
 195
 196       for (j = i + 1; j < K; j++)
 197         {
 198           double Sj = gsl_vector_get (S, j);
 199
 200           if (Sj > S_max)
 201             {
 202               S_max = Sj;
 203               i_max = j;
 204             }
 205         }
 206
 207       if (i_max != i)
 208         {
 209           /* swap eigenvalues */
 210           gsl_vector_swap_elements (S, i, i_max);
 211
 212           /* swap eigenvectors */
 213           gsl_matrix_swap_columns (A, i, i_max);
 214           gsl_matrix_swap_columns (V, i, i_max);
 215         }
 216     }
 217
 218   return GSL_SUCCESS;
 219 }
 220
 221
 222 /* Modified algorithm which is better for M>>N */
 223
 224 int
 225 gsl_linalg_SV_decomp_mod (gsl_matrix * A,
 226                           gsl_matrix * X,
 227                           gsl_matrix * V, gsl_vector * S, gsl_vector * work)
 228 {
 229   size_t i, j;
 230
 231   const size_t M = A->size1;
 232   const size_t N = A->size2;
 233
 234   if (M < N)
 235     {
 236       GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
 237     }
 238   else if (V->size1 != N)
 239     {
 240       GSL_ERROR ("square matrix V must match second dimension of matrix A",
 241                  GSL_EBADLEN);
 242     }
 243   else if (V->size1 != V->size2)
 244     {
 245       GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
 246     }
 247   else if (X->size1 != N)
 248     {
 249       GSL_ERROR ("square matrix X must match second dimension of matrix A",
 250                  GSL_EBADLEN);
 251     }
 252   else if (X->size1 != X->size2)
 253     {
 254       GSL_ERROR ("matrix X must be square", GSL_ENOTSQR);
 255     }
 256   else if (S->size != N)
 257     {
 258       GSL_ERROR ("length of vector S must match second dimension of matrix A",
 259                  GSL_EBADLEN);
 260     }
 261   else if (work->size != N)
 262     {
 263       GSL_ERROR ("length of workspace must match second dimension of matrix A",
 264                  GSL_EBADLEN);
 265     }
 266
 267   if (N == 1)
 268     {
 269       gsl_vector_view column = gsl_matrix_column (A, 0);
 270       double norm = gsl_blas_dnrm2 (&column.vector);
 271
 272       gsl_vector_set (S, 0, norm);
 273       gsl_matrix_set (V, 0, 0, 1.0);
 274
 275       if (norm != 0.0)
 276         {
 277           gsl_blas_dscal (1.0/norm, &column.vector);
 278         }
 279
 280       return GSL_SUCCESS;
 281     }
 282
 283   /* Convert A into an upper triangular matrix R */
 284
 285   for (i = 0; i < N; i++)
 286     {
 287       gsl_vector_view c = gsl_matrix_column (A, i);
 288       gsl_vector_view v = gsl_vector_subvector (&c.vector, i, M - i);
 289       double tau_i = gsl_linalg_householder_transform (&v.vector);
 290
 291       /* Apply the transformation to the remaining columns */
 292
 293       if (i + 1 < N)
 294         {
 295           gsl_matrix_view m =
 296             gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
 297           gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
 298         }
 299
 300       gsl_vector_set (S, i, tau_i);
 301     }
 302
 303   /* Copy the upper triangular part of A into X */
 304
 305   for (i = 0; i < N; i++)
 306     {
 307       for (j = 0; j < i; j++)
 308         {
 309           gsl_matrix_set (X, i, j, 0.0);
 310         }
 311
 312       {
 313         double Aii = gsl_matrix_get (A, i, i);
 314         gsl_matrix_set (X, i, i, Aii);
 315       }
 316
 317       for (j = i + 1; j < N; j++)
 318         {
 319           double Aij = gsl_matrix_get (A, i, j);
 320           gsl_matrix_set (X, i, j, Aij);
 321         }
 322     }
 323
 324   /* Convert A into an orthogonal matrix L */
 325
 326   for (j = N; j > 0 && j--;)
 327     {
 328       /* Householder column transformation to accumulate L */
 329       double tj = gsl_vector_get (S, j);
 330       gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M - j, N - j);
 331       gsl_linalg_householder_hm1 (tj, &m.matrix);
 332     }
 333
 334   /* unpack R into X V S */
 335
 336   gsl_linalg_SV_decomp (X, V, S, work);
 337
 338   /* Multiply L by X, to obtain U = L X, stored in U */
 339
 340   {
 341     gsl_vector_view sum = gsl_vector_subvector (work, 0, N);
 342
 343     for (i = 0; i < M; i++)
 344       {
 345         gsl_vector_view L_i = gsl_matrix_row (A, i);
 346         gsl_vector_set_zero (&sum.vector);
 347
 348         for (j = 0; j < N; j++)
 349           {
 350             double Lij = gsl_vector_get (&L_i.vector, j);
 351             gsl_vector_view X_j = gsl_matrix_row (X, j);
 352             gsl_blas_daxpy (Lij, &X_j.vector, &sum.vector);
 353           }
 354
 355         gsl_vector_memcpy (&L_i.vector, &sum.vector);
 356       }
 357   }
 358
 359   return GSL_SUCCESS;
 360 }
 361
 362
 363 /*  Solves the system A x = b using the SVD factorization
 364  *
 365  *  A = U S V^T
 366  *
 367  *  to obtain x. For M x N systems it finds the solution in the least
 368  *  squares sense.
 369  */
 370
 371 int
 372 gsl_linalg_SV_solve (const gsl_matrix * U,
 373                      const gsl_matrix * V,
 374                      const gsl_vector * S,
 375                      const gsl_vector * b, gsl_vector * x)
 376 {
 377   if (U->size1 != b->size)
 378     {
 379       GSL_ERROR ("first dimension of matrix U must size of vector b",
 380                  GSL_EBADLEN);
 381     }
 382   else if (U->size2 != S->size)
 383     {
 384       GSL_ERROR ("length of vector S must match second dimension of matrix U",
 385                  GSL_EBADLEN);
 386     }
 387   else if (V->size1 != V->size2)
 388     {
 389       GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
 390     }
 391   else if (S->size != V->size1)
 392     {
 393       GSL_ERROR ("length of vector S must match size of matrix V",
 394                  GSL_EBADLEN);
 395     }
 396   else if (V->size2 != x->size)
 397     {
 398       GSL_ERROR ("size of matrix V must match size of vector x", GSL_EBADLEN);
 399     }
 400   else
 401     {
 402       const size_t N = U->size2;
 403       size_t i;
 404
 405       gsl_vector *w = gsl_vector_calloc (N);
 406
 407       gsl_blas_dgemv (CblasTrans, 1.0, U, b, 0.0, w);
 408
 409       for (i = 0; i < N; i++)
 410         {
 411           double wi = gsl_vector_get (w, i);
 412           double alpha = gsl_vector_get (S, i);
 413           if (alpha != 0)
 414             alpha = 1.0 / alpha;
 415           gsl_vector_set (w, i, alpha * wi);
 416         }
 417
 418       gsl_blas_dgemv (CblasNoTrans, 1.0, V, w, 0.0, x);
 419
 420       gsl_vector_free (w);
 421
 422       return GSL_SUCCESS;
 423     }
 424 }
 425
 426 /* This is a the jacobi version */
 427 /* Author:  G. Jungman */
 428
 429 /*
 430  * Algorithm due to J.C. Nash, Compact Numerical Methods for
 431  * Computers (New York: Wiley and Sons, 1979), chapter 3.
 432  * See also Algorithm 4.1 in
 433  * James Demmel, Kresimir Veselic, "Jacobi's Method is more
 434  * accurate than QR", Lapack Working Note 15 (LAWN15), October 1989.
 435  * Available from netlib.
 436  *
 437  * Based on code by Arthur Kosowsky, Rutgers University
 438  *                  kosowsky@physics.rutgers.edu
 439  *
 440  * Another relevant paper is, P.P.M. De Rijk, "A One-Sided Jacobi
 441  * Algorithm for computing the singular value decomposition on a
 442  * vector computer", SIAM Journal of Scientific and Statistical
 443  * Computing, Vol 10, No 2, pp 359-371, March 1989.
 444  *
 445  */
 446
 447 int
 448 gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * Q, gsl_vector * S)
 449 {
 450   if (A->size1 < A->size2)
 451     {
 452       /* FIXME: only implemented  M>=N case so far */
 453
 454       GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
 455     }
 456   else if (Q->size1 != A->size2)
 457     {
 458       GSL_ERROR ("square matrix Q must match second dimension of matrix A",
 459                  GSL_EBADLEN);
 460     }
 461   else if (Q->size1 != Q->size2)
 462     {
 463       GSL_ERROR ("matrix Q must be square", GSL_ENOTSQR);
 464     }
 465   else if (S->size != A->size2)
 466     {
 467       GSL_ERROR ("length of vector S must match second dimension of matrix A",
 468                  GSL_EBADLEN);
 469     }
 470   else
 471     {
 472       const size_t M = A->size1;
 473       const size_t N = A->size2;
 474       size_t i, j, k;
 475
 476       /* Initialize the rotation counter and the sweep counter. */
 477       int count = 1;
 478       int sweep = 0;
 479       int sweepmax = 5*N;
 480
 481       double tolerance = 10 * M * GSL_DBL_EPSILON;
 482
 483       /* Always do at least 12 sweeps. */
 484       sweepmax = GSL_MAX (sweepmax, 12);
 485
 486       /* Set Q to the identity matrix. */
 487       gsl_matrix_set_identity (Q);
 488
 489       /* Store the column error estimates in S, for use during the
 490          orthogonalization */
 491
 492       for (j = 0; j < N; j++)
 493         {
 494           gsl_vector_view cj = gsl_matrix_column (A, j);
 495           double sj = gsl_blas_dnrm2 (&cj.vector);
 496           gsl_vector_set(S, j, GSL_DBL_EPSILON * sj);
 497         }
 498
 499       /* Orthogonalize A by plane rotations. */
 500
 501       while (count > 0 && sweep <= sweepmax)
 502         {
 503           /* Initialize rotation counter. */
 504           count = N * (N - 1) / 2;
 505
 506           for (j = 0; j < N - 1; j++)
 507             {
 508               for (k = j + 1; k < N; k++)
 509                 {
 510                   double a = 0.0;
 511                   double b = 0.0;
 512                   double p = 0.0;
 513                   double q = 0.0;
 514                   double cosine, sine;
 515                   double v;
 516                   double abserr_a, abserr_b;
 517                   int sorted, orthog, noisya, noisyb;
 518
 519                   gsl_vector_view cj = gsl_matrix_column (A, j);
 520                   gsl_vector_view ck = gsl_matrix_column (A, k);
 521
 522                   gsl_blas_ddot (&cj.vector, &ck.vector, &p);
 523                   p *= 2.0 ;  /* equation 9a:  p = 2 x.y */
 524
 525                   a = gsl_blas_dnrm2 (&cj.vector);
 526                   b = gsl_blas_dnrm2 (&ck.vector);
 527
 528                   q = a * a - b * b;
 529                   v = hypot(p, q);
 530
 531                   /* test for columns j,k orthogonal, or dominant errors */
 532
 533                   abserr_a = gsl_vector_get(S,j);
 534                   abserr_b = gsl_vector_get(S,k);
 535
 536                   sorted = (gsl_coerce_double(a) >= gsl_coerce_double(b));
 537                   orthog = (fabs (p) <= tolerance * gsl_coerce_double(a * b));
 538                   noisya = (a < abserr_a);
 539                   noisyb = (b < abserr_b);
 540
 541                   if (sorted && (orthog || noisya || noisyb))
 542                     {
 543                       count--;
 544                       continue;
 545                     }
 546
 547                   /* calculate rotation angles */
 548                   if (v == 0 || !sorted)
 549                     {
 550                       cosine = 0.0;
 551                       sine = 1.0;
 552                     }
 553                   else
 554                     {
 555                       cosine = sqrt((v + q) / (2.0 * v));
 556                       sine = p / (2.0 * v * cosine);
 557                     }
 558
 559                   /* apply rotation to A */
 560                   for (i = 0; i < M; i++)
 561                     {
 562                       const double Aik = gsl_matrix_get (A, i, k);
 563                       const double Aij = gsl_matrix_get (A, i, j);
 564                       gsl_matrix_set (A, i, j, Aij * cosine + Aik * sine);
 565                       gsl_matrix_set (A, i, k, -Aij * sine + Aik * cosine);
 566                     }
 567
 568                   gsl_vector_set(S, j, fabs(cosine) * abserr_a + fabs(sine) * abserr_b);
 569                   gsl_vector_set(S, k, fabs(sine) * abserr_a + fabs(cosine) * abserr_b);
 570
 571                   /* apply rotation to Q */
 572                   for (i = 0; i < N; i++)
 573                     {
 574                       const double Qij = gsl_matrix_get (Q, i, j);
 575                       const double Qik = gsl_matrix_get (Q, i, k);
 576                       gsl_matrix_set (Q, i, j, Qij * cosine + Qik * sine);
 577                       gsl_matrix_set (Q, i, k, -Qij * sine + Qik * cosine);
 578                     }
 579                 }
 580             }
 581
 582           /* Sweep completed. */
 583           sweep++;
 584         }
 585
 586       /*
 587        * Orthogonalization complete. Compute singular values.
 588        */
 589
 590       {
 591         double prev_norm = -1.0;
 592
 593         for (j = 0; j < N; j++)
 594           {
 595             gsl_vector_view column = gsl_matrix_column (A, j);
 596             double norm = gsl_blas_dnrm2 (&column.vector);
 597
 598             /* Determine if singular value is zero, according to the
 599                criteria used in the main loop above (i.e. comparison
 600                with norm of previous column). */
 601
 602             if (norm == 0.0 || prev_norm == 0.0
 603                 || (j > 0 && norm <= tolerance * prev_norm))
 604               {
 605                 gsl_vector_set (S, j, 0.0);     /* singular */
 606                 gsl_vector_set_zero (&column.vector);   /* annihilate column */
 607
 608                 prev_norm = 0.0;
 609               }
 610             else
 611               {
 612                 gsl_vector_set (S, j, norm);    /* non-singular */
 613                 gsl_vector_scale (&column.vector, 1.0 / norm);  /* normalize column */
 614
 615                 prev_norm = norm;
 616               }
 617           }
 618       }
 619
 620       if (count > 0)
 621         {
 622           /* reached sweep limit */
 623           GSL_ERROR ("Jacobi iterations did not reach desired tolerance",
 624                      GSL_ETOL);
 625         }
 626
 627       return GSL_SUCCESS;
 628     }
 629 }