How do I solve $\min_x \max(c_1^Tx, c_2^Tx, \dots, c_k^Tx)$ for $\lVert x \rVert_2 = 1$.

Question

Let $f(x) = \max(c_1^Tx, c_2^Tx, \dots, c_k^Tx)$. where $x, c_1, c_2, \dots, c_k \in \mathbb R^n$. What fast iterative methods are available for finding the (approximate) min of $f$ with the constraint $\lVert x \rVert_2 = 1$?

Notes:

1. $f$ is convex and and non negative, that is, $c_i$ positively span $\mathbb R^n$
2. For my use case $k \gg n$ and (rougly) $3 \le n \le 50$ and $100 \le k \le 10000$. I tried projected subgradients (projected to the sphere) but it can be slow to converge and improving the initial guess doesn't seem to accelerate the method. From the literature, this seems to be equivalent to Riemannian manifold subgradient methods which use the exponential map. I haven't tried bundle methods yet but I am investigating them now. They seem like they might be too slow.

I've also tried approximating $f$ with a smooth maximum LogSumExp and doing gradient descent with projection. I also tried unconstrained gradient descent with a quadratic penalty. The approximation is too inaccurate and the evaluation of exp is too slow for my use case. I'm not interested in SQP because I suspect whatever method used to smooth (e.g. LogSumExp, hyperbolic) will be too costly/inaccurate to evaluate.

Edit: I believe this problem is equivalent to a linear programming problem with quadratic equality constraints. I haven't been able to to find much literature on this type of problem. Perhaps I can use a SDP method but I'm not sure.

Answer 1

EDIT: I have included all of the MATLAB code here: https://gitlab.com/TheRTGuy/point-cloud-support-function-minimization. There I also have other evaluations.

Given the assumption that $f$ is always nonnegative, then this implies that the convex hull of $C=\{c_1,\ldots,c_k\}$ contains the origin. Were this not the case, then $f$ is always negative, and the origin is separable from the point cloud. In this case, minimizing $f$ just amounts to finding the supporting hyperplane that maximizes the distance from the origin, a problem that is solved easily using convex optimization (e.g. solving the dual of the primal problem).

The problem, as some pointed out in the comments, boils down to enumerating the facets of the convex hull of $C$, and finding the one closest to the origin. This is a very hard problem, see https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node22.html#polytope:Vredundancy, since the number of facets is exp. in $n$ and $k$. So, forget about an "efficient" exact solution.

However, in my research I needed a way to compute a tight overapproximation of the convex hull of a point cloud. Here is a simplified adaptation, of the approach I used, to OP's problem:

The approach:

1. Generate vertices of a regular simplex, $\{h_1,\ldots,h_{n+1}\}$, and compute $$b_0 = \begin{pmatrix}\max_{c\in C} h_1^\top c \\ \vdots \\ \max_{c\in C}h_{n+1}^\top c\end{pmatrix}.$$
2. Form the new simplex $S_0 = \{x\in\mathbb{R}^n~\vert~ H_0x \leq b_0\} = \operatorname{Conv}(\{s^1_0,\ldots,s^{n+1}_0\})$, where $$H_0 = \begin{pmatrix}h_1^\top\\\vdots \\h_{n+1}^\top\end{pmatrix},$$ and $s^i_0$ is the $i$-th vertex of the simplex. This simplex tightly bounds the set $C$.
3. Optionally, perturb the simplex by a random rotation.
4. Set $i=1$, and $f_0 = \min_j [b_0]_j$, where $[\cdot]_j$ is the $j$-th component of a vector.
5. For each vertex $s^j_{i-1}$ solve the LP $$ x_j=\arg\min_x ~ g(x)= x^\top s^j_{i-1} \quad \\s.t. \\\quad (c_l - s^j_{i-1})^\top x \leq -1, ~~~\forall l\in\{1,\ldots,k\}, \\ \quad (s^l_{i-1}- s^j_{i-1})^\top x \leq -1, ~~~\forall l\in \{1,\ldots,n+1\}\setminus\{j\},\\ \quad x^\top s^j_{i-1} \geq -1. $$ The LP basically finds a supporting hyperplane of $C$ that separates it from $s^j_{i-1}$. This plane is often a facet of $\operatorname{Conv}(C)$.
6. Set $$ H_i = \begin{pmatrix}\frac{x_1^\top}{\lVert x_1\rVert}\\\vdots\\\frac{x_{n+1}^\top}{\lVert x_1\rVert}\end{pmatrix},\qquad b_i = \begin{pmatrix}\frac{x_1^\top s^1_{i-1}-1}{\lVert x\rVert} \\ \vdots \\ \frac{x_{n+1}^\top s^{n+1}_{i-1}-1}{\lVert x\rVert} \end{pmatrix}. $$ This forms a new simplex $S_i = \{x\in\mathbb{R}^n~\vert~ H_ix \leq b_i\} = \operatorname{Conv}(\{s^1_i,\ldots,s^{n+1}_i\})$.
7. Set $f_i = \min\{\min_j[b_i]_j , f_{i-1}\}$, $i:=i+1$, and repeat steps 5-7 for $N$ iterations, or until $f_i = f_{i-1}$ for some consecutive iterations.
8. Set set the optimal value to $f^* = f_i$.

Evaluation: The implementation I have is for MATLAB.

1. I evaluated the approach for $n=2,\ldots,20$, with $k=1000$. For each $n$ I ran the algorithm 3 times on a uniformly generated set $C$, to which I also apply a random rotation, and time each run. The number of iterations is fixed to $N=5$. In case you are curious, this is the set of MATLAB commands executed to generate the point cloud:

%Randomly rotated uniform
C = 5*rand(n,N)-2.5;
%Indexing pattern for off-diagonal entries of matrix
idx = nchoosek(1:n,2);
dU = zeros(n);
dth = pi*(2*rand(n*(n-1)/2,1) - 1); %Uniform
dU(sub2ind([n,n],idx(:,1),idx(:,2))) = dth;
dU(sub2ind([n,n],idx(:,2),idx(:,1))) = -dth;
C = expm(dU) * C;

2. To evaluate how good the solution is, for each run I generate 1000000 points, $\{x_k\}$, normally distributed on the unit sphere in $\mathbb{R}^n$, and evaluate $f$ at each point. Then, I compute $d = \min_k\{f(x_k)\} - f^*$ and plot it. Positive means that the solution is good, negative means the solution is bad. Of course, if positive, this does not mean the solution is the best. Furthermore, the higher the dimension, the lower the quality of this evaluation. Ideally, we want to compare it with another algorithm, which I don't have.

Results:

The first plot shows the execution times vs. $n$. The second plot shows the difference $d$ vs $n$. The final plot is just an arbitrary run for $n=20$, which shows $f(x_k)$ and $f^*$.

$n=20$" />

Conclusion:

On first glance, the approach seems to perform very well. Only a small number of iteratios is necessary to get a good suboptimal solution. Also, solving LPs is very efficient, even for large $n$ and $k$. However, there is no rigorous tight bound for optimality yet. More rigorous analysis is needed. Furthermore, evaluating the approach with normally distributed points on the unit sphere is obviously not the best way, but at least it gives a sense of optimality. This is clearly seen from the second plot where the gap is increasing. If the OP is interested in further evaluation, then I would be willing to share some code.

The intuition is simple: we try to "squeeze" the points by simplices as much as possible, by using the vertices of previous simplices as anchors. The facets of the simplices often coincide with the faces of $\operatorname{Conv}(\{c_i\})$.

Answer 2

Convert it to a quadratically constrained linear program by the following trick: Add an auxiliary variable $t$ such that: $\max(c_1^Tx, c_2^Tx, \dots, c_k^Tx)\leq t$ Then, this will imply $c_i^Tx\leq t~~ \forall i$.

Now your optimization program can be rewritten as: $$ \text{minimize}_{x, t} ~~t\\ \quad\text{subject to}\\ \quad\|x\|_2= 1\\ \quad \quad \quad \quad \quad \quad \quad c_i^Tx\leq t~~ \forall i\in\{1, \ldots, k\} (*) $$

Up to this point the optimization program is exactly the same as the original one. From a convex optimization perspective relaxing the norm constraint results in an efficiently solvable problem.

If you relax $\|x\|_2= 1$ to $\|x\|_2\leq 1$, then the program is convex and can be efficiently solved using any convex optimization technic such as interior point method. The solution might satisfy the original constraint $\|x\|_2 = 1$ by itself. But there is no guarantee for that as that program is not convex unfortunately. However, you can perform Step two, and project back $x$ to $\|x\|_2 = 1$ by setting $x = x / \|x\|_2$

If the solution turns out to be zero, you should solve Problem (*) using iterative projection method, where in each iteration you project back to $\|x\|_2 = 1$. This will give you an approximate solution to the non-convex problem you have.

Answer 3

We need to solve the following optimization problem (P$0$): \begin{align} (\mathrm{P0})\qquad &\min_{x} \max\{c_1^\mathsf{T}x, \cdots, c_k^\mathsf{T}x\}\\ &\mathrm{s.t.}\quad \|x\|_2 = 1. \end{align}

Assume that $\{c_1, c_2, \cdots, c_k\}$ satisfies $\max\{c_1^\mathsf{T}x, \cdots, c_k^\mathsf{T}x\} > 0$ on $\mathbb{R}^n\backslash \{0\}$.

(P$0$) is equivalent to $\min_{x} \frac{\max\{c_1^\mathsf{T}x, \cdots, c_k^\mathsf{T}x\}}{\|x\|_2}$ and $\max_{x} \frac{\|x\|_2}{\max\{c_1^\mathsf{T}x, \cdots, c_k^\mathsf{T}x\}}$ and \begin{align} &\max_{x} \|x\|_2^2\\ &\mathrm{s.t.}\quad \max\{c_1^\mathsf{T}x, \cdots, c_k^\mathsf{T}x\} = 1 \end{align} and \begin{align} (\mathrm{P1})\qquad &\max_{x} \|x\|_2^2\\ &\mathrm{s.t.}\quad c_i^\mathsf{T}x \le 1, \ i=1, 2, \cdots, k. \end{align} (P$1$) is a Quadratic Programming (QP) problem (https://en.wikipedia.org/wiki/Quadratic_programming).

Some online and software Resources for QP are given in [1]. Some quadratic programming codes are given in [2].

References

[1] https://neos-guide.org/content/quadratic-programming

[2] http://www.numerical.rl.ac.uk/people/nimg/qp/qp.html