Minimize trace of inverse of convex combination of matrices.

Hello! (First question--please forgive me if its unclear.)

I am interested in efficient/approximate optimization techniques for minimizing a norm of a convex combination of symmetric, positive semi-definite matrices. In particular I am wondering what techniques exist to efficiently solve the following problem:

$\alpha^{*}=\arg\min\big\{\operatorname{trace}((\sum_i\alpha_i \mathcal{I}_i)^{-1}):\sum_i\alpha_i=1, \alpha_i\ge0 \big\}$

where $\mathcal{I}_i$ are positive semi-definite matrices and $\operatorname{trace}$ is the sum of the eigenvalues. Edit: Furthermore, assume $\mathcal{I}_i\succ0$ for at least one $i$.

Approximations are acceptable--I am interested in any "standard" approaches, should they exist. Should some other scalar function be more convenient this would be interesting to me too (ie determinant, nuclear norm, spectral radius, etc). The matrices are often sparse and could be extremely large.

The simplest idea I had was to lower bound the sum of matrices using Weyl's theorem, ie, $\lambda_i(A+B)\ge \lambda_j(A)+\lambda_{n+i-j}(B)$ for $j\ge i$ giving me an upper bound to the trace of the inverse of the sum. I am wondering if another bound would be tighter or better still some clever variational approach exists.

The motivation (shouldn't be relevant to the answer) for this problem comes from statistical machine learning. I would like to develop techniques for characterizing maximum composite likelihood estimators (a generalization of the maximum likelihood estimator). Classical statistical analysis yields a convex combination of the expected Hessian matrices and is as above. (The $\mathcal{I}$ denotes the Fisher information matrix.)

There is some similarity to the Eigenvalues of Matrix Sums question although I believe this question is sufficiently different to ask anew.

Thanks,

Josh

• Have you tried the method of Lagrange multipliers? What happens if/when you do it for this problem? (I assume you have the sum of at most n choose 2 symmetric matrices and that the norm is continuous on the range of interest.) Gerhard "Email Me About System Design" Paseman, 2011.06.20 – Gerhard Paseman Jun 20 '11 at 22:44
• In addition to $\alpha^T1=1$, you need $\alpha_i \ge 0$ for the combination to be convex. Also, if all the $A_i=0$, then there is no solution; or do you assume the $A_i > 0$? – Suvrit Jun 21 '11 at 0:31
• @Suvrit Sorry for the ambiguity; indeed $alpha_i\ge 0$. I can only assume that one of the $A_i$ is positive definite. – jvdillon Jun 21 '11 at 13:49

Define,

$$A(\alpha) = \sum_i \alpha_i A_i.$$

In your case (ignoring constraints for the moment), you have

$$\min\quad\mbox{trace}(A(\alpha)^{-1}) = \sum_i e_i^TA(\alpha)^{-1}e_i,$$

which makes it somewhat tricky to optimize.

So, as per your request, here is a somewhat simpler setup that you may find useful. Consider,

Now, look at the following two problems:

$$\min_{\alpha}\quad \|A(\alpha)\|_2,$$

and

$$\min_{\alpha}\quad c^TA(\alpha)^{-1}c.$$

Both of these can be cast into SDPs. For the first, we replace it by $$\begin{eqnarray} \min_{\alpha,t}\quad &t\\\\ &\begin{bmatrix} tI & A(\alpha)\\\\ A(\alpha) & tI \end{bmatrix} \succeq 0 \end{eqnarray}$$

For the second, use a similar Schur-complement based idea to obtain $$\begin{eqnarray} \min_{\alpha,t}\quad &t\\\\ &\begin{bmatrix} A(\alpha) & c\\\\ c^T & t \end{bmatrix} \succeq 0 \end{eqnarray}$$

You can add more constraints on $\alpha$ if you want. However, because these things are SDPs, they won't be too easy to solve for really large-matrices.

• you'll need to do this for $c=e_{1}$, $c=e_{2}$, ..., so you'll end up with a block diagonal SDP constraint. – Brian Borchers Jun 21 '11 at 5:38
• Hi Suvrit--thanks for taking the time to think about this. First question: Does minimizing the L2 norm make sense? Perhaps you mean to maximize this? Anyway, that's an easy enough fix. Second question: Regarding SDP approaches: I have toyed with SDP solvers (and similar formulations as you have given; I also tried logdet) for toy models but it is far too unscalable to be useful. I would much prefer an approximate solution. Do you know of techniques to approximate SDP problems? (I am imagining a majorization/minimization procedure.) – jvdillon Jun 21 '11 at 14:01
• If you just want to approximate the solution that can be done easily---one way is, as you mention, MM; more simply, you can just lower-bound the objective, and then solve a linear problem. However, I did not mention that, because it leads to a "trivial" solution. As for "minimization", since the leading singular value (here eigenvalue) is a convex function, yes, you do want to minimize, not maximize. – Suvrit Jun 21 '11 at 21:50

Since the functional is convex, I would try modified gradient descent. The main problem with this approach is that you'll really have to compute the inverse matrix a few times, which may be out of question due to the "extremely large size" condition. If so, just disregard this answer.

The algorithm is fairly simple: if you have $m$ matrices $A_j$ start with $a_j=1/m$ for all $j$. Compute $A(a)=\sum_j a_j A_j$. Compute $b_j=\operatorname{trace} A_jA^{-2}$. Now trace the path in the simplex $a_j\ge 0,\sum_j a_j=1$ that goes from $a$ in the direction of $b$ until it reaches some vertex (so you start moving in the direction of the projection of $b$ to the plane $\sum_j x_j=0$ until you reach a face, then move along the face in the direction that is the projection of $b$ to that face and so on). Let $a(t)$ be the point on that path at distance $t$ (along the path) and $T$ be the path length (you can use the $\ell^2$ or the $\ell^1$ distance, it doesn't seem to matter). Now start computing $\operatorname{trace}A(a(2^{-k}T))^{-1}$ ($k\ge 0$) until $2^{-k}T<\varepsilon$. If the minimal one you attained is less than $\operatorname{trace}A(a)^{-1}$, move $a$ to the corresponding point and repeat.

To keep track of how close you are to the optimal value, just compare $\max_j b_j$ and $\operatorname{trace}A(a)^{-1}$. Their difference dominates the error.

I tried it with thirty 100 by 100 random matrices and it converged pretty fast but, as I said, if you are in no position to attempt a single matrix inversion or multiplication, this is pretty useless.

It is perhaps useful to restate your problem as the obviously equivalent problem of minimizing $-\gamma_1/\gamma_0$ where $x^n+\sum_{i=0}^{n-1}\gamma_ix^i$ is the characteristic polynomial of a matrix in the set $\{\sum \alpha_i\mathcal I_i,\sum \alpha_i=1\}$. Usual numerical methods applied to the rational function (in terms of parameters) $-\gamma_1/\gamma_0$ work perhaps well.

It is perhaps also useful to have information on local maxima of $\det(A)=\pm \alpha_0$ for $A$ in your set. Solutions for your problem are perhaps not too far away.

• Hi Rolanda--thanks for answering. Unfortunately I think the matrices are too complicated to use the characteristic polynomial. – jvdillon Jun 21 '11 at 19:12