Suppose we have a set of rank deficient covariance matrices. How can I know the similarities and differences between those set of matrices?
MathOverflow is a question and answer site for professional mathematicians. It only takes a minute to sign up.Sign up to join this community
Per the comments, the objective function has real(trace(.)), not trace(.). I will not bother with dividing by N (also mentioned in comments), since that doesn't affect the argmin.
This can be formulated as a non-convex Quadratically-Constrained Quadratic Programming (non-convex QCQP) problem, which is much more difficult to solve than a convex QCQP. I will illustrate formulation and solution by numerical nonlinear optimization using YALMIP under MATLAB. In particular, I will solve the problem for a single R, as provided by @Rawan in the comments. This can be extended to multiple $R_k$, as discussed.
R=[1.00000000000000 + 0.00000000000000i,0.00743491184002355 + 0.963010458121173i,-0.856972159053805 + 0.0118693873167640i;0.00743491184002355 - 0.963010458121173i,1.00000000000000 + 0.00000000000000i,0.00743491184002355 + 0.963010458121173i;-0.856972159053805 - 0.0118693873167640i,0.00743491184002355 - 0.963010458121173i,1.00000000000000 + 0.00000000000000i] R = 1.000000000000000 + 0.000000000000000i 0.007434911840024 + 0.963010458121173i -0.856972159053805 + 0.011869387316764i 0.007434911840024 - 0.963010458121173i 1.000000000000000 + 0.000000000000000i 0.007434911840024 + 0.963010458121173i -0.856972159053805 - 0.011869387316764i 0.007434911840024 - 0.963010458121173i 1.000000000000000 + 0.000000000000000i
In order to accommodate the inverse in the objective function using YALMIP, I will replace the inverse term in each summand in the objective function by $Y_k$, and add a constraint forcing it to be the requisite inverse. Namely $Y_k(X^ HR_k X + I) = I$.
Xinit is formed as the normalized eigenvectors of R in order of largest to smallest eigenvalue per the OP. And can be extended to multiple $R_k$ per OP. Xinit is used as a starting value for the numerical optimization. A corresponding Yinit_k is calculated as the requisite matrix inverse corresponding to Xinit for the kth term. This is also part of the starting value for the nonlinear optimization, given that Yinit_k is an optimization variable. Y and Yinit will be used in the example below, given K = 1.
Xinit = -0.570058480786439 + 0.008229635484230i -0.706935175940842 + 0.015200995175087i -0.418052075226221 + 0.014479321306010i 0.004269613754453 + 0.591534422975448i -0.004813610768336 + 0.000051746765831i -0.013956113277691 - 0.806133270749517i 0.570117881158577 + 0.000000000000000i -0.707098588060266 + 0.000000000000000i 0.418302747237492 + 0.000000000000000i Yinit = 0.259316197658287 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i -0.000000000000000 + 0.000000000000000i 0.000000000000000 - 0.000000000000000i 0.874928251150862 - 0.000000000000000i 0.000000000000000 - 0.000000000000000i -0.000000000000000 - 0.000000000000000i 0.000000000000000 + 0.000000000000000i 0.999247369410125 + 0.000000000000000i
The value of the objective function at the starting value in this example is 0.8665. The optimization below reduces it to 7.5634e-04 with (locally) optimal X =
1.0e+03 * -1.029044383364627 + 0.014855774377308i -0.381043583110040 + 0.008193455157269i -0.000431665630392 + 0.000014950829646i 0.007707318110758 + 1.067811805360625i -0.002594573816645 + 0.000027891910840i -0.000014410583735 - 0.000832384401651i 1.029151610292002 + 0.000000000000217i -0.381131663518623 - 0.000000000011879i 0.000431924465344 + 0.000000000000002i
In the YALMIP optimize command, the first argument is the constraints, and the second argument is the objective function.
N=3; k=3; % problem dimensions X = sdpvar(N,k,'full','complex'); % declare X as a full complex N by k matrix assign(X,Xinit); % assign Xinit as starting value for X Y = sdpvar(N,k,'full','complex'); % declare Y as a full complex N by k matrix assign(Y,Yinit); % assign Yinit as starting value for Y optimize(Y*(X'*R*X+eye(N)) == eye(N),real(trace(R-R*X*Y*X'*R)),sdpsettings('usex0',1))
value(X) now has the (locally) optimal value of X, presuming the solver successfully solved the problem.
On larger problems, this could be computationally demanding.
Note that the problem is non-convex. Therefore, there may be, and generally are for this problem, local minima which are not globally optimal. So good starting value for X will help a local optimizer a lot. Using a descent method, the solution obtained should be no worse than, and may be much better than, the starting value. If a global optimizer is used, this may not be relevant, but a good starting value may speed things up. A (rigorous) global optimizer will likely be intractable for desired size versions of this problem.
Edit: In response to a new comment by the OP, the optimization can be changed to add a new constraint trace(X*X') <= upper_bound as follows: Use
optimize([Y*(X'*R*X+eye(N)) == eye(N),trace(X*X') <= upper_bound],real(trace(R-R*X*Y*X'*R)),sdpsettings('usex0',1))
where upper_bound can be whatever you want it to be.