How to maximize the product of scalar products? 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.
 A: 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]]]

