2
$\begingroup$

Suppose we have a set of rank deficient covariance matrices. How can I know the similarities and differences between those set of matrices?

Regards,

Improve this question
$\endgroup$
3
  • $\begingroup$ You could do numerical nonlinear optimization of MSE w.r.t. X, using your eignevector-based X as starting value. if that is fairly near to the optimum, it might not be too bad. If you "need to" get rid of the inverse, you can replace $(X^{H}R_{k}X+{I})^{-1}$ by $Y_k$ with $Y_k$ being another matrix variable for each k, and adding the constraints $Y_k(X^{H}R_{k}X+{I}) = I$. How large is N? $\endgroup$
    – Mark L. Stone
    Commented Mar 23, 2018 at 1:29
  • $\begingroup$ I was suggesting using your existing approach to find a value of X which can be used as a starting value for a numerical nonlinear optimization. If there are extra requirements on the properties or structure of X beyond being an N by N REAL matrix which minimizes MSE, you need to say what they are. Are all matrices actually real, not complex? If so, why use $X^H$ rather than $X^T$? N=500 will be much harder to solve than X=50. Either use numerical differentiation or automatic differentiation (perhaps, ADiMat in order to handle matrix calculations) or Matrix Cookbook. $\endgroup$
    – Mark L. Stone
    Commented Mar 23, 2018 at 11:42
  • $\begingroup$ Are covariance matrices complex, not real? If so, then MSE is complex not real, so doesn't make sense as an objective function. $\endgroup$
    – Mark L. Stone
    Commented Mar 23, 2018 at 12:24

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ You can find a single X which optimizes over several $R_k$. You need to sum objective function terms, such as I showed, for each $R_k$. The YALMIP solution I showed will require a separate $Y_k$ matrix variable and matrix constraint for each $R_k$ as I mentioned. You specified a problem which is not so easy, so unless you have a very clever method, the solution may require a lot of computation. You don't necessarily have to solve to optimality if you are content with an approximate solution. You can get improvement in the objective function on each iteration once feasibiility is attained. $\endgroup$
    – Mark L. Stone
    Commented Mar 23, 2018 at 20:53
  • $\begingroup$ I selected X, call it Xinit, as you described. But you can make it whatever you want. In my answer I wrote "A corresponding Yinit_k is calculated as the requisite matrix inverse corresponding to Xinit for the kth term." More explicitly, that is $Yinit_k = (X^ HR_k Xinit + I)^{-1}$, which makes it satisfy the constraint $Yinit_k(X^ HR_k Xinit + I) = I$, and therefore be feasible with respect to it. $\endgroup$
    – Mark L. Stone
    Commented Mar 23, 2018 at 21:08
  • $\begingroup$ Note, I answered the original question. The question has now been edited to bear no relation to the question I answered. $\endgroup$
    – Mark L. Stone
    Commented Aug 30, 2018 at 14:32

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .