Given a matrix $A$, I'm interested in finding the vectors $x,y$ maximizing $$x^T A y$$ where $\sum_i x_i^p = \sum_i y_i^q = 1$. Here $p,q\in[0,\infty)$ are some given values.
This is similar to the normal singular vector problem in the case $p=q=2$.
I considered the Lagrangian $x^T A y - \lambda_1 \|x\|_p^p - \lambda_2 \|x\|_p^p$, which gives the equations:
$$\begin{align} A y &= \lambda_1 p x^{p-1},\\ x^T A &= \lambda_2 q y^{q-1}. \end{align}$$
This suggests an iterative algorithm, similar to the power iteration method for singular values, where we repeatedly compute $x_{n+1} = (A y_n)^{1/(p-1)}$ and $y_{n+1} = (x_n^T A)^{1/(q-1)}$ and normalize.
This algorithm appears experimentally to work when $p$ and $q$ are sufficiently large. ($p,q\ge 1$ might not be quite enough.)
I wonder if anyone knows any literature on computing said values. Are there any efficient algorithms for computing $x$ and $y$? Does my algorithm work for some range of $p,q$? I'm particularly interested in the case where $A$ is a stochastic matrix.
Update: I found http://www.njohnston.ca/2016/01/how-to-compute-hard-to-compute-matrix-norms/ which suggests an algorithm very similar to the above, except that we alternate in updating $x$ and $y$, rather than doing it simultaneously. For some reason, that change makes a big difference, and I now get near-immediate convergence for all the random $A$s I try.
It seems that the algorithm (Boyd's algorithm) is known to converge to the optimal value whenever $p,q\ge 1$: https://arxiv.org/pdf/1001.2613.pdf