# Minimizing expression over symmetric matrices

I'm trying to solve the following maximization problem over space of symmetric matrices A and positive definite H

$$R=\max_{A\in S(R^d)}\frac{\text{tr}(HA)^2+2\text{tr}(HAHA)}{\text{tr}(AHA)}$$

So far I've used regular Maximize over matrix coefficients, but it only works for small matrices. Is there a way to massage this problem into some other optimization tool that's more matrix specific?

Clear[a];

d = 2;
A = Array[a[Min[#1, #2], Max[#1, #2]] &, {d, d}];
vars = DeleteDuplicates[Flatten[A]];
H = DiagonalMatrix[{1., 2.}];
objective = (Tr[H.A]^2 + 2 Tr[H.A.H.A])/Tr[A.H.A];
Maximize[objective, vars] // First (* 6.56155 *)
$$$$


That's an eigenvalue problem! How to see that? Well, let's define some example data:

d = 3;
H = N@DiagonalMatrix[Range[d]];
A = # + #\[Transpose] &@RandomReal[{-1, 1}, {d, d}];


Both numerator and denominator can be written as bilinear forms of the vector $$\operatorname{vec}(A)$$:

Avec = Flatten[A];


To this end, we need the following matrices and vectors:

Hvec = Flatten[H];
(*representing multiplication by H from the left*)

HL = KroneckerProduct[H, IdentityMatrix[d]];
(*representing multiplication by H from the right*)

HR = KroneckerProduct[IdentityMatrix[d], Transpose[H]];
(*representing the bilinear form of the numerator*)

M1 = 2 HR.HL + KroneckerProduct[Hvec, Hvec];
(*representing the bilinear form of the denominator*)
M2 = HL;


Testing for correctness:

Avec.M1.Avec == Tr[A.H]^2 + 2 Tr[H.A.H.A]
Avec.M2.Avec == Tr[A.H.A]


True

True

Thus your optimatization problem is equivalent to maximizing a Rayleigh quotient:

$$\max_{A \in S(\mathbb{R}^d)} \frac{\operatorname{vec}(A)^\top \, M_1 \, \operatorname{vec}(A)}{\operatorname{vec}(A)^\top \, M_2 \, \operatorname{vec}(A)}$$

Pulling back the bilinear forms M1 and M2 to the subspace of symmetric matrices is fairly easy:

idx = Partition[Range[d^2], d];
p = DeleteCases[Flatten[UpperTriangularize[idx]], 0];
q = DeleteCases[Flatten[UpperTriangularize[Transpose[idx]]], 0];
L = SparseArray[
Transpose[{Join[p, q],
Join[Range[Length[p]], Range[Length[p]]]}] -> 1.,
{d^2, Length[p]}, 0.
];
B1 = L\[Transpose].M1.L;
B2 = L\[Transpose].M2.L;


Test:

Avec[[p]].B1.Avec[[p]] == Avec.M1.Avec
Avec[[p]].B2.Avec[[p]] == Avec.M2.Avec


So we arrive at the following maximization problem

$$\max_{v \in \mathbb{R}^{d (d+1)/2)}} \frac{v^\top \, B_1 \, v}{v^\top \, B_2 \, v}$$

and this an generalized eigenproblem for the matrix pair $$(B_1,B_2)$$. In fact, we are looking for the greatest generalized eigenvalue and that can be obtain as follows:

Eigenvalues[{B1, B2}, 1]


30.0574

However, this method has super high numerical complexity. The matrices $$B_1$$ and $$B_2$$ are of size $$n \times n$$ where $$n = \frac{d (d+1)}{2}$$. And most generalized eigenvalue algorithms require to invert $$B_2$$ (or to solve linear system with it by utilizing a $$LU$$-factorization). This takes $$O(n^3)$$ time, so the overall complexity is $$O(d^6)$$. However, if $$H$$ is diagonal, then $$B_2$$ is also diagonal and the inversion of $$B_1$$ has complexity only $$O(n) = O(d^2)$$.

I think one can bring it down to $$O(d^3)$$ for generic symmetric positive-definite $$H$$ by diagonalizing $$H$$ first (before building $$B_1$$ and $$B_2$$. In fact, I believe that one can arrange everything so that $$B_2$$ is the identity matrix; then we would have reduced the problem to a classical eigenvalue problem for a single matrix $$B$$. And that can be solved very efficiently with the power method, requiring only matrix-vector multiplications with $$B$$. And since $$B$$ has quite a particular structure (a rank one matrix + HR.HL), there is quite a lot of potential to make this faster.

• For computing the eigenvalues of a symmetric-definite pencil, you indeed can convert it into the usual symmetric eigenproblem by performing the Cholesky decomposition of B2; as it is diagonal in this case, one can just do With[{dd = DiagonalMatrix[1/Sqrt[Diagonal[B2]]]}, Eigenvalues[dd.B1.dd, 1]]. (The general procedure is described here, or see this for some Mathematica code.) – J. M. is in limbo Dec 1 '19 at 15:06
• Thanks for in-depth answer! Actually the problem is rotationally invariant and H is covariance matrix so I can always replace H with a diagonal – Yaroslav Bulatov Dec 1 '19 at 22:14
• You're welcome! And as J.M. mentioned, one can indeed use a Cholesky-factorization of H` to recast everything into a classical eigenvalue problem. And since only the highest eigenvalue is required, the power method should work fine. The first paper cited by J.M. might be really worth a read. – Henrik Schumacher Dec 1 '19 at 22:31