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.

  • $\begingroup$ The most obvious suggestion, I think, is to maximize $\ln f(x)=\sum_{j=1}^N \ln \langle x,v_j\rangle$ instead. That doesn't trivialize it but does make it additive. $\endgroup$ Commented Feb 15, 2021 at 23:09
  • $\begingroup$ Perhaps I'm being slower than usual, but what is the norm you're using on $\otimes^N\Bbb R^d$? $\endgroup$ Commented Feb 15, 2021 at 23:19
  • $\begingroup$ @TedShifrin : I had in mind the Euclidean one, coming with the usual scalar product extension $\langle x\otimes u, y\otimes v\rangle_{E\otimes F}=\langle x,y\rangle_E\langle u,v\rangle_F$ $\endgroup$
    – Student
    Commented Feb 16, 2021 at 8:00
  • 1
    $\begingroup$ I think you mean $x\in\Bbb R^d$ in your first equation as $\langle x,v_j\rangle$ is only defined for equal dimensions. $\endgroup$
    – TheSimpliFire
    Commented Feb 16, 2021 at 13:38
  • $\begingroup$ @TheSimpliFire: Indeed, typo corrected, thanks $\endgroup$
    – Student
    Commented Feb 16, 2021 at 14:31

1 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.


This nonlinear system can be solved iteratively as follows

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


$$ \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) } $$


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]
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]]]
  • $\begingroup$ This is indeed a natural way to approximate a solution, thanks, but I fear this is numerically intractable even for small dimensions: I gave a quick try on Python, but the inverse differential matrix seems quickly singular for general initial conditions. $\endgroup$
    – Student
    Commented Feb 17, 2021 at 11:34
  • $\begingroup$ I will prepare ASAP to include, a MATHEMATICA script to cope with the numerical calculations. $\endgroup$
    – Cesareo
    Commented Feb 17, 2021 at 12:33
  • 1
    $\begingroup$ @Student Included the script. $\endgroup$
    – Cesareo
    Commented Feb 17, 2021 at 14:13

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .