1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003, 2004, 2007 Ferdinando Ametrano
5 Copyright (C) 2006 Yiping Chen
6 Copyright (C) 2007 Neil Firth
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include <ql/math/comparison.hpp>
23#include <ql/math/matrixutilities/choleskydecomposition.hpp>
24#include <ql/math/matrixutilities/pseudosqrt.hpp>
25#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
26#include <ql/math/optimization/conjugategradient.hpp>
27#include <ql/math/optimization/constraint.hpp>
28#include <ql/math/optimization/problem.hpp>
29#include <utility>
30
31namespace QuantLib {
32
33 namespace {
34
35 #if defined(QL_EXTRA_SAFETY_CHECKS)
36 void checkSymmetry(const Matrix& matrix) {
37 Size size = matrix.rows();
38 QL_REQUIRE(size == matrix.columns(),
39 "non square matrix: " << size << " rows, " <<
40 matrix.columns() << " columns");
41 for (Size i=0; i<size; ++i)
42 for (Size j=0; j<i; ++j)
43 QL_REQUIRE(close(matrix[i][j], matrix[j][i]),
44 "non symmetric matrix: " <<
45 "[" << i << "][" << j << "]=" << matrix[i][j] <<
46 ", [" << j << "][" << i << "]=" << matrix[j][i]);
47 }
48 #endif
49
50 void normalizePseudoRoot(const Matrix& matrix,
51 Matrix& pseudo) {
52 Size size = matrix.rows();
53 QL_REQUIRE(size == pseudo.rows(),
54 "matrix/pseudo mismatch: matrix rows are " << size <<
55 " while pseudo rows are " << pseudo.columns());
56 Size pseudoCols = pseudo.columns();
57
58 // row normalization
59 for (Size i=0; i<size; ++i) {
60 Real norm = 0.0;
61 for (Size j=0; j<pseudoCols; ++j)
62 norm += pseudo[i][j]*pseudo[i][j];
63 if (norm>0.0) {
64 Real normAdj = std::sqrt(x: matrix[i][i]/norm);
65 for (Size j=0; j<pseudoCols; ++j)
66 pseudo[i][j] *= normAdj;
67 }
68 }
69
70
71 }
72
73 //cost function for hypersphere and lower-diagonal algorithm
74 class HypersphereCostFunction : public CostFunction {
75 private:
76 Size size_;
77 bool lowerDiagonal_;
78 Matrix targetMatrix_;
79 Array targetVariance_;
80 mutable Matrix currentRoot_, tempMatrix_, currentMatrix_;
81 public:
82 HypersphereCostFunction(const Matrix& targetMatrix,
83 Array targetVariance,
84 bool lowerDiagonal)
85 : size_(targetMatrix.rows()), lowerDiagonal_(lowerDiagonal),
86 targetMatrix_(targetMatrix), targetVariance_(std::move(targetVariance)),
87 currentRoot_(size_, size_), tempMatrix_(size_, size_), currentMatrix_(size_, size_) {}
88 Array values(const Array&) const override {
89 QL_FAIL("values method not implemented");
90 }
91 Real value(const Array& x) const override {
92 Size i,j,k;
93 std::fill(first: currentRoot_.begin(), last: currentRoot_.end(), value: 1.0);
94 if (lowerDiagonal_) {
95 for (i=0; i<size_; i++) {
96 for (k=0; k<size_; k++) {
97 if (k>i) {
98 currentRoot_[i][k]=0;
99 } else {
100 for (j=0; j<=k; j++) {
101 if (j == k && k!=i)
102 currentRoot_[i][k] *=
103 std::cos(x: x[i*(i-1)/2+j]);
104 else if (j!=i)
105 currentRoot_[i][k] *=
106 std::sin(x: x[i*(i-1)/2+j]);
107 }
108 }
109 }
110 }
111 } else {
112 for (i=0; i<size_; i++) {
113 for (k=0; k<size_; k++) {
114 for (j=0; j<=k; j++) {
115 if (j == k && k!=size_-1)
116 currentRoot_[i][k] *=
117 std::cos(x: x[j*size_+i]);
118 else if (j!=size_-1)
119 currentRoot_[i][k] *=
120 std::sin(x: x[j*size_+i]);
121 }
122 }
123 }
124 }
125 Real temp, error=0;
126 tempMatrix_ = transpose(m: currentRoot_);
127 currentMatrix_ = currentRoot_ * tempMatrix_;
128 for (i=0;i<size_;i++) {
129 for (j=0;j<size_;j++) {
130 temp = currentMatrix_[i][j]*targetVariance_[i]
131 *targetVariance_[j]-targetMatrix_[i][j];
132 error += temp*temp;
133 }
134 }
135 return error;
136 }
137 };
138
139 // Optimization function for hypersphere and lower-diagonal algorithm
140 Matrix hypersphereOptimize(const Matrix& targetMatrix,
141 const Matrix& currentRoot,
142 const bool lowerDiagonal) {
143 Size i,j,k,size = targetMatrix.rows();
144 Matrix result(currentRoot);
145 Array variance(size, 0);
146 for (i=0; i<size; i++){
147 variance[i]=std::sqrt(x: targetMatrix[i][i]);
148 }
149 if (lowerDiagonal) {
150 Matrix approxMatrix(result*transpose(m: result));
151 result = CholeskyDecomposition(m: approxMatrix, flexible: true);
152 for (i=0; i<size; i++) {
153 for (j=0; j<size; j++) {
154 result[i][j]/=std::sqrt(x: approxMatrix[i][i]);
155 }
156 }
157 } else {
158 for (i=0; i<size; i++) {
159 for (j=0; j<size; j++) {
160 result[i][j]/=variance[i];
161 }
162 }
163 }
164
165 ConjugateGradient optimize;
166 EndCriteria endCriteria(100, 10, 1e-8, 1e-8, 1e-8);
167 HypersphereCostFunction costFunction(targetMatrix, variance,
168 lowerDiagonal);
169 NoConstraint constraint;
170
171 // hypersphere vector optimization
172
173 if (lowerDiagonal) {
174 Array theta(size * (size-1)/2);
175 const Real eps=1e-16;
176 for (i=1; i<size; i++) {
177 for (j=0; j<i; j++) {
178 theta[i*(i-1)/2+j]=result[i][j];
179 if (theta[i*(i-1)/2+j]>1-eps)
180 theta[i*(i-1)/2+j]=1-eps;
181 if (theta[i*(i-1)/2+j]<-1+eps)
182 theta[i*(i-1)/2+j]=-1+eps;
183 for (k=0; k<j; k++) {
184 theta[i*(i-1)/2+j] /= std::sin(x: theta[i*(i-1)/2+k]);
185 if (theta[i*(i-1)/2+j]>1-eps)
186 theta[i*(i-1)/2+j]=1-eps;
187 if (theta[i*(i-1)/2+j]<-1+eps)
188 theta[i*(i-1)/2+j]=-1+eps;
189 }
190 theta[i*(i-1)/2+j] = std::acos(x: theta[i*(i-1)/2+j]);
191 if (j==i-1) {
192 if (result[i][i]<0)
193 theta[i*(i-1)/2+j]=-theta[i*(i-1)/2+j];
194 }
195 }
196 }
197 Problem p(costFunction, constraint, theta);
198 optimize.minimize(P&: p, endCriteria);
199 theta = p.currentValue();
200 std::fill(first: result.begin(),last: result.end(),value: 1.0);
201 for (i=0; i<size; i++) {
202 for (k=0; k<size; k++) {
203 if (k>i) {
204 result[i][k]=0;
205 } else {
206 for (j=0; j<=k; j++) {
207 if (j == k && k!=i)
208 result[i][k] *=
209 std::cos(x: theta[i*(i-1)/2+j]);
210 else if (j!=i)
211 result[i][k] *=
212 std::sin(x: theta[i*(i-1)/2+j]);
213 }
214 }
215 }
216 }
217 } else {
218 Array theta(size * (size-1));
219 const Real eps=1e-16;
220 for (i=0; i<size; i++) {
221 for (j=0; j<size-1; j++) {
222 theta[j*size+i]=result[i][j];
223 if (theta[j*size+i]>1-eps)
224 theta[j*size+i]=1-eps;
225 if (theta[j*size+i]<-1+eps)
226 theta[j*size+i]=-1+eps;
227 for (k=0;k<j;k++) {
228 theta[j*size+i] /= std::sin(x: theta[k*size+i]);
229 if (theta[j*size+i]>1-eps)
230 theta[j*size+i]=1-eps;
231 if (theta[j*size+i]<-1+eps)
232 theta[j*size+i]=-1+eps;
233 }
234 theta[j*size+i] = std::acos(x: theta[j*size+i]);
235 if (j==size-2) {
236 if (result[i][j+1]<0)
237 theta[j*size+i]=-theta[j*size+i];
238 }
239 }
240 }
241 Problem p(costFunction, constraint, theta);
242 optimize.minimize(P&: p, endCriteria);
243 theta=p.currentValue();
244 std::fill(first: result.begin(),last: result.end(),value: 1.0);
245 for (i=0; i<size; i++) {
246 for (k=0; k<size; k++) {
247 for (j=0; j<=k; j++) {
248 if (j == k && k!=size-1)
249 result[i][k] *= std::cos(x: theta[j*size+i]);
250 else if (j!=size-1)
251 result[i][k] *= std::sin(x: theta[j*size+i]);
252 }
253 }
254 }
255 }
256
257 for (i=0; i<size; i++) {
258 for (j=0; j<size; j++) {
259 result[i][j]*=variance[i];
260 }
261 }
262 return result;
263 }
264
265 // Matrix infinity norm. See Golub and van Loan (2.3.10) or
266 // <http://en.wikipedia.org/wiki/Matrix_norm>
267 Real normInf(const Matrix& M) {
268 Size rows = M.rows();
269 Size cols = M.columns();
270 Real norm = 0.0;
271 for (Size i=0; i<rows; ++i) {
272 Real colSum = 0.0;
273 for (Size j=0; j<cols; ++j)
274 colSum += std::fabs(x: M[i][j]);
275 norm = std::max(a: norm, b: colSum);
276 }
277 return norm;
278 }
279
280 // Take a matrix and make all the diagonal entries 1.
281 Matrix projectToUnitDiagonalMatrix(const Matrix& M) {
282 Size size = M.rows();
283 QL_REQUIRE(size == M.columns(),
284 "matrix not square");
285
286 Matrix result(M);
287 for (Size i=0; i<size; ++i)
288 result[i][i] = 1.0;
289
290 return result;
291 }
292
293 // Take a matrix and make all the eigenvalues non-negative
294 Matrix projectToPositiveSemidefiniteMatrix(Matrix& M) {
295 Size size = M.rows();
296 QL_REQUIRE(size == M.columns(),
297 "matrix not square");
298
299 Matrix diagonal(size, size, 0.0);
300 SymmetricSchurDecomposition jd(M);
301 for (Size i=0; i<size; ++i)
302 diagonal[i][i] = std::max<Real>(a: jd.eigenvalues()[i], b: 0.0);
303
304 Matrix result =
305 jd.eigenvectors()*diagonal*transpose(m: jd.eigenvectors());
306 return result;
307 }
308
309 // implementation of the Higham algorithm to find the nearest
310 // correlation matrix.
311 Matrix highamImplementation(const Matrix& A, const Size maxIterations, const Real& tolerance) {
312
313 Size size = A.rows();
314 Matrix R, Y(A), X(A), deltaS(size, size, 0.0);
315
316 Matrix lastX(X);
317 Matrix lastY(Y);
318
319 for (Size i=0; i<maxIterations; ++i) {
320 R = Y - deltaS;
321 X = projectToPositiveSemidefiniteMatrix(M&: R);
322 deltaS = X - R;
323 Y = projectToUnitDiagonalMatrix(M: X);
324
325 // convergence test
326 if (std::max(a: normInf(M: X-lastX)/normInf(M: X),
327 b: std::max(a: normInf(M: Y-lastY)/normInf(M: Y),
328 b: normInf(M: Y-X)/normInf(M: Y)))
329 <= tolerance)
330 {
331 break;
332 }
333 lastX = X;
334 lastY = Y;
335 }
336
337 // ensure we return a symmetric matrix
338 for (Size i=0; i<size; ++i)
339 for (Size j=0; j<i; ++j)
340 Y[i][j] = Y[j][i];
341
342 return Y;
343 }
344 }
345
346
347 Matrix pseudoSqrt(const Matrix& matrix, SalvagingAlgorithm::Type sa) {
348 Size size = matrix.rows();
349
350 #if defined(QL_EXTRA_SAFETY_CHECKS)
351 checkSymmetry(matrix);
352 #else
353 QL_REQUIRE(size == matrix.columns(),
354 "non square matrix: " << size << " rows, " <<
355 matrix.columns() << " columns");
356 #endif
357
358 // spectral (a.k.a Principal Component) analysis
359 SymmetricSchurDecomposition jd(matrix);
360 Matrix diagonal(size, size, 0.0);
361
362 // salvaging algorithm
363 Matrix result(size, size);
364 bool negative;
365 switch (sa) {
366 case SalvagingAlgorithm::None:
367 // eigenvalues are sorted in decreasing order
368 QL_REQUIRE(jd.eigenvalues()[size-1]>=-1e-16,
369 "negative eigenvalue(s) ("
370 << std::scientific << jd.eigenvalues()[size-1]
371 << ")");
372 result = CholeskyDecomposition(m: matrix, flexible: true);
373 break;
374 case SalvagingAlgorithm::Spectral:
375 // negative eigenvalues set to zero
376 for (Size i=0; i<size; i++)
377 diagonal[i][i] =
378 std::sqrt(x: std::max<Real>(a: jd.eigenvalues()[i], b: 0.0));
379
380 result = jd.eigenvectors() * diagonal;
381 normalizePseudoRoot(matrix, pseudo&: result);
382 break;
383 case SalvagingAlgorithm::Hypersphere:
384 // negative eigenvalues set to zero
385 negative=false;
386 for (Size i=0; i<size; ++i){
387 diagonal[i][i] =
388 std::sqrt(x: std::max<Real>(a: jd.eigenvalues()[i], b: 0.0));
389 if (jd.eigenvalues()[i]<0.0) negative=true;
390 }
391 result = jd.eigenvectors() * diagonal;
392 normalizePseudoRoot(matrix, pseudo&: result);
393
394 if (negative)
395 result = hypersphereOptimize(targetMatrix: matrix, currentRoot: result, lowerDiagonal: false);
396 break;
397 case SalvagingAlgorithm::LowerDiagonal:
398 // negative eigenvalues set to zero
399 negative=false;
400 for (Size i=0; i<size; ++i){
401 diagonal[i][i] =
402 std::sqrt(x: std::max<Real>(a: jd.eigenvalues()[i], b: 0.0));
403 if (jd.eigenvalues()[i]<0.0) negative=true;
404 }
405 result = jd.eigenvectors() * diagonal;
406
407 normalizePseudoRoot(matrix, pseudo&: result);
408
409 if (negative)
410 result = hypersphereOptimize(targetMatrix: matrix, currentRoot: result, lowerDiagonal: true);
411 break;
412 case SalvagingAlgorithm::Higham: {
413 int maxIterations = 40;
414 Real tol = 1e-6;
415 result = highamImplementation(A: matrix, maxIterations, tolerance: tol);
416 result = CholeskyDecomposition(m: result, flexible: true);
417 }
418 break;
419 default:
420 QL_FAIL("unknown salvaging algorithm");
421 }
422
423 return result;
424 }
425
426
427 Matrix rankReducedSqrt(const Matrix& matrix,
428 Size maxRank,
429 Real componentRetainedPercentage,
430 SalvagingAlgorithm::Type sa) {
431 Size size = matrix.rows();
432
433 #if defined(QL_EXTRA_SAFETY_CHECKS)
434 checkSymmetry(matrix);
435 #else
436 QL_REQUIRE(size == matrix.columns(),
437 "non square matrix: " << size << " rows, " <<
438 matrix.columns() << " columns");
439 #endif
440
441 QL_REQUIRE(componentRetainedPercentage>0.0,
442 "no eigenvalues retained");
443
444 QL_REQUIRE(componentRetainedPercentage<=1.0,
445 "percentage to be retained > 100%");
446
447 QL_REQUIRE(maxRank>=1,
448 "max rank required < 1");
449
450 // spectral (a.k.a Principal Component) analysis
451 SymmetricSchurDecomposition jd(matrix);
452 Array eigenValues = jd.eigenvalues();
453
454 // salvaging algorithm
455 switch (sa) {
456 case SalvagingAlgorithm::None:
457 // eigenvalues are sorted in decreasing order
458 QL_REQUIRE(eigenValues[size-1]>=-1e-16,
459 "negative eigenvalue(s) ("
460 << std::scientific << eigenValues[size-1]
461 << ")");
462 break;
463 case SalvagingAlgorithm::Spectral:
464 // negative eigenvalues set to zero
465 for (Size i=0; i<size; ++i)
466 eigenValues[i] = std::max<Real>(a: eigenValues[i], b: 0.0);
467 break;
468 case SalvagingAlgorithm::Higham:
469 {
470 int maxIterations = 40;
471 Real tolerance = 1e-6;
472 Matrix adjustedMatrix = highamImplementation(A: matrix, maxIterations, tolerance);
473 jd = SymmetricSchurDecomposition(adjustedMatrix);
474 eigenValues = jd.eigenvalues();
475 }
476 break;
477 default:
478 QL_FAIL("unknown or invalid salvaging algorithm");
479 }
480
481 // factor reduction
482 Real enough = componentRetainedPercentage *
483 std::accumulate(first: eigenValues.begin(),
484 last: eigenValues.end(), init: Real(0.0));
485 if (componentRetainedPercentage == 1.0) {
486 // numerical glitches might cause some factors to be discarded
487 enough *= 1.1;
488 }
489 // retain at least one factor
490 Real components = eigenValues[0];
491 Size retainedFactors = 1;
492 for (Size i=1; components<enough && i<size; ++i) {
493 components += eigenValues[i];
494 retainedFactors++;
495 }
496 // output is granted to have a rank<=maxRank
497 retainedFactors=std::min(a: retainedFactors, b: maxRank);
498
499 Matrix diagonal(size, retainedFactors, 0.0);
500 for (Size i=0; i<retainedFactors; ++i)
501 diagonal[i][i] = std::sqrt(x: eigenValues[i]);
502 Matrix result = jd.eigenvectors() * diagonal;
503
504 normalizePseudoRoot(matrix, pseudo&: result);
505 return result;
506 }
507}
508

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