Computing Nearest Correlation Matrix via Low-Rank ODE’s Based Technique

: For n -dimensional real-valued matrix A , the computation of nearest correlation matrix; that is, a symmetric, positive semi-deﬁnite, unit diagonal and off-diagonal entries between − 1 and 1 is a problem that arises in the ﬁnance industry where the correlations exist between the stocks. The proposed methodology presented in this article computes the admissible perturbation matrix and a perturbation level to shift the negative spectrum of perturbed matrix to become non-negative or strictly positive. The solution to optimization problems constructs a gradient system of ordinary differential equations that turn over the desired perturbation matrix. Numerical testing provides enough evidence for the shifting of the negative spectrum and the computation of nearest correlation matrix.


Introduction
The correlation matrices A i ∈ R n,n are real, symmetric A t i = A i , positive semi-definite λ j (A i ) ≥ 0, ∀j, with λ j being the spectrum of A i and having a unit diagonal, diag(A i ) = 1. These correlation matrices appear when one is interested in computing correlations between pairs of random variables. For example, in the finance industry, the correlation between stocks measured over a fixed amount of time, and the missing data can compromise the correlations and hence will direct to a non-positive semi-definite matrix. Similarly, a practitioner can explore the effect on a portfolio that assigns the correlations between assets differently from a certain amount of historical values, which can easily destroy the semi-definiteness of the matrix (or matrices); for more details, refer to [1][2][3].
The approximation of correlation matrices has led much intention and interest in finding a nearest correlation matrix X that minimizes the matrix Frobenius norm of A − X in matrix nearness problem min{ A − X F : X = X t , X ≥ 0, diag(X) = 1}, where for X, Y being as the symmetric matrices and the quantity X − Y is positive semi-definite matrix and X F = trace(X t X) as R n,n represents Hilbert space with X, Y = trace(X t Y). The constraints in above minimization problem are closed and convex sets so that the matrix nearness problem has a unique solution [4].
The matrix nearness problems have been extensively studied over the last thirty-two years. Much of the literature is available about some ad hoc mathematical methods that are not guaranteed to give the best solution to the problem. The method given by Knal and ten Burge [5] decomposes the matrix X as X = Y t Y and then minimizes the objective function over Y s columns. Lurie and Goldberg [6] minimize the quantity A − R t R 2 F with R being as an upper triangular matrix having a unit 2-norm for it's columns by using the Gauss-Newton method. In [7], Higham proposed an alternating projection algorithm using convex analysis, which has a linear convergence towards the solution of matrix nearness problems.
A more general matrix nearness problem than the above problem was studied in [8] by Malick, where he had used convex set and general linear constraints to replace positive semi-definiteness and unit diag(X). Malick had applied the quasi-Newton method to dual problem after dualize the linear constraints on the diag(X). A quadratically convergent Newton method for solving matrix nearness problem was studied by Qi and Sun [9], where they proposed a dual problem to matrix nearness problem and then proposed its solution. The numerical method developed by Toh [10] demands the solution of dense linear systems with dimension n 2 2 , mathematical construction of preconditioners and then apply them to compute the nearest correlation matrix X to solve the matrix nearness problem. A globally convergent Newton's method proposed by Qi and Sun [9] computes the nearest correlation matrix X for the above matrix nearness problem.
In this work, we propose a low-rank ODE's based method to compute a nearest correlation matrix X for a given n-dimensional matrix A. Our method works in two ways: first, it allows us to shift the smallest negative eigenvalue λ 1 of the perturbed matrix A + E, such that it becomes non-negative, that is, λ 1 ≥ 0. We compute the perturbation matrix E with diag(E) = 0 by constructing and then solving an associated optimization problem. Secondly, the construction and solution to optimization problem allow to shift all negative eigenvalues λ i , ∀ i = n − 1 of given matrix A such that its spectrum becomes non-negative, that is, λ i ≥ 0, ∀ i = 1 : n − 1.

Overview of the Article
In Section 2, we provide the preliminaries of our article. We give the definitions of symmetric, positive semi-definite matrix, a correlation matrix, and nearest correlation matrix. Section 3 of our article is devoted to the problem under consideration. We briefly discuss the aim of matrix nearness problem. Section 4 of our article is dedicated to the computation of a gradient system of ordinary differential equations to localize the smallest negative eigenvalue λ 1 from the spectrum of an admissible perturbed matrix (A + E(t))) for given A ∈ R n,n .
The localization of eigenvalues λ 1 , λ 2 simultaneously from the spectrum of admissible perturbation matrix is addressed in Section 5. Finally, numerical experimentation and conclusion are presented in Sections 6 and 7, respectively.

Definition 3.
A square symmetric matrix A ∈ R n,n is called correlation matrix if (i, j)-entry is the correlation between the columns i and j of A.

Matrix Nearness Problem
The matrix nearness problem aims the computation of a matrix X, which is symmetric, positive semi-definite, unit diagonal and having off diagonal entries between −1 and 1 for a given matrix A ∈ R n,n . The matrix A is not necessary to have it's entries belonging to interval [−1, 1]. In fact, it can have real random entries. It is possible that all the eigenvalues of A are negative, or some of them are negative while others are positive. The negative eigenvalues are the basic hindrance to make A a nearest correlation matrix. Moreover, one must remain careful about the structure of A. The structure of A must have all properties, as defined in the preliminaries Section. For this purpose, we aim to compute the perturbation matrix E = E(t), ∀t having zero diagonal and a Frobenius norm bounded above by 1. In the very next Section, we present our methodology to compute the matrix E, which allows us to shift all negative eigenvalues of given A to make them non-negative.

Localization of Smallest Negative Eigenvalue λ 1
In this section, we aim to localize the smallest negative eigenvalue λ 1 of given matrix A ∈ R n,n . This involves shifting the smallest negative eigenvalue of perturbation matrix (A + E) with > 0, a small positive parameter. The matrix E has structure with diag(E) = 0. Moreover, i,j ≤ 1, ∀t.

Construction of Perturbation Matrix
We compute the perturbation matrix E(t) and then determine the direction Z =Ė(t). The computation ofĖ(t) indicates that how fast the smallest negative eigenvalue λ 1 (t) grows, that is, with η(t) being an eigenvector for smallest negative eigenvalue By differentiating Equation (1) Multiplying both sides with η(t) * (t) gives us As we know that, This further implies that, Thus, we have Using Equation (2), we have that

Formulation of Optimization Problem
The optimization problem given below allows one to determine the direction Z =Ė(t) so that the solution of the system of ODE's indicates the maximum growth of the smallest negative eigenvalue λ(t).
with η 1 ∈ R n,1 is eigenvector corresponding to smallest negative eigenvalue λ 1 . The notation * represents complex conjugate transpose. The solution of the maximization problem addressed in Equation (5) is given as the following.

Lemma 4.2.1.
Consider that E(t) is a non-zero matrix with an admissible perturbation matrix Let η 1 , η * 1 are non-zero eigen vectors associated with smallest negative eigenvalue λ 1 (t). Then solution Z to maximization problem in Equation (5) is of the form having Proj(·) as projection of direction matrix Z onto manifold of matrices E(t).
Proof. The proof is based on the idea of computation of orthogonal projection on manifold of matrices, and for this we refer to [11].

A Gradient System of ODE's
The solution matrix of the maximization problem addressed in Equation (5) allows following gradient system of ODEs on manifold of matrices E(t),

Characterization of Gradient System of ODE's
The solution of gradient system of ODE's addressed in Equation (7) possesses the following characteristic properties: Our goal in this Section is to transfer simultaneously the smallest negative eigenvalue λ 1 (t) and the negative eigenvalue λ 2 (t) which is next to λ 1 (t) from spectrum of perturbation matrix (A + E(t)) so that λ 1 (t), λ 2 (t) becomes strictly positive.

Optimization Problem
The following maximization problem allows one to determine direction matrix Z =Ė(t) so that a solution associated with gradient system of ODEs obtained after solving the problem, which indicates maximum growth for λ 1 (t) and λ 2 (t) respectively. The maximization problem, in order to have a simultaneous maximum growth for both λ 1 (t) and λ 2 (t), is given as Now, our aim is to give the solution corresponding to maximization problem addressed in Equation (8).

A Gradient System of ODE's
The optimal solution to maximization addressed in Equation (8) is given by gradient system of ODE's asĖ The gradient system of ODEs addressed in Equation (9) as a function K(t) with For parameter 0 ≤ , we have that 0 = diag(A − I) ≤ . Next, we fix e 11 (t), e 12 (t), . . . , e nn (t) and for parameter >> 0 , we obtain The admissible perturbation matrix A + E(0)) takes the form As, e 11 (t) = 1−a 11 , e 12 (t) = 1−a 22 , . . . , e nn (t) = 1−a nn . The initial perturbation matrix E(0) is given as The initial perturbation matrix E(0) is decomposed into K T (0) and K(0), the upper and lower triangular matrices respectively as, In a similar manner, the perturbation matrix E(t) has the structure The maximization of eigenvalues λ 1 (t) and λ 2 (t) demands the computation of the perturbation matrixĖ(t) and projection of Z onto manifold of family of matrices E(t). The change in the admissible perturbation matrix E(t) is given aṡ The computationK(t) cause maximum growth of λ min (t) in the maximization problem as, The solution to maximization problem addressed in Equation (12) is given by ODEĖ(t) as, having and has the form The matrix Proj(ηη * ) is decomposed into U(t), an upper triangular matrix and lower triangular matrix L T (t) as  (8) has the forṁ

Remark 1. The upper triangular matrix U(t) = (A + E(t)). The optimal solution to maximization problem addressed in Equation
The solutionĖ(t) addressed in Equation (13) is determined with the help of Euler's method, that is, Finally, the eigenvalue problem computes the non-negative spectrum for the perturbed system and it show that the matrix X is a nearest correlation matrix as it is symmetric, positive semi-definite, unit diagonal and off diagonal entries lies in the interval −1 and 1.

Numerical Experimentation
This section presents numerical experimentation for the computation of nearest correlation matrices for a given matrix A. For simplicity, we take A ∈ R n,n .

Example 1. We consider a ten-dimensional matrix A as
The eigenvalues of A are {−5.7249, −2.8009, −0.7529, 0, 1, 1, 2, 2.9796, 5.0768, 7.2223} which clearly contains the negative eigenvalues. The perturbation matrix E has a zero diagonal and help us to shift the negative eigenvalues of perturbed matrix (A + E). The perturbation matrix E is computed as We compute matrix B 1 , which alters the eigenvalues of given matrix A; at the initial stage, this matrix also has three negative eigenvalues but different from those corresponding to A. The matrix B 1 and its eigenvalues are computed as  and clearly is not a nearest correlation matrix. Matrix B 2 is also not a nearest correlation matrix but it has shifted the eigenvalues of A + E such that only one eigenvalue is negative from the spectrum. The matrix B 2 and its spectrum are given as The perturbation matrix E has a zero diagonal and helps us to shift the negative eigenvalues of perturbed matrix (A + E). The perturbation matrix E is computed as We compute matrix B 1 , which alters the eigenvalues of given matrix A; at the initial stage, this matrix also has three negative eigenvalues but different from those corresponding to A. The matrix B 1 and its eigenvalues are computed as

Conclusions
In this article, we have presented a low-rank ordinary differential equations-based technique to compute the nearest correlation matrix; that is, symmetric, positive semidefinite, unit diagonal and off-diagonal entries between −1 and 1 for a given n-dimensional real-valued matrix. The proposed methodology is based on transforming the negative spectrum of the original matrix, which involves the computation of perturbation level and an admissible perturbation having zero diagonal. The numerical experimentation turns over the desired perturbations for the computation of nearest correlation matrix for randomly chosen real-valued matrices.