Minimizing expression over symmetric matrices

Question

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 *)
```

Answer 1

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.