1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Neil Firth
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18
19 Adapted from the TNT project
20 http://math.nist.gov/tnt/download.html
21
22 This software was developed at the National Institute of Standards
23 and Technology (NIST) by employees of the Federal Government in the
24 course of their official duties. Pursuant to title 17 Section 105
25 of the United States Code this software is not subject to copyright
26 protection and is in the public domain. NIST assumes no responsibility
27 whatsoever for its use by other parties, and makes no guarantees,
28 expressed or implied, about its quality, reliability, or any other
29 characteristic.
30
31 We would appreciate acknowledgement if the software is incorporated in
32 redistributable libraries or applications.
33*/
34
35
36#include <ql/math/matrixutilities/svd.hpp>
37
38namespace QuantLib {
39
40 namespace {
41
42 /* returns hypotenuse of real (non-complex) scalars a and b by
43 avoiding underflow/overflow
44 using (a * sqrt( 1 + (b/a) * (b/a))), rather than
45 sqrt(a*a + b*b).
46 */
47 Real hypot(const Real &a, const Real &b) {
48 if (a == 0) {
49 return std::fabs(x: b);
50 } else {
51 Real c = b/a;
52 return std::fabs(x: a) * std::sqrt(x: 1 + c*c);
53 }
54 }
55
56 }
57
58
59 SVD::SVD(const Matrix& M) {
60
61 using std::swap;
62
63 Matrix A;
64
65 /* The implementation requires that rows > columns.
66 If this is not the case, we decompose M^T instead.
67 Swapping the resulting U and V gives the desired
68 result for M as
69
70 M^T = U S V^T (decomposition of M^T)
71
72 M = (U S V^T)^T (transpose)
73
74 M = (V^T^T S^T U^T) ((AB)^T = B^T A^T)
75
76 M = V S U^T (idempotence of transposition,
77 symmetry of diagonal matrix S)
78
79 */
80
81 if (M.rows() >= M.columns()) {
82 A = M;
83 transpose_ = false;
84 } else {
85 A = transpose(m: M);
86 transpose_ = true;
87 }
88
89 m_ = A.rows();
90 n_ = A.columns();
91
92 // we're sure that m_ >= n_
93
94 s_ = Array(n_);
95 U_ = Matrix(m_,n_, 0.0);
96 V_ = Matrix(n_,n_);
97 Array e(n_);
98 Array work(m_);
99 Integer i, j, k;
100
101 // Reduce A to bidiagonal form, storing the diagonal elements
102 // in s and the super-diagonal elements in e.
103
104 Integer nct = std::min(a: m_-1,b: n_);
105 Integer nrt = std::max(a: 0,b: n_-2);
106 for (k = 0; k < std::max(a: nct,b: nrt); k++) {
107 if (k < nct) {
108
109 // Compute the transformation for the k-th column and
110 // place the k-th diagonal in s[k].
111 // Compute 2-norm of k-th column without under/overflow.
112 s_[k] = 0;
113 for (i = k; i < m_; i++) {
114 s_[k] = hypot(a: s_[k],b: A[i][k]);
115 }
116 if (s_[k] != 0.0) {
117 if (A[k][k] < 0.0) {
118 s_[k] = -s_[k];
119 }
120 for (i = k; i < m_; i++) {
121 A[i][k] /= s_[k];
122 }
123 A[k][k] += 1.0;
124 }
125 s_[k] = -s_[k];
126 }
127 for (j = k+1; j < n_; j++) {
128 if ((k < nct) && (s_[k] != 0.0)) {
129
130 // Apply the transformation.
131
132 Real t = 0;
133 for (i = k; i < m_; i++) {
134 t += A[i][k]*A[i][j];
135 }
136 t = -t/A[k][k];
137 for (i = k; i < m_; i++) {
138 A[i][j] += t*A[i][k];
139 }
140 }
141
142 // Place the k-th row of A into e for the
143 // subsequent calculation of the row transformation.
144
145 e[j] = A[k][j];
146 }
147 if (k < nct) {
148
149 // Place the transformation in U for subsequent back
150 // multiplication.
151
152 for (i = k; i < m_; i++) {
153 U_[i][k] = A[i][k];
154 }
155 }
156 if (k < nrt) {
157
158 // Compute the k-th row transformation and place the
159 // k-th super-diagonal in e[k].
160 // Compute 2-norm without under/overflow.
161 e[k] = 0;
162 for (i = k+1; i < n_; i++) {
163 e[k] = hypot(a: e[k],b: e[i]);
164 }
165 if (e[k] != 0.0) {
166 if (e[k+1] < 0.0) {
167 e[k] = -e[k];
168 }
169 for (i = k+1; i < n_; i++) {
170 e[i] /= e[k];
171 }
172 e[k+1] += 1.0;
173 }
174 e[k] = -e[k];
175 if ((k+1 < m_) && (e[k] != 0.0)) {
176
177 // Apply the transformation.
178
179 for (i = k+1; i < m_; i++) {
180 work[i] = 0.0;
181 }
182 for (j = k+1; j < n_; j++) {
183 for (i = k+1; i < m_; i++) {
184 work[i] += e[j]*A[i][j];
185 }
186 }
187 for (j = k+1; j < n_; j++) {
188 Real t = -e[j]/e[k+1];
189 for (i = k+1; i < m_; i++) {
190 A[i][j] += t*work[i];
191 }
192 }
193 }
194
195 // Place the transformation in V for subsequent
196 // back multiplication.
197
198 for (i = k+1; i < n_; i++) {
199 V_[i][k] = e[i];
200 }
201 }
202 }
203
204 // Set up the final bidiagonal matrix or order n.
205
206 if (nct < n_) {
207 s_[nct] = A[nct][nct];
208 }
209 if (nrt+1 < n_) {
210 e[nrt] = A[nrt][n_-1];
211 }
212 e[n_-1] = 0.0;
213
214 // generate U
215
216 for (j = nct; j < n_; j++) {
217 for (i = 0; i < m_; i++) {
218 U_[i][j] = 0.0;
219 }
220 U_[j][j] = 1.0;
221 }
222 for (k = nct-1; k >= 0; --k) {
223 if (s_[k] != 0.0) {
224 for (j = k+1; j < n_; ++j) {
225 Real t = 0;
226 for (i = k; i < m_; i++) {
227 t += U_[i][k]*U_[i][j];
228 }
229 t = -t/U_[k][k];
230 for (i = k; i < m_; i++) {
231 U_[i][j] += t*U_[i][k];
232 }
233 }
234 for (i = k; i < m_; i++ ) {
235 U_[i][k] = -U_[i][k];
236 }
237 U_[k][k] = 1.0 + U_[k][k];
238 for (i = 0; i < k-1; i++) {
239 U_[i][k] = 0.0;
240 }
241 } else {
242 for (i = 0; i < m_; i++) {
243 U_[i][k] = 0.0;
244 }
245 U_[k][k] = 1.0;
246 }
247 }
248
249 // generate V
250
251 for (k = n_-1; k >= 0; --k) {
252 if ((k < nrt) && (e[k] != 0.0)) {
253 for (j = k+1; j < n_; ++j) {
254 Real t = 0;
255 for (i = k+1; i < n_; i++) {
256 t += V_[i][k]*V_[i][j];
257 }
258 t = -t/V_[k+1][k];
259 for (i = k+1; i < n_; i++) {
260 V_[i][j] += t*V_[i][k];
261 }
262 }
263 }
264 for (i = 0; i < n_; i++) {
265 V_[i][k] = 0.0;
266 }
267 V_[k][k] = 1.0;
268 }
269
270 // Main iteration loop for the singular values.
271
272 Integer p = n_, pp = p-1;
273 Integer iter = 0;
274 Real eps = std::pow(x: 2.0,y: -52.0);
275 while (p > 0) {
276 Integer k;
277 Integer kase;
278
279 // Here is where a test for too many iterations would go.
280
281 // This section of the program inspects for
282 // negligible elements in the s and e arrays. On
283 // completion the variables kase and k are set as follows.
284
285 // kase = 1 if s(p) and e[k-1] are negligible and k<p
286 // kase = 2 if s(k) is negligible and k<p
287 // kase = 3 if e[k-1] is negligible, k<p, and
288 // s(k), ..., s(p) are not negligible (qr step).
289 // kase = 4 if e(p-1) is negligible (convergence).
290
291 for (k = p-2; k >= -1; --k) {
292 if (k == -1) {
293 break;
294 }
295 if (std::fabs(x: e[k]) <= eps*(std::fabs(x: s_[k]) +
296 std::fabs(x: s_[k+1]))) {
297 e[k] = 0.0;
298 break;
299 }
300 }
301 if (k == p-2) {
302 kase = 4;
303 } else {
304 Integer ks;
305 for (ks = p-1; ks >= k; --ks) {
306 if (ks == k) {
307 break;
308 }
309 Real t = (ks != p ? Real(std::fabs(x: e[ks])) : 0.) +
310 (ks != k+1 ? Real(std::fabs(x: e[ks-1])) : 0.);
311 if (std::fabs(x: s_[ks]) <= eps*t) {
312 s_[ks] = 0.0;
313 break;
314 }
315 }
316 if (ks == k) {
317 kase = 3;
318 } else if (ks == p-1) {
319 kase = 1;
320 } else {
321 kase = 2;
322 k = ks;
323 }
324 }
325 k++;
326
327 // Perform the task indicated by kase.
328
329 switch (kase) { // NOLINT(bugprone-switch-missing-default-case)
330
331 // Deflate negligible s(p).
332
333 case 1: {
334 Real f = e[p-2];
335 e[p-2] = 0.0;
336 for (j = p-2; j >= k; --j) {
337 Real t = hypot(a: s_[j],b: f);
338 Real cs = s_[j]/t;
339 Real sn = f/t;
340 s_[j] = t;
341 if (j != k) {
342 f = -sn*e[j-1];
343 e[j-1] = cs*e[j-1];
344 }
345 for (i = 0; i < n_; i++) {
346 t = cs*V_[i][j] + sn*V_[i][p-1];
347 V_[i][p-1] = -sn*V_[i][j] + cs*V_[i][p-1];
348 V_[i][j] = t;
349 }
350 }
351 }
352 break;
353
354 // Split at negligible s(k).
355
356 case 2: {
357 Real f = e[k-1];
358 e[k-1] = 0.0;
359 for (j = k; j < p; j++) {
360 Real t = hypot(a: s_[j],b: f);
361 Real cs = s_[j]/t;
362 Real sn = f/t;
363 s_[j] = t;
364 f = -sn*e[j];
365 e[j] = cs*e[j];
366 for (i = 0; i < m_; i++) {
367 t = cs*U_[i][j] + sn*U_[i][k-1];
368 U_[i][k-1] = -sn*U_[i][j] + cs*U_[i][k-1];
369 U_[i][j] = t;
370 }
371 }
372 }
373 break;
374
375 // Perform one qr step.
376
377 case 3: {
378
379 // Calculate the shift.
380 Real scale = std::max(
381 a: std::max(
382 a: std::max(
383 a: std::max(a: std::fabs(x: s_[p-1]),
384 b: std::fabs(x: s_[p-2])),
385 b: std::fabs(x: e[p-2])),
386 b: std::fabs(x: s_[k])),
387 b: std::fabs(x: e[k]));
388 Real sp = s_[p-1]/scale;
389 Real spm1 = s_[p-2]/scale;
390 Real epm1 = e[p-2]/scale;
391 Real sk = s_[k]/scale;
392 Real ek = e[k]/scale;
393 Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
394 Real c = (sp*epm1)*(sp*epm1);
395 Real shift = 0.0;
396 if ((b != 0.0) || (c != 0.0)) {
397 shift = std::sqrt(x: b*b + c);
398 if (b < 0.0) {
399 shift = -shift;
400 }
401 shift = c/(b + shift);
402 }
403 Real f = (sk + sp)*(sk - sp) + shift;
404 Real g = sk*ek;
405
406 // Chase zeros.
407
408 for (j = k; j < p-1; j++) {
409 Real t = hypot(a: f,b: g);
410 Real cs = f/t;
411 Real sn = g/t;
412 if (j != k) {
413 e[j-1] = t;
414 }
415 f = cs*s_[j] + sn*e[j];
416 e[j] = cs*e[j] - sn*s_[j];
417 g = sn*s_[j+1];
418 s_[j+1] = cs*s_[j+1];
419 for (i = 0; i < n_; i++) {
420 t = cs*V_[i][j] + sn*V_[i][j+1];
421 V_[i][j+1] = -sn*V_[i][j] + cs*V_[i][j+1];
422 V_[i][j] = t;
423 }
424 t = hypot(a: f,b: g);
425 cs = f/t;
426 sn = g/t;
427 s_[j] = t;
428 f = cs*e[j] + sn*s_[j+1];
429 s_[j+1] = -sn*e[j] + cs*s_[j+1];
430 g = sn*e[j+1];
431 e[j+1] = cs*e[j+1];
432 if (j < m_-1) {
433 for (i = 0; i < m_; i++) {
434 t = cs*U_[i][j] + sn*U_[i][j+1];
435 U_[i][j+1] = -sn*U_[i][j] + cs*U_[i][j+1];
436 U_[i][j] = t;
437 }
438 }
439 }
440 e[p-2] = f;
441 iter = iter + 1;
442 }
443 break;
444
445 // Convergence.
446
447 case 4: {
448
449 // Make the singular values positive.
450
451 if (s_[k] <= 0.0) {
452 s_[k] = (s_[k] < 0.0 ? Real(-s_[k]) : 0.0);
453 for (i = 0; i <= pp; i++) {
454 V_[i][k] = -V_[i][k];
455 }
456 }
457
458 // Order the singular values.
459
460 while (k < pp) {
461 if (s_[k] >= s_[k+1]) {
462 break;
463 }
464 swap(a&: s_[k], b&: s_[k+1]);
465 if (k < n_-1) {
466 for (i = 0; i < n_; i++) {
467 swap(a&: V_[i][k], b&: V_[i][k+1]);
468 }
469 }
470 if (k < m_-1) {
471 for (i = 0; i < m_; i++) {
472 swap(a&: U_[i][k], b&: U_[i][k+1]);
473 }
474 }
475 k++;
476 }
477 iter = 0;
478 --p;
479 }
480 break;
481 }
482 }
483 }
484
485 const Matrix& SVD::U() const {
486 return (transpose_ ? V_ : U_);
487 }
488
489 const Matrix& SVD::V() const {
490 return (transpose_ ? U_ : V_);
491 }
492
493 const Array& SVD::singularValues() const {
494 return s_;
495 }
496
497 Matrix SVD::S() const {
498 Matrix S(n_,n_);
499 for (Size i = 0; i < Size(n_); i++) {
500 for (Size j = 0; j < Size(n_); j++) {
501 S[i][j] = 0.0;
502 }
503 S[i][i] = s_[i];
504 }
505 return S;
506 }
507
508 Real SVD::norm2() const {
509 return s_[0];
510 }
511
512 Real SVD::cond() const {
513 return s_[0]/s_[n_-1];
514 }
515
516 Size SVD::rank() const {
517 constexpr auto eps = QL_EPSILON;
518 Real tol = m_*s_[0]*eps;
519 Size r = 0;
520 for (Real i : s_) {
521 if (i > tol) {
522 r++;
523 }
524 }
525 return r;
526 }
527
528 Array SVD::solveFor(const Array& b) const{
529 Matrix W(n_, n_, 0.0);
530 const Size numericalRank = this->rank();
531 for (Size i=0; i<numericalRank; i++)
532 W[i][i] = 1./s_[i];
533
534 Matrix inverse = V()* W * transpose(m: U());
535 return inverse * b;
536 }
537
538}
539
540

source code of quantlib/ql/math/matrixutilities/svd.cpp