14910176578

Committed 08 May 2025 03:28PM UTC coverage: 73.315% (+0.02%) from 73.3%

Build # 14910176578

Build Type

Pull #2195

github

Committed by

web-flow

Commit Message

Merge 3a61f499c into 5d972fb7b

Pull Request Pull Request #2195: Added `Handle<Quote>` for spread in `OISRateHelper`

Run Details

32 of 33 new or added lines in 2 files covered. (96.97%)

277 existing lines in 25 files now uncovered.

56277 of 76761 relevant lines covered (73.31%)

8687029.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

43.19

/ql/math/matrixutilities/pseudosqrt.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003, 2004, 2007 Ferdinando Ametrano
 Copyright (C) 2006 Yiping Chen
 Copyright (C) 2007 Neil Firth

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <algorithm>
#include <ql/math/comparison.hpp>
#include <ql/math/matrixutilities/choleskydecomposition.hpp>
#include <ql/math/matrixutilities/pseudosqrt.hpp>
#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
#include <ql/math/optimization/conjugategradient.hpp>
#include <ql/math/optimization/constraint.hpp>
#include <ql/math/optimization/problem.hpp>
#include <utility>

namespace QuantLib {

    namespace {

        #if defined(QL_EXTRA_SAFETY_CHECKS)
        void checkSymmetry(const Matrix& matrix) {
            Size size = matrix.rows();
            QL_REQUIRE(size == matrix.columns(),
                       "non square matrix: " << size << " rows, " <<
                       matrix.columns() << " columns");
            for (Size i=0; i<size; ++i)
                for (Size j=0; j<i; ++j)
                    QL_REQUIRE(close(matrix[i][j], matrix[j][i]),
                               "non symmetric matrix: " <<
                               "[" << i << "][" << j << "]=" << matrix[i][j] <<
                               ", [" << j << "][" << i << "]=" << matrix[j][i]);
        }
        #endif

        void normalizePseudoRoot(const Matrix& matrix,
                                 Matrix& pseudo) {
            Size size = matrix.rows();
            QL_REQUIRE(size == pseudo.rows(),
                       "matrix/pseudo mismatch: matrix rows are " << size <<
                       " while pseudo rows are " << pseudo.columns());
            Size pseudoCols = pseudo.columns();

            // row normalization
            for (Size i=0; i<size; ++i) {
                Real norm = 0.0;
                for (Size j=0; j<pseudoCols; ++j)
                    norm += pseudo[i][j]*pseudo[i][j];
                if (norm>0.0) {
                    Real normAdj = std::sqrt(matrix[i][i]/norm);
                    for (Size j=0; j<pseudoCols; ++j)
                        pseudo[i][j] *= normAdj;
                }
            }


        }

        //cost function for hypersphere and lower-diagonal algorithm
        class HypersphereCostFunction : public CostFunction {
          private:
            Size size_;
            bool lowerDiagonal_;
            Matrix targetMatrix_;
            Array targetVariance_;
            mutable Matrix currentRoot_, tempMatrix_, currentMatrix_;
          public:
            HypersphereCostFunction(const Matrix& targetMatrix,
                                    Array targetVariance,
                                    bool lowerDiagonal)
            : size_(targetMatrix.rows()), lowerDiagonal_(lowerDiagonal),
              targetMatrix_(targetMatrix), targetVariance_(std::move(targetVariance)),
              currentRoot_(size_, size_), tempMatrix_(size_, size_), currentMatrix_(size_, size_) {}
            Array values(const Array&) const override {
                QL_FAIL("values method not implemented");
            }
            Real value(const Array& x) const override {
                Size i,j,k;
                std::fill(currentRoot_.begin(), currentRoot_.end(), 1.0);
                if (lowerDiagonal_) {
                    for (i=0; i<size_; i++) {
                        for (k=0; k<size_; k++) {
                            if (k>i) {
                                currentRoot_[i][k]=0;
                            } else {
                                for (j=0; j<=k; j++) {
                                    if (j == k && k!=i)
                                        currentRoot_[i][k] *=
                                            std::cos(x[i*(i-1)/2+j]);
                                    else if (j!=i)
                                        currentRoot_[i][k] *=
                                            std::sin(x[i*(i-1)/2+j]);
                                }
                            }
                        }
                    }
                } else {
                    for (i=0; i<size_; i++) {
                        for (k=0; k<size_; k++) {
                            for (j=0; j<=k; j++) {
                                if (j == k && k!=size_-1)
                                    currentRoot_[i][k] *=
                                        std::cos(x[j*size_+i]);
                                else if (j!=size_-1)
                                    currentRoot_[i][k] *=
                                        std::sin(x[j*size_+i]);
                            }
                        }
                    }
                }
                Real temp, error=0;
                tempMatrix_ = transpose(currentRoot_);
                currentMatrix_ = currentRoot_ * tempMatrix_;
                for (i=0;i<size_;i++) {
                    for (j=0;j<size_;j++) {
                        temp = currentMatrix_[i][j]*targetVariance_[i]
                          *targetVariance_[j]-targetMatrix_[i][j];
                        error += temp*temp;
                    }
                }
                return error;
            }
        };

        // Optimization function for hypersphere and lower-diagonal algorithm
        Matrix hypersphereOptimize(const Matrix& targetMatrix,
                                   const Matrix& currentRoot,
                                   const bool lowerDiagonal) {
            Size i,j,k,size = targetMatrix.rows();
            Matrix result(currentRoot);
            Array variance(size, 0);
            for (i=0; i<size; i++){
                variance[i]=std::sqrt(targetMatrix[i][i]);
            }
            if (lowerDiagonal) {
                Matrix approxMatrix(result*transpose(result));
                result = CholeskyDecomposition(approxMatrix, true);
                for (i=0; i<size; i++) {
                    for (j=0; j<size; j++) {
                        result[i][j]/=std::sqrt(approxMatrix[i][i]);
                    }
                }
            } else {
                for (i=0; i<size; i++) {
                    for (j=0; j<size; j++) {
                        result[i][j]/=variance[i];
                    }
                }
            }

            ConjugateGradient optimize;
            EndCriteria endCriteria(100, 10, 1e-8, 1e-8, 1e-8);
            HypersphereCostFunction costFunction(targetMatrix, variance,
                                                 lowerDiagonal);
            NoConstraint constraint;

            // hypersphere vector optimization

            if (lowerDiagonal) {
                Array theta(size * (size-1)/2);
                const Real eps=1e-16;
                for (i=1; i<size; i++) {
                    for (j=0; j<i; j++) {
                        theta[i*(i-1)/2+j]=result[i][j];
                        if (theta[i*(i-1)/2+j]>1-eps)
                            theta[i*(i-1)/2+j]=1-eps;
                        if (theta[i*(i-1)/2+j]<-1+eps)
                            theta[i*(i-1)/2+j]=-1+eps;
                        for (k=0; k<j; k++) {
                            theta[i*(i-1)/2+j] /= std::sin(theta[i*(i-1)/2+k]);
                            if (theta[i*(i-1)/2+j]>1-eps)
                                theta[i*(i-1)/2+j]=1-eps;
                            if (theta[i*(i-1)/2+j]<-1+eps)
                                theta[i*(i-1)/2+j]=-1+eps;
                        }
                        theta[i*(i-1)/2+j] = std::acos(theta[i*(i-1)/2+j]);
                        if (j==i-1) {
                            if (result[i][i]<0)
                                theta[i*(i-1)/2+j]=-theta[i*(i-1)/2+j];
                        }
                    }
                }
                Problem p(costFunction, constraint, theta);
                optimize.minimize(p, endCriteria);
                theta = p.currentValue();
                std::fill(result.begin(),result.end(),1.0);
                for (i=0; i<size; i++) {
                    for (k=0; k<size; k++) {
                        if (k>i) {
                            result[i][k]=0;
                        } else {
                            for (j=0; j<=k; j++) {
                                if (j == k && k!=i)
                                    result[i][k] *=
                                        std::cos(theta[i*(i-1)/2+j]);
                                else if (j!=i)
                                    result[i][k] *=
                                        std::sin(theta[i*(i-1)/2+j]);
                            }
                        }
                    }
                }
            } else {
                Array theta(size * (size-1));
                const Real eps=1e-16;
                for (i=0; i<size; i++) {
                    for (j=0; j<size-1; j++) {
                        theta[j*size+i]=result[i][j];
                        if (theta[j*size+i]>1-eps)
                            theta[j*size+i]=1-eps;
                        if (theta[j*size+i]<-1+eps)
                            theta[j*size+i]=-1+eps;
                        for (k=0;k<j;k++) {
                            theta[j*size+i] /= std::sin(theta[k*size+i]);
                            if (theta[j*size+i]>1-eps)
                                theta[j*size+i]=1-eps;
                            if (theta[j*size+i]<-1+eps)
                                theta[j*size+i]=-1+eps;
                        }
                        theta[j*size+i] = std::acos(theta[j*size+i]);
                        if (j==size-2) {
                            if (result[i][j+1]<0)
                                theta[j*size+i]=-theta[j*size+i];
                        }
                    }
                }
                Problem p(costFunction, constraint, theta);
                optimize.minimize(p, endCriteria);
                theta=p.currentValue();
                std::fill(result.begin(),result.end(),1.0);
                for (i=0; i<size; i++) {
                    for (k=0; k<size; k++) {
                        for (j=0; j<=k; j++) {
                            if (j == k && k!=size-1)
                                result[i][k] *= std::cos(theta[j*size+i]);
                            else if (j!=size-1)
                                result[i][k] *= std::sin(theta[j*size+i]);
                        }
                    }
                }
            }

            for (i=0; i<size; i++) {
                for (j=0; j<size; j++) {
                    result[i][j]*=variance[i];
                }
            }
            return result;
        }

        // Matrix infinity norm. See Golub and van Loan (2.3.10) or
        // <http://en.wikipedia.org/wiki/Matrix_norm>
        Real normInf(const Matrix& M) {
            Size rows = M.rows();
            Size cols = M.columns();
            Real norm = 0.0;
            for (Size i=0; i<rows; ++i) {
                Real colSum = 0.0;
                for (Size j=0; j<cols; ++j)
                    colSum += std::fabs(M[i][j]);
                norm = std::max(norm, colSum);
            }
            return norm;
        }

        // Take a matrix and make all the diagonal entries 1.
        Matrix projectToUnitDiagonalMatrix(const Matrix& M) {
            Size size = M.rows();
            QL_REQUIRE(size == M.columns(),
                       "matrix not square");

            Matrix result(M);
            for (Size i=0; i<size; ++i)
                result[i][i] = 1.0;

            return result;
        }

        // Take a matrix and make all the eigenvalues non-negative
        Matrix projectToPositiveSemidefiniteMatrix(Matrix& M) {
            Size size = M.rows();
            QL_REQUIRE(size == M.columns(),
                       "matrix not square");

            Matrix diagonal(size, size, 0.0);
            SymmetricSchurDecomposition jd(M);
            for (Size i=0; i<size; ++i)
                diagonal[i][i] = std::max<Real>(jd.eigenvalues()[i], 0.0);

            Matrix result =
                jd.eigenvectors()*diagonal*transpose(jd.eigenvectors());
            return result;
        }

        // implementation of the Higham algorithm to find the nearest
        // correlation matrix.
        Matrix highamImplementation(const Matrix& A, const Size maxIterations, const Real& tolerance) {

            Size size = A.rows();
            Matrix R, Y(A), X(A), deltaS(size, size, 0.0);

            Matrix lastX(X);
            Matrix lastY(Y);

            for (Size i=0; i<maxIterations; ++i) {
                R = Y - deltaS;
                X = projectToPositiveSemidefiniteMatrix(R);
                deltaS = X - R;
                Y = projectToUnitDiagonalMatrix(X);

                // convergence test
                if (std::max({
                            normInf(X-lastX)/normInf(X),
                            normInf(Y-lastY)/normInf(Y),
                            normInf(Y-X)/normInf(Y)
                        })
                    <= tolerance)
                {
                    break;
                }
                lastX = X;
                lastY = Y;
            }

            // ensure we return a symmetric matrix
            for (Size i=0; i<size; ++i)
                for (Size j=0; j<i; ++j)
                    Y[i][j] = Y[j][i];

            return Y;
        }
    }


    Matrix pseudoSqrt(const Matrix& matrix, SalvagingAlgorithm::Type sa) {
        Size size = matrix.rows();

        #if defined(QL_EXTRA_SAFETY_CHECKS)
        checkSymmetry(matrix);
        #else
        QL_REQUIRE(size == matrix.columns(),
                   "non square matrix: " << size << " rows, " <<
                   matrix.columns() << " columns");
        #endif

        // spectral (a.k.a Principal Component) analysis
        SymmetricSchurDecomposition jd(matrix);
        Matrix diagonal(size, size, 0.0);

        // salvaging algorithm
        Matrix result(size, size);
        bool negative;
        switch (sa) {
          case SalvagingAlgorithm::None:
            // eigenvalues are sorted in decreasing order
            QL_REQUIRE(jd.eigenvalues()[size-1]>=-1e-16,
                       "negative eigenvalue(s) ("
                       << std::scientific << jd.eigenvalues()[size-1]
                       << ")");
            result = CholeskyDecomposition(matrix, true);
            break;
          case SalvagingAlgorithm::Spectral:
            // negative eigenvalues set to zero
            for (Size i=0; i<size; i++)
                diagonal[i][i] =
                    std::sqrt(std::max<Real>(jd.eigenvalues()[i], 0.0));

            result = jd.eigenvectors() * diagonal;
            normalizePseudoRoot(matrix, result);
            break;
          case SalvagingAlgorithm::Hypersphere:
            // negative eigenvalues set to zero
            negative=false;
            for (Size i=0; i<size; ++i){
                diagonal[i][i] =
                    std::sqrt(std::max<Real>(jd.eigenvalues()[i], 0.0));
                if (jd.eigenvalues()[i]<0.0) negative=true;
            }
            result = jd.eigenvectors() * diagonal;
            normalizePseudoRoot(matrix, result);

            if (negative)
                result = hypersphereOptimize(matrix, result, false);
            break;
          case SalvagingAlgorithm::LowerDiagonal:
            // negative eigenvalues set to zero
            negative=false;
            for (Size i=0; i<size; ++i){
                diagonal[i][i] =
                    std::sqrt(std::max<Real>(jd.eigenvalues()[i], 0.0));
                if (jd.eigenvalues()[i]<0.0) negative=true;
            }
            result = jd.eigenvectors() * diagonal;

            normalizePseudoRoot(matrix, result);

            if (negative)
                result = hypersphereOptimize(matrix, result, true);
            break;
          case SalvagingAlgorithm::Higham: {
              int maxIterations = 40;
              Real tol = 1e-6;
              result = highamImplementation(matrix, maxIterations, tol);
              result = CholeskyDecomposition(result, true);
            }
            break;
          case SalvagingAlgorithm::Principal: {
              QL_REQUIRE(jd.eigenvalues().back()>=-10*QL_EPSILON,
                         "negative eigenvalue(s) ("
                         << std::scientific << jd.eigenvalues().back()
                         << ")");

              Array sqrtEigenvalues(size);
              std::transform(
                  jd.eigenvalues().begin(), jd.eigenvalues().end(),
                  sqrtEigenvalues.begin(),
                  [](Real lambda) -> Real {
                      return std::sqrt(std::max<Real>(lambda, 0.0));
                  }
              );

              for (Size i=0; i < size; ++i)
                  std::transform(
                       sqrtEigenvalues.begin(), sqrtEigenvalues.end(),
                       jd.eigenvectors().row_begin(i),
                       diagonal.column_begin(i),
                       std::multiplies<>()
                   );

              result = jd.eigenvectors()*diagonal;
              result = 0.5*(result + transpose(result));
            }
            break;
          default:
            QL_FAIL("unknown salvaging algorithm");
        }

        return result;
    }


    Matrix rankReducedSqrt(const Matrix& matrix,
                           Size maxRank,
                           Real componentRetainedPercentage,
                           SalvagingAlgorithm::Type sa) {
        Size size = matrix.rows();

        #if defined(QL_EXTRA_SAFETY_CHECKS)
        checkSymmetry(matrix);
        #else
        QL_REQUIRE(size == matrix.columns(),
                   "non square matrix: " << size << " rows, " <<
                   matrix.columns() << " columns");
        #endif

        QL_REQUIRE(componentRetainedPercentage>0.0,
                   "no eigenvalues retained");

        QL_REQUIRE(componentRetainedPercentage<=1.0,
                   "percentage to be retained > 100%");

        QL_REQUIRE(maxRank>=1,
                   "max rank required < 1");

        // spectral (a.k.a Principal Component) analysis
        SymmetricSchurDecomposition jd(matrix);
        Array eigenValues = jd.eigenvalues();

        // salvaging algorithm
        switch (sa) {
          case SalvagingAlgorithm::None:
            // eigenvalues are sorted in decreasing order
            QL_REQUIRE(eigenValues[size-1]>=-1e-16,
                       "negative eigenvalue(s) ("
                       << std::scientific << eigenValues[size-1]
                       << ")");
            break;
          case SalvagingAlgorithm::Spectral:
            // negative eigenvalues set to zero
            for (Size i=0; i<size; ++i)
                eigenValues[i] = std::max<Real>(eigenValues[i], 0.0);
            break;
          case SalvagingAlgorithm::Higham:
              {
                  int maxIterations = 40;
                  Real tolerance = 1e-6;
                  Matrix adjustedMatrix = highamImplementation(matrix, maxIterations, tolerance);
                  jd = SymmetricSchurDecomposition(adjustedMatrix);
                  eigenValues = jd.eigenvalues();
              }
              break;
          default:
            QL_FAIL("unknown or invalid salvaging algorithm");
        }

        // factor reduction
        Real enough = componentRetainedPercentage *
                      std::accumulate(eigenValues.begin(),
                                      eigenValues.end(), Real(0.0));
        if (componentRetainedPercentage == 1.0) {
            // numerical glitches might cause some factors to be discarded
            enough *= 1.1;
        }
        // retain at least one factor
        Real components = eigenValues[0];
        Size retainedFactors = 1;
        for (Size i=1; components<enough && i<size; ++i) {
            components += eigenValues[i];
            retainedFactors++;
        }
        // output is granted to have a rank<=maxRank
        retainedFactors=std::min(retainedFactors, maxRank);

        Matrix diagonal(size, retainedFactors, 0.0);
        for (Size i=0; i<retainedFactors; ++i)
            diagonal[i][i] = std::sqrt(eigenValues[i]);
        Matrix result = jd.eigenvectors() * diagonal;

        normalizePseudoRoot(matrix, result);
        return result;
    }
}

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

lballabio / QuantLib / 14910176578

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous