A correlation matrix is a symmetric matrix with unit diagonal and nonnegative eigenvalues. In 2000 I was approached by a London fund management company who wanted to find the nearest correlation matrix (NCM) in the Frobenius norm to an almost correlation matrix: a symmetric matrix having a significant number of (small) negative eigenvalues. This problem arises when the data from which the correlations are constructed is asynchronous or incomplete, or when models are stress-tested by artificially adjusting individual correlations. Solving the NCM problem (or obtaining a true correlation matrix some other way) is important in order to avoid subsequent calculations breaking down due to negative variances or volatilities, for example.
The convexity properties of the problem mean that there is a unique nearest correlation matrix, which is hence a global minimizer. In the 1990s several algorithms had been proposed for computing it, but none was guaranteed to work. Prompted by the approach from the company, I investigated the problem. I proved some results characterizing the solution and derived an alternating projections algorithm for computing it 1. The algorithm repeatedly projects onto the set of matrices with unit diagonal and the cone of symmetric positive semidefinite matrices. It is guaranteed to converge to the minimum, but does so at a linear rate. An important feature of the algorithm is that other projections can be added on. Thus, for example, if we want to leave the trailing principal submatrix of order three unchanged, we simply restore it at the end of each iteration 2, 3.
The alternating projections algorithm is widely used, but can be slow to converge, especially for large matrices 4. In 2006, Qi and Sun 5 derived a Newton method for the NCM problem. They work with the dual of the original problem, which is unconstrained. The objective function of the dual is not twice continuously differentiable, but by using the theory of strongly semismooth matrix functions Qi and Sun show that Newton’s method nevertheless has global quadratic convergence.
Ruediger Borsdorf and I, building on work in his M.Sc. thesis 3, built an algorithm that solves the Newton equations using minres with a Jacobi preconditioner (a nontrivial task since the coefficient matrix is not explicitly available), and has some other refinements described in 6. This algorithm has been implemented in the NAG Library 7.
In subsequent work, Borsdorf, Marcos Raydan and I 8 9 used the spectral projected gradient method (SPGM) to solve the k-factor NCM, in which the correlation matrix is constrained to have the form of a diagonal matrix plus a rank-k matrix. This problem variant arises in multifactor normal copula models, collateralized debt obligations (CDOs), and multivariate time series. One existing previous algorithm can fail to converge or solve the problem, but the SPGM has guaranteed convergence to a stationary point. This algorithm has also been implemented in the NAG Library.
The NCM problem has proved to be of very wide interest beyond the world of finance, as indicated by the fact that 1 is now my third best cited paper on the Web of Science. Recent applications in which the problem arises include reconstructing 20th century sea levels, genetic evaluations for thoroughbred horse breeding, modelling public health data sets, modelling storm damage of buildings, and a Kriging model for reservoirs.
I regularly receive emails asking for software implementing algorithms for the NCM problem. I thought it would be useful to summarize what is available. In general, the Newton method is preferred, but the alternating projections method is more flexible as regards incorporating additional constraints.
Original NCM Problem
- Alternating Projections Method
- Newton Method
- MATLAB: Defeng Sun, various codes.
- NAG Library (Fortran/SMP, C, NAG Toolbox for MATLAB, .NET, Java, Excel, R):
- g02aaf (Fortran), g02aac (C), g02aa.m (MATLAB), g02aa (R),
- g02abf (Fortran), g02abc (C), g02ab.m (MATLAB), g02ab (R): allows for weighted Frobenius norm based on two-side scaling and lower bounds on the eigenvalues.
- g02ajf (Fortran), g02ajf.m (MATLAB) to appear in Mark 24: allows for componentwise weighted Frobenius norm and lower bounds on the eigenvalues.
A MATLAB Alternating Projections Function
I thought it would be useful to provide my own MATLAB function
nearcorr.m implementing the alternating projections algorithm. The listing is below. To see how it compares with the NAG code
g02aa.m I ran the test code
%NEARCORR_TEST Compare g02aa and nearcorr. rng(10) % Seed random number generators. n = 100; A = gallery('randcorr',n); % Random correlation matrix. E = randn(n)*1e-1; A = A + (E + E')/2; % Perturb it. tol = 1e-10; % A = cor1399; tol = 1e-4; fprintf('g02aa:\n') maxits = int64(-1); % For linear equation solver. maxit = int64(-1); % For Newton iteration. tic [~,X1,iter1,feval,nrmgrd,ifail] = g02aa(A,'errtol',tol,'maxits',maxits, ... 'maxit',maxit); toc fprintf(' Newton steps taken: %d\n', iter1); fprintf(' Norm of gradient of last Newton step: %6.4f\n', nrmgrd); if ifail > 0, fprintf(' g02aa failed with ifail = %g\n', ifail), end fprintf('nearcorr:\n') tic [X2,iter2] = nearcorr(A,tol,,,,,1); toc fprintf(' Number of iterations: %d\n', iter2); fprintf(' Normwise relative difference between computed solutions:') fprintf('%9.2e\n', norm(X1-X2,1)/norm(X1,1))
Running under Windows 7 on an Ivy Bridge Core i7 processor @4.4Ghz I obtained the following results, where the "real-life" matrix is based on stock data:
|1. Random (100), tol = 1e-10||g02aa||0.023||4|
|2. Random (500), tol = 1e-10||g02aa||0.48||4|
|3. Real-life (1399), tol = 1e-4||g02aa||6.8||5|
The results show that while
nearcorr can be fast for small dimensions, the number of iterations, and hence its run time, tends to increase with the dimension and it can be many times slower than the Newton method. This is a stark illustration of the difference between quadratic convergence and linear (with problem-dependent constant) convergence. Here is my MATLAB function
function [X,iter] = nearcorr(A,tol,flag,maxits,n_pos_eig,w,prnt) %NEARCORR Nearest correlation matrix. % X = NEARCORR(A,TOL,FLAG,MAXITS,N_POS_EIG,W,PRNT) % finds the nearest correlation matrix to the symmetric matrix A. % TOL is a convergence tolerance, which defaults to 16*EPS. % If using FLAG == 1, TOL must be a 2-vector, with first component % the convergence tolerance and second component a tolerance % for defining "sufficiently positive" eigenvalues. % FLAG = 0: solve using full eigendecomposition (EIG). % FLAG = 1: treat as "highly non-positive definite A" and solve % using partial eigendecomposition (EIGS). % MAXITS is the maximum number of iterations (default 100, but may % need to be increased). % N_POS_EIG (optional) is the known number of positive eigenvalues of A. % W is a vector defining a diagonal weight matrix diag(W). % PRNT = 1 for display of intermediate output. % By N. J. Higham, 13/6/01, updated 30/1/13. % Reference: N. J. Higham, Computing the nearest correlation % matrix---A problem from finance. IMA J. Numer. Anal., % 22(3):329-343, 2002. if ~isequal(A,A'), error('A must by symmetric.'), end if nargin < 2 || isempty(tol), tol = length(A)*eps*[1 1]; end if nargin < 3 || isempty(flag), flag = 0; end if nargin < 4 || isempty(maxits), maxits = 100; end if nargin < 6 || isempty(w), w = ones(length(A),1); end if nargin < 7, prnt = 1; end n = length(A); if flag >= 1 if nargin < 5 || isempty(n_pos_eig) [V,D] = eig(A); d = diag(D); n_pos_eig = sum(d >= tol(2)*d(n)); end if prnt, fprintf('n = %g, n_pos_eig = %g\n', n, n_pos_eig), end end X = A; Y = A; iter = 1; rel_diffX = inf; rel_diffY = inf; rel_diffXY = inf; dS = zeros(size(A)); w = w(:); Whalf = sqrt(w*w'); while max([rel_diffX rel_diffY rel_diffXY]) > tol(1) Xold = X; R = X - dS; R_wtd = Whalf.*R; if flag == 0 X = proj_spd(R_wtd); elseif flag == 1 [X,np] = proj_spd_eigs(R_wtd,n_pos_eig,tol(2)); end X = X ./ Whalf; dS = X - R; Yold = Y; Y = proj_unitdiag(X); rel_diffX = norm(X-Xold,'fro')/norm(X,'fro'); rel_diffY = norm(Y-Yold,'fro')/norm(Y,'fro'); rel_diffXY = norm(Y-X,'fro')/norm(Y,'fro'); if prnt fprintf('%2.0f: %9.2e %9.2e %9.2e', ... iter, rel_diffX, rel_diffY, rel_diffXY) if flag >= 1, fprintf(' np = %g\n',np), else fprintf('\n'), end end iter = iter + 1; if iter > maxits error(['Stopped after ' num2str(maxits) ' its. Try increasing MAXITS.']) end X = Y; end %%%%%%%%%%%%%%%%%%%%%%%% function A = proj_spd(A) %PROJ_SPD if ~isequal(A,A'), error('Not symmetric!'), end [V,D] = eig(A); A = V*diag(max(diag(D),0))*V'; A = (A+A')/2; % Ensure symmetry. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [A,n_pos_eig_found] = proj_spd_eigs(A,n_pos_eig,tol) %PROJ_SPD_EIGS if ~isequal(A,A'), error('Not symmetric!'), end k = n_pos_eig + 10; % 10 is safety factor. if k > length(A), k = n_pos_eig; end opts.disp = 0; [V,D] = eigs(A,k,'LA',opts); d = diag(D); j = (d > tol*max(d)); n_pos_eig_found = sum(j); A = V(:,j)*D(j,j)*V(:,j)'; % Build using only the selected eigenpairs. A = (A+A')/2; % Ensure symmetry. %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A = proj_unitdiag(A) %PROJ_SPD n = length(A); A(1:n+1:n^2) = 1;
1 Nicholas J. Higham, Computing the Nearest Correlation Matrix—A Problem from Finance, IMA J. Numer. Anal. 22, 329–343, 2002.
2 Craig Lucas, Computing Nearest Covariance and Correlation Matrices, M.Sc. Thesis, University of Manchester, 2001.
3 Ruediger Borsdorf, A Newton Algorithm for the Nearest Correlation Matrix, M.Sc. Thesis, University of Manchester, 2007.
5 Hou-Duo Qi and Defeng Sun, A Quadratically Convergent Newton Method for Computing the Nearest Correlation Matrix, SIAM J. Matrix Anal. Appl. 28, 360-385, 2006
6 Ruediger Borsdorf and Nicholas J. Higham, A Preconditioned Newton Algorithm for the Nearest Correlation Matrix, IMA J. Numer. Anal. 30, 94-107, 2010.
8 Ruediger Borsdorf, Nicholas Higham and Marcos Raydan, Computing a Nearest Correlation Matrix with Factor Structure, SIAM J. Matrix Anal., Appl. 31, 2603-2622, 2010
9 Ruediger Borsdorf, Structured Matrix Nearness Problems: Theory and Algorithms, Ph.D. Thesis, University of Manchester, 2012.