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

• 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. Commented Feb 15, 2021 at 23:09
• Perhaps I'm being slower than usual, but what is the norm you're using on $\otimes^N\Bbb R^d$? Commented Feb 15, 2021 at 23:19
• @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$ Commented Feb 16, 2021 at 8:00
• 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. Commented Feb 16, 2021 at 13:38
• @TheSimpliFire: Indeed, typo corrected, thanks Commented Feb 16, 2021 at 14:31

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$$

$$\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];

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"]