11
$\begingroup$

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 *)
```
$\endgroup$
13
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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.) $\endgroup$ – J. M. is in limbo Dec 1 '19 at 15:06
  • $\begingroup$ 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 $\endgroup$ – Yaroslav Bulatov Dec 1 '19 at 22:14
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 1 '19 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.