1 /* linalg/svd.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007, 2010 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
gsl_linalg_SV_decomp(gsl_matrix * A,gsl_matrix * V,gsl_vector * S,gsl_vector * work)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 int rescale = 0;
164 double scale = 1;
165 double norm = 0;
166
167 /* Find the maximum absolute values of the diagonal and subdiagonal */
168
169 for (i = 0; i < n_block; i++)
170 {
171 double s_i = gsl_vector_get (&S_block.vector, i);
172 double a = fabs(s_i);
173 if (a > norm) norm = a;
174 }
175
176 for (i = 0; i < n_block - 1; i++)
177 {
178 double f_i = gsl_vector_get (&f_block.vector, i);
179 double a = fabs(f_i);
180 if (a > norm) norm = a;
181 }
182
183 /* Temporarily scale the submatrix if necessary */
184
185 if (norm > GSL_SQRT_DBL_MAX)
186 {
187 scale = (norm / GSL_SQRT_DBL_MAX);
188 rescale = 1;
189 }
190 else if (norm < GSL_SQRT_DBL_MIN && norm > 0)
191 {
192 scale = (norm / GSL_SQRT_DBL_MIN);
193 rescale = 1;
194 }
195
196 if (rescale)
197 {
198 gsl_blas_dscal(1.0 / scale, &S_block.vector);
199 gsl_blas_dscal(1.0 / scale, &f_block.vector);
200 }
201
202 /* Perform the implicit QR step */
203
204 qrstep (&S_block.vector, &f_block.vector, &U_block.matrix, &V_block.matrix);
205 /* remove any small off-diagonal elements */
206
207 chop_small_elements (&S_block.vector, &f_block.vector);
208
209 /* Undo the scaling if needed */
210
211 if (rescale)
212 {
213 gsl_blas_dscal(scale, &S_block.vector);
214 gsl_blas_dscal(scale, &f_block.vector);
215 }
216 }
217
218 }
219 }
220
221 /* Make singular values positive by reflections if necessary */
222
223 for (j = 0; j < K; j++)
224 {
225 double Sj = gsl_vector_get (S, j);
226
227 if (Sj < 0.0)
228 {
229 for (i = 0; i < N; i++)
230 {
231 double Vij = gsl_matrix_get (V, i, j);
232 gsl_matrix_set (V, i, j, -Vij);
233 }
234
235 gsl_vector_set (S, j, -Sj);
236 }
237 }
238
239 /* Sort singular values into decreasing order */
240
241 for (i = 0; i < K; i++)
242 {
243 double S_max = gsl_vector_get (S, i);
244 size_t i_max = i;
245
246 for (j = i + 1; j < K; j++)
247 {
248 double Sj = gsl_vector_get (S, j);
249
250 if (Sj > S_max)
251 {
252 S_max = Sj;
253 i_max = j;
254 }
255 }
256
257 if (i_max != i)
258 {
259 /* swap eigenvalues */
260 gsl_vector_swap_elements (S, i, i_max);
261
262 /* swap eigenvectors */
263 gsl_matrix_swap_columns (A, i, i_max);
264 gsl_matrix_swap_columns (V, i, i_max);
265 }
266 }
267
268 return GSL_SUCCESS;
269 }
270
271
272 /* Modified algorithm which is better for M>>N */
273
274 int
gsl_linalg_SV_decomp_mod(gsl_matrix * A,gsl_matrix * X,gsl_matrix * V,gsl_vector * S,gsl_vector * work)275 gsl_linalg_SV_decomp_mod (gsl_matrix * A,
276 gsl_matrix * X,
277 gsl_matrix * V, gsl_vector * S, gsl_vector * work)
278 {
279 size_t i, j;
280
281 const size_t M = A->size1;
282 const size_t N = A->size2;
283
284 if (M < N)
285 {
286 GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
287 }
288 else if (V->size1 != N)
289 {
290 GSL_ERROR ("square matrix V must match second dimension of matrix A",
291 GSL_EBADLEN);
292 }
293 else if (V->size1 != V->size2)
294 {
295 GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
296 }
297 else if (X->size1 != N)
298 {
299 GSL_ERROR ("square matrix X must match second dimension of matrix A",
300 GSL_EBADLEN);
301 }
302 else if (X->size1 != X->size2)
303 {
304 GSL_ERROR ("matrix X must be square", GSL_ENOTSQR);
305 }
306 else if (S->size != N)
307 {
308 GSL_ERROR ("length of vector S must match second dimension of matrix A",
309 GSL_EBADLEN);
310 }
311 else if (work->size != N)
312 {
313 GSL_ERROR ("length of workspace must match second dimension of matrix A",
314 GSL_EBADLEN);
315 }
316
317 if (N == 1)
318 {
319 gsl_vector_view column = gsl_matrix_column (A, 0);
320 double norm = gsl_blas_dnrm2 (&column.vector);
321
322 gsl_vector_set (S, 0, norm);
323 gsl_matrix_set (V, 0, 0, 1.0);
324
325 if (norm != 0.0)
326 {
327 gsl_blas_dscal (1.0/norm, &column.vector);
328 }
329
330 return GSL_SUCCESS;
331 }
332
333 /* Convert A into an upper triangular matrix R */
334
335 for (i = 0; i < N; i++)
336 {
337 gsl_vector_view c = gsl_matrix_column (A, i);
338 gsl_vector_view v = gsl_vector_subvector (&c.vector, i, M - i);
339 double tau_i = gsl_linalg_householder_transform (&v.vector);
340
341 /* Apply the transformation to the remaining columns */
342
343 if (i + 1 < N)
344 {
345 gsl_matrix_view m =
346 gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
347 gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
348 }
349
350 gsl_vector_set (S, i, tau_i);
351 }
352
353 /* Copy the upper triangular part of A into X */
354
355 for (i = 0; i < N; i++)
356 {
357 for (j = 0; j < i; j++)
358 {
359 gsl_matrix_set (X, i, j, 0.0);
360 }
361
362 {
363 double Aii = gsl_matrix_get (A, i, i);
364 gsl_matrix_set (X, i, i, Aii);
365 }
366
367 for (j = i + 1; j < N; j++)
368 {
369 double Aij = gsl_matrix_get (A, i, j);
370 gsl_matrix_set (X, i, j, Aij);
371 }
372 }
373
374 /* Convert A into an orthogonal matrix L */
375
376 for (j = N; j-- > 0;)
377 {
378 /* Householder column transformation to accumulate L */
379 double tj = gsl_vector_get (S, j);
380 gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M - j, N - j);
381 gsl_linalg_householder_hm1 (tj, &m.matrix);
382 }
383
384 /* unpack R into X V S */
385
386 gsl_linalg_SV_decomp (X, V, S, work);
387
388 /* Multiply L by X, to obtain U = L X, stored in U */
389
390 {
391 gsl_vector_view sum = gsl_vector_subvector (work, 0, N);
392
393 for (i = 0; i < M; i++)
394 {
395 gsl_vector_view L_i = gsl_matrix_row (A, i);
396 gsl_vector_set_zero (&sum.vector);
397
398 for (j = 0; j < N; j++)
399 {
400 double Lij = gsl_vector_get (&L_i.vector, j);
401 gsl_vector_view X_j = gsl_matrix_row (X, j);
402 gsl_blas_daxpy (Lij, &X_j.vector, &sum.vector);
403 }
404
405 gsl_vector_memcpy (&L_i.vector, &sum.vector);
406 }
407 }
408
409 return GSL_SUCCESS;
410 }
411
412
413 /* Solves the system A x = b using the SVD factorization
414 *
415 * A = U S V^T
416 *
417 * to obtain x. For M x N systems it finds the solution in the least
418 * squares sense.
419 */
420
421 int
gsl_linalg_SV_solve(const gsl_matrix * U,const gsl_matrix * V,const gsl_vector * S,const gsl_vector * b,gsl_vector * x)422 gsl_linalg_SV_solve (const gsl_matrix * U,
423 const gsl_matrix * V,
424 const gsl_vector * S,
425 const gsl_vector * b, gsl_vector * x)
426 {
427 if (U->size1 != b->size)
428 {
429 GSL_ERROR ("first dimension of matrix U must size of vector b",
430 GSL_EBADLEN);
431 }
432 else if (U->size2 != S->size)
433 {
434 GSL_ERROR ("length of vector S must match second dimension of matrix U",
435 GSL_EBADLEN);
436 }
437 else if (V->size1 != V->size2)
438 {
439 GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
440 }
441 else if (S->size != V->size1)
442 {
443 GSL_ERROR ("length of vector S must match size of matrix V",
444 GSL_EBADLEN);
445 }
446 else if (V->size2 != x->size)
447 {
448 GSL_ERROR ("size of matrix V must match size of vector x", GSL_EBADLEN);
449 }
450 else
451 {
452 const size_t N = U->size2;
453 size_t i;
454
455 gsl_vector *w = gsl_vector_calloc (N);
456
457 gsl_blas_dgemv (CblasTrans, 1.0, U, b, 0.0, w);
458
459 for (i = 0; i < N; i++)
460 {
461 double wi = gsl_vector_get (w, i);
462 double alpha = gsl_vector_get (S, i);
463 if (alpha != 0)
464 alpha = 1.0 / alpha;
465 gsl_vector_set (w, i, alpha * wi);
466 }
467
468 gsl_blas_dgemv (CblasNoTrans, 1.0, V, w, 0.0, x);
469
470 gsl_vector_free (w);
471
472 return GSL_SUCCESS;
473 }
474 }
475
476 /*
477 gsl_linalg_SV_leverage()
478 Compute statistical leverage values of a matrix:
479
480 h = diag(A (A^T A)^{-1} A^T)
481
482 Inputs: U - U matrix in SVD decomposition
483 h - (output) vector of leverages
484
485 Return: success or error
486 */
487
488 int
gsl_linalg_SV_leverage(const gsl_matrix * U,gsl_vector * h)489 gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h)
490 {
491 const size_t M = U->size1;
492
493 if (M != h->size)
494 {
495 GSL_ERROR ("first dimension of matrix U must match size of vector h",
496 GSL_EBADLEN);
497 }
498 else
499 {
500 size_t i;
501
502 for (i = 0; i < M; ++i)
503 {
504 gsl_vector_const_view v = gsl_matrix_const_row(U, i);
505 double hi;
506
507 gsl_blas_ddot(&v.vector, &v.vector, &hi);
508 gsl_vector_set(h, i, hi);
509 }
510
511 return GSL_SUCCESS;
512 }
513 } /* gsl_linalg_SV_leverage() */
514
515 /* This is a the jacobi version */
516 /* Author: G. Jungman */
517
518 /*
519 * Algorithm due to J.C. Nash, Compact Numerical Methods for
520 * Computers (New York: Wiley and Sons, 1979), chapter 3.
521 * See also Algorithm 4.1 in
522 * James Demmel, Kresimir Veselic, "Jacobi's Method is more
523 * accurate than QR", Lapack Working Note 15 (LAWN15), October 1989.
524 * Available from netlib.
525 *
526 * Based on code by Arthur Kosowsky, Rutgers University
527 * kosowsky@physics.rutgers.edu
528 *
529 * Another relevant paper is, P.P.M. De Rijk, "A One-Sided Jacobi
530 * Algorithm for computing the singular value decomposition on a
531 * vector computer", SIAM Journal of Scientific and Statistical
532 * Computing, Vol 10, No 2, pp 359-371, March 1989.
533 *
534 */
535
536 int
gsl_linalg_SV_decomp_jacobi(gsl_matrix * A,gsl_matrix * Q,gsl_vector * S)537 gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * Q, gsl_vector * S)
538 {
539 if (A->size1 < A->size2)
540 {
541 /* FIXME: only implemented M>=N case so far */
542
543 GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
544 }
545 else if (Q->size1 != A->size2)
546 {
547 GSL_ERROR ("square matrix Q must match second dimension of matrix A",
548 GSL_EBADLEN);
549 }
550 else if (Q->size1 != Q->size2)
551 {
552 GSL_ERROR ("matrix Q must be square", GSL_ENOTSQR);
553 }
554 else if (S->size != A->size2)
555 {
556 GSL_ERROR ("length of vector S must match second dimension of matrix A",
557 GSL_EBADLEN);
558 }
559 else
560 {
561 const size_t M = A->size1;
562 const size_t N = A->size2;
563 size_t i, j, k;
564
565 /* Initialize the rotation counter and the sweep counter. */
566 int count = 1;
567 int sweep = 0;
568 int sweepmax = 5*N;
569
570 double tolerance = 10 * M * GSL_DBL_EPSILON;
571
572 /* Always do at least 12 sweeps. */
573 sweepmax = GSL_MAX (sweepmax, 12);
574
575 /* Set Q to the identity matrix. */
576 gsl_matrix_set_identity (Q);
577
578 /* Store the column error estimates in S, for use during the
579 orthogonalization */
580
581 for (j = 0; j < N; j++)
582 {
583 gsl_vector_view cj = gsl_matrix_column (A, j);
584 double sj = gsl_blas_dnrm2 (&cj.vector);
585 gsl_vector_set(S, j, GSL_DBL_EPSILON * sj);
586 }
587
588 /* Orthogonalize A by plane rotations. */
589
590 while (count > 0 && sweep <= sweepmax)
591 {
592 /* Initialize rotation counter. */
593 count = N * (N - 1) / 2;
594
595 for (j = 0; j < N - 1; j++)
596 {
597 for (k = j + 1; k < N; k++)
598 {
599 double a = 0.0;
600 double b = 0.0;
601 double p = 0.0;
602 double q = 0.0;
603 double cosine, sine;
604 double v;
605 double abserr_a, abserr_b;
606 int sorted, orthog, noisya, noisyb;
607
608 gsl_vector_view cj = gsl_matrix_column (A, j);
609 gsl_vector_view ck = gsl_matrix_column (A, k);
610
611 gsl_blas_ddot (&cj.vector, &ck.vector, &p);
612 p *= 2.0 ; /* equation 9a: p = 2 x.y */
613
614 a = gsl_blas_dnrm2 (&cj.vector);
615 b = gsl_blas_dnrm2 (&ck.vector);
616
617 q = a * a - b * b;
618 v = hypot(p, q);
619
620 /* test for columns j,k orthogonal, or dominant errors */
621
622 abserr_a = gsl_vector_get(S,j);
623 abserr_b = gsl_vector_get(S,k);
624
625 sorted = (GSL_COERCE_DBL(a) >= GSL_COERCE_DBL(b));
626 orthog = (fabs (p) <= tolerance * GSL_COERCE_DBL(a * b));
627 noisya = (a < abserr_a);
628 noisyb = (b < abserr_b);
629
630 if (sorted && (orthog || noisya || noisyb))
631 {
632 count--;
633 continue;
634 }
635
636 /* calculate rotation angles */
637 if (v == 0 || !sorted)
638 {
639 cosine = 0.0;
640 sine = 1.0;
641 }
642 else
643 {
644 cosine = sqrt((v + q) / (2.0 * v));
645 sine = p / (2.0 * v * cosine);
646 }
647
648 /* apply rotation to A */
649 for (i = 0; i < M; i++)
650 {
651 const double Aik = gsl_matrix_get (A, i, k);
652 const double Aij = gsl_matrix_get (A, i, j);
653 gsl_matrix_set (A, i, j, Aij * cosine + Aik * sine);
654 gsl_matrix_set (A, i, k, -Aij * sine + Aik * cosine);
655 }
656
657 gsl_vector_set(S, j, fabs(cosine) * abserr_a + fabs(sine) * abserr_b);
658 gsl_vector_set(S, k, fabs(sine) * abserr_a + fabs(cosine) * abserr_b);
659
660 /* apply rotation to Q */
661 for (i = 0; i < N; i++)
662 {
663 const double Qij = gsl_matrix_get (Q, i, j);
664 const double Qik = gsl_matrix_get (Q, i, k);
665 gsl_matrix_set (Q, i, j, Qij * cosine + Qik * sine);
666 gsl_matrix_set (Q, i, k, -Qij * sine + Qik * cosine);
667 }
668 }
669 }
670
671 /* Sweep completed. */
672 sweep++;
673 }
674
675 /*
676 * Orthogonalization complete. Compute singular values.
677 */
678
679 {
680 double prev_norm = -1.0;
681
682 for (j = 0; j < N; j++)
683 {
684 gsl_vector_view column = gsl_matrix_column (A, j);
685 double norm = gsl_blas_dnrm2 (&column.vector);
686
687 /* Determine if singular value is zero, according to the
688 criteria used in the main loop above (i.e. comparison
689 with norm of previous column). */
690
691 if (norm == 0.0 || prev_norm == 0.0
692 || (j > 0 && norm <= tolerance * prev_norm))
693 {
694 gsl_vector_set (S, j, 0.0); /* singular */
695 gsl_vector_set_zero (&column.vector); /* annihilate column */
696
697 prev_norm = 0.0;
698 }
699 else
700 {
701 gsl_vector_set (S, j, norm); /* non-singular */
702 gsl_vector_scale (&column.vector, 1.0 / norm); /* normalize column */
703
704 prev_norm = norm;
705 }
706 }
707 }
708
709 if (count > 0)
710 {
711 /* reached sweep limit */
712 GSL_ERROR ("Jacobi iterations did not reach desired tolerance",
713 GSL_ETOL);
714 }
715
716 return GSL_SUCCESS;
717 }
718 }
719