15
$\begingroup$

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.

$\endgroup$
13
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Xander Henderson
    Feb 10, 2021 at 23:16
  • 7
    $\begingroup$ You have made more than 20 edits in the week or so since you asked this question. These kinds of frequent edits create a kind of moving target for answerers, and repeatedly bump your question to the top of the front page. Please think very carefully about what it is that you want to ask, and then type up that question and let it be. $\endgroup$
    – Xander Henderson
    Feb 10, 2021 at 23:18
  • $\begingroup$ The norm of the vectors $c_i$ plays a huge role in this problem. In fact, when they are all the same the problem is exactly solvable (it becomes a search for the smallest geodesic distance which is $\mathcal{O}(k^2)$ comparisons)- when they differ the problem becomes much more difficult to parse. If they are all comparable in size the landscape has $\mathcal{O}(k)$ minima- I think a good solution is simulated annealing, which doesn't depend on convexity, but it will require several evaluations of the starting function and an adaptive algorithm. I don't know if these suit your purposes. $\endgroup$ Feb 11, 2021 at 22:27
  • $\begingroup$ In other words, you are looking to the closest to the origin face of a convex polyhedron with up to 10000 vertices represented as the convex hull of given points in dimension up to 50 and want an exact solution in under 100 iterations of something. I'm afraid that is a bit too optimistic... $\endgroup$
    – fedja
    Feb 12, 2021 at 1:39
  • $\begingroup$ @DinosaurEgg If they are all comparable in size the landscape has $O(k)$ minima With the implicit constant independent of $n$ (or, at least growing reasonably slow in $n$)? I have my doubts about that but I may be wrong. How do you show it? $\endgroup$
    – fedja
    Feb 12, 2021 at 1:45

3 Answers 3

4
+350
$\begingroup$

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;
  1. 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^*$.

Timing

Solution quality

Example for <span class=$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\})$.

$\endgroup$
2
  • $\begingroup$ May I ask why the downvote? $\endgroup$
    – V.S.e.H.
    Feb 15, 2021 at 14:06
  • 1
    $\begingroup$ It is a nice solution (already +1 yesterday). $\endgroup$
    – River Li
    Feb 17, 2021 at 1:05
2
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ This is not the quadratic programming in the usual sense, is this? Usually, it minimizes the quadratic form with positive definite matrix. Maximization makes the solution not always unique. $\endgroup$
    – Hans
    Feb 26, 2021 at 21:24
  • $\begingroup$ @Hans The quadratic programming (QP) problem involves minimizing a quadratic function subject to linear constraints. Here it is non-convex quadratic programming problem. You can see wikipedia and the references in my answer. $\endgroup$
    – River Li
    Feb 26, 2021 at 23:16
  • $\begingroup$ All the (effective) quadratic algorithms deal with the convex kind or where the matrix is positive semidefinite. Is there an effective algorithm dealing with non-positive-semidefinite quadratic programming and finding the global optimizing point? $\endgroup$
    – Hans
    Feb 27, 2021 at 8:18
  • $\begingroup$ @Hans I think it is difficult to find global optimal solution. For the OP, $3 \le n \le 50$ and $100 \le k \le 10000$, Matlab quadprog, for example, can quickly find suboptimal solutions. $\endgroup$
    – River Li
    Feb 27, 2021 at 11:16
  • 1
    $\begingroup$ @Hans By the way, there are necessary and sufficient conditions for local minimizer of quadratic programming (not necessarily convex) in: Antal Majthay, “Optimality conditions for quadratic programming”,Mathematical Programming volume 1, pages359–365(1971). But I don't know necessary and sufficient conditions for global minimizer. $\endgroup$
    – River Li
    Feb 28, 2021 at 2:30
1
$\begingroup$

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.

$\endgroup$
8
  • 2
    $\begingroup$ Relaxing $\lVert x\rVert_2=1$ to $\lVert x \rVert_2\leq 1$ gives zero as solution if the point set is centered around the origin. $\endgroup$
    – V.S.e.H.
    Feb 13, 2021 at 15:18
  • $\begingroup$ @bodil The problem is not convex unfortunately. I'll update the answer to clarify that. $\endgroup$
    – Alt
    Feb 13, 2021 at 15:19
  • $\begingroup$ But that is what OP wants, his point cloud is centered, hence the assumption $f$ is always nonnegative. $\endgroup$
    – V.S.e.H.
    Feb 13, 2021 at 15:20
  • $\begingroup$ The point cloud is slightly off center so its not as bad as completely centered. This method gives $0$ as a solution if $f$ is completely positive $\endgroup$
    – sudowoodo
    Feb 13, 2021 at 19:08
  • 1
    $\begingroup$ @LinAlg OP's requirements are unrealistic. Even a standard LP takes at least 1-2 seconds for a 50x10000 problem, let alone a heuristic such as projected subgradients. Thus, any accelerations are out of the scope of mathematics at this point, and the problem is more related to implementation optimizations (e.g. HPS, hardware acceleration, parallel processing, etc). Not sure what initial guess OP uses, but I think unless he becomes a bit more transparent about the data and problem, then I don't think he'll get a useful answer. $\endgroup$
    – V.S.e.H.
    Feb 14, 2021 at 19:18

You must log in to answer this question.

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