1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2009 Jose Aparicio
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
20#include <ql/math/matrixutilities/factorreduction.hpp>
21#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
22#include <vector>
23
24namespace QuantLib {
25
26 std::vector<Real> factorReduction(Matrix mtrx,
27 Size maxIters) {
28 static Real tolerance = 1.e-6;
29
30 QL_REQUIRE(mtrx.rows() == mtrx.columns(),
31 "Input matrix is not square");
32
33 const Size n = mtrx.columns();
34
35 #if defined(QL_EXTRA_SAFETY_CHECKS)
36 // check symmetry
37 for (Size iRow=0; iRow<mtrx.rows(); iRow++)
38 for (Size iCol=0; iCol<iRow; iCol++)
39 QL_REQUIRE(mtrx[iRow][iCol] == mtrx[iCol][iRow],
40 "input matrix is not symmetric");
41 QL_REQUIRE(*std::max_element(mtrx.begin(), mtrx.end()) <=1 &&
42 *std::min_element(mtrx.begin(), mtrx.end()) >=-1,
43 "input matrix data is not correlation data");
44 #endif
45
46 // Initial guess value
47 std::vector<Real> previousCorrels(n, 0.);
48 for(Size iCol=0; iCol<n; iCol++) {
49 for(Size iRow=0; iRow<n; iRow++)
50 previousCorrels[iCol] +=
51 mtrx[iRow][iCol] * mtrx[iRow][iCol];
52 previousCorrels[iCol] =
53 std::sqrt(x: (previousCorrels[iCol]-1.)/(n-1.));
54 }
55
56 // iterative solution
57 Size iteration = 0;
58 Real distance;
59 do {
60 // patch Matrix diagonal
61 for(Size iCol=0; iCol<n; iCol++)
62 mtrx[iCol][iCol] =
63 previousCorrels[iCol];
64 // compute eigenvector decomposition
65 SymmetricSchurDecomposition ssDec(mtrx);
66 //const Matrix& eigenVect = ssDec.eigenvectors();
67 const Array& eigenVals = ssDec.eigenvalues();
68 Array::const_iterator itMaxEval =
69 std::max_element(first: eigenVals.begin(), last: eigenVals.end());
70 // We do not need the max value, only the position of the
71 // corresponding eigenvector
72 Size iMax = std::distance(first: eigenVals.begin(), last: itMaxEval);
73 std::vector<Real> newCorrels, distances;
74 for(Size iCol=0; iCol<n; iCol++) {
75 Real thisCorrel = mtrx[iMax][iCol];
76 newCorrels.push_back(x: thisCorrel);
77 // strictly is:
78 // abs(\sqrt{\rho}- \sqrt{\rho_{old}})/\sqrt{\rho_{old}}
79 distances.push_back(
80 x: std::abs(x: thisCorrel - previousCorrels[iCol])/
81 previousCorrels[iCol]);
82 }
83 previousCorrels = newCorrels;
84 distance = *std::max_element(first: distances.begin(), last: distances.end());
85 }while(distance > tolerance && ++iteration <= maxIters );
86
87 // test it did not go up to the max iteration and the matrix can
88 // be reduced to one factor.
89 QL_REQUIRE(iteration < maxIters,
90 "convergence not reached after " <<
91 iteration << " iterations");
92
93 #if defined(QL_EXTRA_SAFETY_CHECKS)
94 QL_REQUIRE(*std::max_element(previousCorrels.begin(),
95 previousCorrels.end()) <=1 &&
96 *std::min_element(previousCorrels.begin(),
97 previousCorrels.end()) >=-1,
98 "matrix can not be decomposed to a single factor dependence");
99 #endif
100
101 return previousCorrels;
102 }
103
104}
105

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