# NAG FL Interfaceg02akf (corrmat_​nearest_​rank)

## 1Purpose

g02akf computes the nearest correlation matrix of maximum prescribed rank, in the Frobenius norm, to a given square, input matrix.

## 2Specification

Fortran Interface
 Subroutine g02akf ( g, ldg, n, rank, x, ldx, f, nsub,
 Integer, Intent (In) :: ldg, n, rank, maxits, maxit, ldx Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: nsub Real (Kind=nag_wp), Intent (In) :: errtol, ranktol Real (Kind=nag_wp), Intent (Inout) :: g(ldg,n), x(ldx,n) Real (Kind=nag_wp), Intent (Out) :: f, rankerr
#include <nag.h>
 void g02akf_ (double g[], const Integer *ldg, const Integer *n, const Integer *rank, const double *errtol, const double *ranktol, const Integer *maxits, const Integer *maxit, double x[], const Integer *ldx, double *f, double *rankerr, Integer *nsub, Integer *ifail)
The routine may be called by the names g02akf or nagf_correg_corrmat_nearest_rank.

## 3Description

g02akf finds the nearest correlation matrix $X$ of maximum prescribed rank $r$ to an approximate correlation matrix $G$ in the Frobenius norm.
The solver is based on the Majorized Penalty Approach (MPA) proposed by Gao and Sun (2010). One of the key elements in this type of method is that the subproblems are similar to the nearest correlation matrix problem without rank constraint, and can be solved efficiently by g02aaf. The total number of subproblems solved is controlled by the arguments maxit and maxits. The algorithm behaviour and solver accuracy can be modified by these and other input arguments. The default values for these arguments are chosen to work well in the general case but it is recommended that you tune them to your particular problem. For a detailed description of the algorithm see Section 11.
Bai S, Qi H–D and Xiu N (2015) Constrained best Euclidean distance embedding on a sphere: A matrix optimization approach SIAM J. Optim. 25(1) 439–467
Gao Y and Sun D (2010) A majorized penalty approach for calibrating rank constrained correlation matrix problems Technical report Department of Mathematics, National University of Singapore
Qi H–D and Yuan X (2014) Computing the nearest Euclidean distance matrix with low embedding dimensions Mathematical Programming 147(1–2) 351–389

## 5Arguments

1: $\mathbf{g}\left({\mathbf{ldg}},{\mathbf{n}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: $\stackrel{~}{G}$, the initial matrix.
On exit: a symmetric matrix $G=\frac{1}{2}\left(\stackrel{~}{G}+{\stackrel{~}{G}}^{\mathrm{T}}\right)$ with diagonal elements set to $1.0$.
2: $\mathbf{ldg}$Integer Input
On entry: the first dimension of the array g as declared in the (sub)program from which g02akf is called.
Constraint: ${\mathbf{ldg}}\ge {\mathbf{n}}$.
3: $\mathbf{n}$Integer Input
On entry: the order of the matrix $G$.
Constraint: ${\mathbf{n}}>0$.
4: $\mathbf{rank}$Integer Input
On entry: $r$, the upper bound for the rank of $X$.
Constraint: $0<{\mathbf{rank}}\le {\mathbf{n}}$.
5: $\mathbf{errtol}$Real (Kind=nag_wp) Input
On entry: the termination tolerance for the convergence measure of the objective function value.
If ${\mathbf{errtol}}\le 0.0$, then a value of $\text{1.0E−5}$ is used. See Section 11 for further details.
6: $\mathbf{ranktol}$Real (Kind=nag_wp) Input
On entry: the feasibility tolerance for the rank constraint.
If ${\mathbf{ranktol}}\le 0.0$, is used. See Section 11 for further details.
7: $\mathbf{maxits}$Integer Input
On entry: specifies the maximum number of iterations used for the majorization approach to solve penalized problems at each level of penalty parameter.
If ${\mathbf{maxits}}\le 0$, then a value of $100$ is used.
8: $\mathbf{maxit}$Integer Input
On entry: specifies the maximum number of iterations for the penalty method, i.e., the maximum level of penalty parameter.
If ${\mathbf{maxit}}\le 0$, then a value of $200$ is used.
9: $\mathbf{x}\left({\mathbf{ldx}},{\mathbf{n}}\right)$Real (Kind=nag_wp) array Output
On exit: $X$, the nearest correlation matrix of rank $r$.
10: $\mathbf{ldx}$Integer Input
On entry: the first dimension of the array x as declared in the (sub)program from which g02akf is called.
Constraint: ${\mathbf{ldx}}\ge {\mathbf{n}}$.
11: $\mathbf{f}$Real (Kind=nag_wp) Output
On exit: the difference between $X$ and $G$ given by $\frac{1}{2}{‖X-G‖}_{F}^{2}$.
12: $\mathbf{rankerr}$Real (Kind=nag_wp) Output
On exit: the rank error of $X$, defined as ${\sum }_{i=r+1}^{n}\left({\lambda }_{i}\right)$, given that ${\lambda }_{i}$ denote eigenvalues of $X$ sorted in non-increasing order.
13: $\mathbf{nsub}$Integer Output
On exit: the total number of majorized problems that have been solved, i.e., the total number of calls for g02aaf.
14: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $-\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}>0$.
${\mathbf{ifail}}=2$
On entry, ${\mathbf{ldg}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{ldg}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=3$
On entry, ${\mathbf{ldx}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{ldx}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=4$
On entry, ${\mathbf{rank}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: $0<{\mathbf{rank}}\le {\mathbf{n}}$.
${\mathbf{ifail}}=5$
Majorized penalty approach fails to converge in maxit level of penalty iterations.
${\mathbf{ifail}}=6$
Convergence is limited by machine precision. The objective function value or rank is decreasing very slowly.
The array returned in x may still be of interest.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

The returned accuracy is controlled by errtol and ranktol, and is limited by machine precision.

## 8Parallelism and Performance

g02akf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02akf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

Arrays are internally allocated by g02akf. The total size of these arrays is $11×{\mathbf{n}}+3×{\mathbf{n}}×{\mathbf{n}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(2×{\mathbf{n}}×{\mathbf{n}}+6×{\mathbf{n}}+1,120+9×{\mathbf{n}}\right)$ real elements and $5×{\mathbf{n}}+3$ integer elements.

## 10Example

This example finds the nearest correlation matrix with target rank $2$ to:
 $G = 2 -1 0 0 -1 2 -1 0 0 -1 2 -1 0 0 -1 2$

### 10.1Program Text

Program Text (g02akfe.f90)

### 10.2Program Data

Program Data (g02akfe.d)

### 10.3Program Results

Program Results (g02akfe.r)

## 11Algorithmic Details

g02akf is aimed at solving the rank constrained nearest correlation matrix problem formulated as follows:
 $min X∈𝕊n 12 X-G 2 s.t. Xii =1 , i=1,…,n, X⪰0, rankX≤r,$ (1)
where $G$ is the input approximate correlation matrix, $r$ is the upper bound of rank of the output nearest correlation matrix $X$, the expression $X⪰0$ stands for a constraint on eigenvalues of $X$ in the space of $n$ by $n$ symmetric matrices ${𝕊}^{n}$, namely, all the eigenvalues should be non-negative, i.e., the matrix should be positive semidefinite, $‖·‖$ denotes the Frobenius norm. Note that the rank constraint is given as an inequality, which means if the intrinsic rank of the input matrix is already less than or equal to $r$, the solver will calculate a nearest correlation matrix without increasing the rank.
This section contains a short description of the algorithm used in g02akf which is based on the Majorized Penalty Approach (MPA) by Gao and Sun (2010). Further details on accuracy and stopping criterion are also included.

### 11.1Penalty approach

Let $\lambda \left(X\right)$ be the vector of the eigenvalues of $X$ (arranged in the non-increasing order) and $P\in {ℝ}^{n×n}$ be the corresponding matrix of orthonormal eigenvectors of $X$. The equivalent relationship for a positive semidefinite matrix $X$ is as follows.
 $rankX≤r⇔λr+1X+⋯+λnX=0.$
Therefore, problem (1) can be equivalently written as
 $minX∈𝕊n fX≔12X-G2 s.t. Xii=1, i=1,…n, X⪰0, λr+1X+⋯+λnX=0.$ (2)
Introducing the penalty parameter $c>0$, we can obtain the following penalized problem by taking a trade-off between the rank constraint and the least square distance $f\left(X\right)$.
 $minX∈𝕊n fX + cλr+1X+⋯+λnX s.t. Xii=1, i=1,…,n, X⪰0.$ (3)
Define
 $pX ≔ ∑ i=r+1 n λiX = I,X - ∑ i=1 r λiX,$
where $I$ is the identity matrix, $〈·,·〉$ is the standard trace inner product in ${𝕊}^{n}$. We can rewrite problem (3) as
 $minX∈𝕊n fX+c⁢pX s.t. Xii=1, i=1,…,n, X⪰0.$ (4)
The penalty parameter $c$ is updated according to the progress of rank reduction. The input argument maxit controls the maximum number of updates of $c$.
The penalized problem (4) is not equivalent to the original problem (1) and the relationship can be described as follows. If the rank of the minimizer ${X}_{c}^{*}$ to problem (4) is not larger than $r$, then ${X}_{c}^{*}$ is a global optimal solution to problem (1), otherwise an $\epsilon$-optimal solution to problem (1) is guaranteed given that parameter $c$ satisfies some bound constraint. Please see Gao and Sun (2010) for more details.

### 11.2Majorization approach

The focus now is on solving the penalty problem (4). Since $p\left(X\right)$ is nonsmooth and concave, we majorize it by the linear function defined by its subgradient. For given ${X}^{k}$ (the current iteration) and ${U}^{k}\in \partial p\left({X}^{k}\right)$, we have
 $pX ≤ mkpX ≔ pXk + Uk,X-Xk ∀X.$
Now, instead of solving the nonconvex problem (4), we solve the following convex model:
 $minX∈𝕊n fX + c⁢mkpX s.t. Xii = 1 , i=1,…,n, X⪰0.$ (5)
The framework combines ideas of a penalty method with majorization, it can be described as follows:
 Majorized Penalty Algorithm (MPA) 1.Select a penalty parameter $c>0$ and a feasible point ${X}^{0}$, set $k≔0$. 2.Solve subproblem (5) to get ${X}^{k+1}={\mathrm{argmin}}_{X}\left\{f\left(X\right)+c{m}_{k}^{p}\left(X\right)|{X}_{ii}=1\text{, ​}i=1,\dots ,n\text{, ​}X⪰0\right\}$. 3.If ${X}^{k+1}={X}^{k}$, stop; otherwise, update penalty parameter $c$, set $k≔k+1$ and go to step 2.
Let $\overline{G}=G-c{U}^{k}$, the subproblem (5) is actually a nearest correlation matrix problem with input $\overline{G}$ without rank constraint, which can be solved efficiently by g02aaf. Input maxits controls the maximum number of iterations used in solving one problem (5) with fixed $c$.
For the convergence result of the algorithm, see Gao and Sun (2010). For other implementations as well as parameter setting strategies, see Qi and Yuan (2014) and Bai et al. (2015).

### 11.3Stopping criteria

The algorithm shown in Table 1 is stopped when all the stopping criteria are satisfied to the requested accuracy, these are:
 $2fXk - 2fXk-1 max100, 2fXk-1 ≤ ε1, (relative precision)$
 $∑ i=r+1 n λiXk ≤ ε2. (rank feasibility)$
Here ${\epsilon }_{1}$ and ${\epsilon }_{2}$ may be set by arguments errtol and ranktol respectively in order to achieve various returned accuracy. The above quantity to measure rank feasibility does not scale well with the magnitude of ${X}^{k}$. To rectify this drawback, we also build in the third stopping criterion to control the percentage of the first $r$ eigenvalues of ${X}^{k}$ out of all the eigenvalues:
 $∑ i=1 r λiXk / ∑ i=1 n λiXk ≥ 90% .$