NAG CPP Interface
nagcpp::correg::corrmat_nearest_rank (g02ak)

Settings help

CPP Name Style:



1 Purpose

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

2 Specification

#include "g02/nagcpp_g02ak.hpp"
template <typename G, typename X>

void function corrmat_nearest_rank(G &&g, const types::f77_integer rank, X &&x, double &f, double &rankerr, types::f77_integer &nsub, OptionalG02AK opt)
template <typename G, typename X>

void function corrmat_nearest_rank(G &&g, const types::f77_integer rank, X &&x, double &f, double &rankerr, types::f77_integer &nsub)

3 Description

corrmat_nearest_rank 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 (no CPP interface). 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.

4 References

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

5 Arguments

1: g(n,n) double array Input/Output
On entry: G~, the initial matrix.
On exit: a symmetric matrix G=12(G~+G~T) with diagonal elements set to 1.0.
2: rank types::f77_integer Input
On entry: r, the upper bound for the rank of X.
Constraint: 0<rankn.
3: x(n,n) double array Output
On exit: X, the nearest correlation matrix of rank r.
4: f double Output
On exit: the difference between X and G given by 12 X-G F 2 .
5: rankerr double Output
On exit: the rank error of X, defined as i=r+1 n (λi), given that λi denote eigenvalues of X sorted in non-increasing order.
6: nsub types::f77_integer Output
On exit: the total number of majorized problems that have been solved, i.e., the total number of calls for g02aaf (no CPP interface).
7: opt OptionalG02AK Input/Output
Optional parameter container, derived from Optional.
Container for:
errtoldouble
This optional parameter may be set using the method OptionalG02AK::errtol and accessed via OptionalG02AK::get_errtol.
Default: 0.0
On entry: the termination tolerance for the convergence measure of the objective function value.
If errtol0.0, then a value of 1.0E−5 is used. See Section 11 for further details.
ranktoldouble
This optional parameter may be set using the method OptionalG02AK::ranktol and accessed via OptionalG02AK::get_ranktol.
Default: 0.0
On entry: the feasibility tolerance for the rank constraint.
If ranktol0.0, machine precision is used. See Section 11 for further details.
maxitstypes::f77_integer
This optional parameter may be set using the method OptionalG02AK::maxits and accessed via OptionalG02AK::get_maxits.
Default: 0
On entry: specifies the maximum number of iterations used for the majorization approach to solve penalized problems at each level of penalty parameter.
If maxits0, then a value of 100 is used.
maxittypes::f77_integer
This optional parameter may be set using the method OptionalG02AK::maxit and accessed via OptionalG02AK::get_maxit.
Default: 0
On entry: specifies the maximum number of iterations for the penalty method, i.e., the maximum level of penalty parameter.
If maxit0, then a value of 200 is used.

5.1Additional Quantities

1: n
The order of the matrix G

6 Exceptions and Warnings

Errors or warnings detected by the function:
All errors and warnings have an associated numeric error code field, errorid, stored either as a member of the thrown exception object (see errorid), or as a member of opt.ifail, depending on how errors and warnings are being handled (see Error Handling for more details).
Raises: ErrorException
errorid=1
On entry, n = value.
Constraint: n > 0.
errorid=4
On entry, rank = value and n = value.
Constraint: 0<rankn.
errorid=5
Majorized penalty approach fails to converge in maxit level of penalty
iterations.
errorid=10601
On entry, argument value must be a value x value array.
Supplied argument has value dimensions.
errorid=10601
On entry, argument value must be a value x value array.
Supplied argument was a value x value array.
errorid=10601
On entry, argument value must be a value x value array.
Not all of the sizes for the supplied array could be ascertained.
errorid=10602
On entry, the raw data component of value is null.
errorid=10603
On entry, unable to ascertain a value for value.
errorid=10604
On entry, the data in value is stored in value Major Order.
The data was expected to be in value Major Order.
errorid=-99
An unexpected error has been triggered by this routine.
errorid=-399
Your licence key may have expired or may not have been installed correctly.
errorid=-999
Dynamic memory allocation failed.
Raises: WarningException
errorid=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.

7 Accuracy

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

8 Parallelism and Performance

Please see the description for the underlying computational routine in this section of the FL Interface documentation.

9 Further Comments

Arrays are internally allocated by corrmat_nearest_rank. The total size of these arrays is 11×n+3×n×n+max(2×n×n+6×n+1,120+9×n) real elements and 5×n+3 integer elements.

10 Example

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.1 Example Program

Source FileDataResults
ex_g02ak.cppNoneex_g02ak.r

11 Algorithmic Details

corrmat_nearest_rank 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, X0, rank(X)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 X0 stands for a constraint on eigenvalues of X in the space of n×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 corrmat_nearest_rank 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.1 Penalty approach

Let λ(X) be the vector of the eigenvalues of X (arranged in the non-increasing order) and Pn×n be the corresponding matrix of orthonormal eigenvectors of X. The equivalent relationship for a positive semidefinite matrix X is as follows.
rank(X)rλr+1(X)++λn(X)=0.  
Therefore, problem (1) can be equivalently written as
minX𝕊n f(X)12X-G2 s.t. Xii=1,  i=1,n, X0, λr+1(X)++λn(X)=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(X).
minX𝕊n f(X) + c(λr+1(X)++λn(X)) s.t. Xii=1,  i=1,,n, X0. (3)
Define
p(X) i=r+1 n λi(X) = I,X - i=1 r λi(X),  
where I is the identity matrix, ·,· is the standard trace inner product in 𝕊n. We can rewrite problem (3) as
minX𝕊n f(X)+cp(X) s.t. Xii=1,  i=1,,n, X0. (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 Xc* to problem (4) is not larger than r, then Xc* is a global optimal solution to problem (1), otherwise an ε-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.2 Majorization approach

The focus now is on solving the penalty problem (4). Since p(X) is nonsmooth and concave, we majorize it by the linear function defined by its subgradient. For given Xk (the current iteration) and Ukp(Xk), we have
p(X) mkp(X) p(Xk) + Uk,X-Xk X.  
Now, instead of solving the nonconvex problem (4), we solve the following convex model:
minX𝕊n f(X) + cmkp(X) s.t. Xii = 1 ,   i=1,,n, X0. (5)
The framework combines ideas of a penalty method with majorization, it can be described as follows:
Majorized Penalty Algorithm (MPA)
  1. 1.Select a penalty parameter c>0 and a feasible point X0, set k0.
  2. 2.Solve subproblem (5) to get Xk+1 = argminX { f(X)+ cmkp(X) | Xii=1, ​ i=1,,n, ​ X0 }.
  3. 3.If Xk+1=Xk, stop; otherwise, update penalty parameter c, set kk+1 and go to step 2.
Let G¯=G-cUk, the subproblem (5) is actually a nearest correlation matrix problem with input G¯ without rank constraint, which can be solved efficiently by g02aaf (no CPP interface). 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.3 Stopping criteria

The algorithm shown in Table 1 is stopped when all the stopping criteria are satisfied to the requested accuracy, these are:
2f(Xk)-2f(Xk-1) max{100, 2f(Xk-1) } ε1,   (relative precision)  
i=r+1 n λi(Xk) ε2.   (rank feasibility)  
Here ε1 and ε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 Xk. To rectify this drawback, we also build in the third stopping criterion to control the percentage of the first r eigenvalues of Xk out of all the eigenvalues:
i=1 r λi(Xk) / i=1 n λi(Xk) 90% .