How to maximize the product of scalar products?

Question

Given $N$ unit vectors $v_1,\ldots,v_N$ of $\mathbb R^d$, I was curious about finding an explicit maximizer for the product mapping $$ f(x) := \prod_{j=1}^N \langle x,v_j\rangle,\qquad x\in\mathbb R^d, $$ assuming that, say, $\|x\|=1$.

If one replaces the product by a sum, then the Cauchy-Schwarz inequality would directly yield that the maximum is reached at the unit vector proportional to the sum $\sum_jv_j$, but for the product I was not able to find such a simple solution.

I've tried to use the tensor product formalism to write $f(x) = \langle x^{\otimes N},\otimes_{j=1}^N v_j\rangle$, so that we see that maximizing $f$ boils down to minimize $$ g(x):=\|x^{\otimes N}-\otimes_{j=1}^N v_j\|^2, $$ and thus the solution is the first tensor component of the "orthogonal projection" of $\otimes_{j=1}^N v_j$ onto the subset $\Delta := \{x^{\otimes N}: x\in\mathbb R^d\}\subset(\mathbb R^d)^{\otimes N}$. Unfortunately, $\Delta$ is not a vector space and is not even convex, hence the quotation marks. I don't know how to continue (and gradient derivations did not bring me anywhere). Any ideas?

I'm pretty sure that clever people have worked on how to approximate a tensor product by a tensor product of the same vector, but I guess I didn't give the good keywords to Google.

Answer 1

Try Lagrange multipliers. As $f(x)=\prod_{j=1}^N\langle x,v_k\rangle$ making

$$ L = f(x)+\lambda (\|x\|^2-1) $$

we have

$$ \nabla L = f(x)\left(\sum_{k=1}^N\frac{v_k}{\langle x, v_k\rangle}\right)+2\lambda x=0 $$

now calling $\mu = \frac{2\lambda}{f(x)}$ we have the stationary points as the solutions for

$$ \cases{ \sum_{k=1}^N\frac{v_k}{\langle x, v_k\rangle}+\mu x=0\\ \|x\|^2=1 } $$

or $\mu = -N$

so we follow with

$$ \sum_{k=1}^N\frac{v_k}{\langle x, v_k\rangle}-N x=0 $$

a set of $d$ equations with $d$ unknowns.

NOTE

This nonlinear system can be solved iteratively as follows

$$ x_{k+1}= x_k -\left(\nabla g(x_k)\right)^{-1}g(x_k) $$

where

$$ \cases{ g(x) = \sum_{k=1}^N\frac{v_k}{\langle x, v_k\rangle}-N x\\ \nabla g(x) = -\left(\sum_{k=1}^N\frac{v_k\otimes v_k}{\langle x,v_k\rangle^2}+N I_d\right) } $$

NOTE

Attached a MATHEMATICA script to perform the calculations. The iterative process is fired maxtries because the stationary points can be many.

Clear[dif, grad, iG, Y, Y1]
SeedRandom[1]
dim = 10;
n = 10;
v = RandomReal[{-1, 1}, {n, dim}];
V = Table[v[[k]]/Norm[v[[k]]], {k, 1, n}];
X = Array[x, dim];
obj = Product[(V.X)[[k]], {k, 1, n}];

dif[X_List] := Sum[V[[k]]/(V[[k]].X), {k, 1, n}] - n X;
grad[X_List] := -Sum[Outer[Times, V[[k]], V[[k]]]/(X.V[[k]])^2, {k, 1, n}] - n IdentityMatrix[dim];
iG[X_List] := Quiet@Inverse[grad[X]]

kmax = 20;
maxtries = 1000;
SOLS = {};

For[try = 1, try <= maxtries, try++,
  Y = RandomReal[{-1, 1}, dim];
  aerror = 10^-20;

  For[k = 1, k <= kmax, k++,
    Y1 = Y - iG[Y].dif[Y];
    If[Max[Abs[Y - Y1]] < aerror, AppendTo[SOLS, Y1]; Break[]];
    Y = Y1
  ]
]

precise = 10^15;
sols = N[Union[Round[precise SOLS] / precise]];
best = Sort[Table[{obj, sols[[k]]} /. Thread[X -> sols[[k]]], {k, 1, Length[sols]}]][[Length[sols]]]

We can compare it with the result obtained from a nonlinear solver

solX = NMaximize[{obj, X.X == 1}, X, Method -> "DifferentialEvolution"]
obj /. Thread[X -> best[[2]]]