Continuous Linear Programming: Estimating a Solution

Question

I have a "continuous" linear programming problem that involves maximizing a linear function over a curved convex space. In typical LP problems, the convex space is a polytope, but in this case the convex space is piecewise curved -- that is, it has faces, edges, and vertices, but the edges aren't straight and the faces aren't flat. Instead of being specified by a finite number of linear inequalities, I have a continuously infinite number. I'm interested in estimating solutions numerically, and my current method is to approximate the surface by a polytope, which means discretizing the continuously infinite number of constraints into a very large finite number of constraints. Unfortunately, typical linear programming algorithms run in something like cubic-time in the number of constraints, so I'm getting a huge performance hit as I make the discretization finer. Firstly, I'm interested to know if this kind of problem has been studied before, and what's been done. Secondly, I'm looking for good strategies for approaching my problem numerically (good LP packages, suggested algorithms, optimizations, etc.).

For concreteness, here is a simplified version of the problem I'm trying to solve:

I have $N$ fixed functions $f_i:[0,\infty]\to \mathbb{R}$. I want to find $x_i$ $(i=1,\dots,N)$ that minimize $\sum_{i=1}^N x_i f_i(0)$, subject to the constraints:

$\sum_{i=1}^N x_i f_i(1) = 1$, and

$\sum_{i=1}^N x_i f_i(y) \geq 0$ for all $y>2$

More succinctly, if we define the function $F(y)=\sum_{i=1}^N x_i f_i(y)$, then I want to minimize $F(0)$ subject to the condition that $F(1)=1$, and $F(y)$ is positive on the entire interval $[2,\infty)$. Note that this latter positivity condition is really an infinite number of linear constraints on the $x_i$'s, one for each $y$. A specific $y_0$ restricts me to the half-space $F(y_0) \geq 0$ in the space of $x_i$'s. As I vary $y_0$ between 2 and infinity, these half-spaces change continuously, carving out a curved convex shape. The geometry of this shape depends implicitly (and in a complicated way) on the functions $f_i$.

The reason I suspect there should be an approach that's better than just discretizing the number of constraints is that continuity of the $f_i$'s implies a kind of local structure on the space of constraints that becomes invisible under discretization. If we sit on the boundary of our convex space (so that at least N constraints are saturated, corresponding to some $y_k$), and we want to move along the boundary, then generically only those constraints corresponding to small neighborhoods of the $y_k$ are important. Sometimes when the function $F(y)$ develops a new zero, new $y$ can become important, but this is nongeneric.

NOTE: I asked this question first on stackoverflow.net, and was told it was a nonstandard enough CS problem that I should ask about it here.

Answer 1

Even though the $f_i$ are not polynomials I'll give the answer in that case because it is very nice and it seems like there is some interest. I have to stress in advance though that the answer exploits the resulting algebraic structure in a fundamental way, and so is unlikely to extend to the case when the $f_i$ are not polynomials.

First of all, a semidefinite program (SDP) is an optimization problem with matrix variables, linear objective, and positive semidefiniteness constraints on symmetric (real) matrices in addition to the standard linear (in)equalities allowed in linear programming. They are a generalization of linear programs (LP) and are vastly more expressive. LPs are the case when the matrices are constrained to be diagonal. SDP can also be viewed as a noncommutative version of LP.

The relationship with semi-infinite programming suggested by Gilead is I think the fact that one can view the constraint "A is positive semidefinite" as $x^TAx\geq 0$ for all $x$, which is an infinite number of constraints. On the other hand, one can view any convex constraint in this way, because any closed convex set can be described by (infinitely many) linear inequalities.

Theoretically SDPs are important both because many problems can be written as SDPs and because they can be solved using interior point methods in polynomial time, almost as efficiently as LPs in theory. In practice, the technology is much newer than that for LPs so one cannot solve SDPs which are nearly as big using off-the-shelf software, but those days seem to be getting closer.

To see how to turn your problem into an SDP if the $f_i$ were polynomials, let $x_i$ be your decision variables. Note that $f_i(1)$ is just a constant, so your first constraint is just a linear equation on the $x_i$ and that is no problem in an SDP. For the others, we need to do a little work.

Let $H$ denote the operator which sends a symmetric matrix to its sums along antidiagonals, so

$H: \begin{bmatrix}a & b \\\\ b & c \end{bmatrix}\mapsto\begin{bmatrix}a & 2b & c\end{bmatrix}$, and so on for bigger matrices. If we identify polynomials with their sequences of coefficients, then $p = q^2$ as polynomials if and only if $p = H(qq^T)$ as vectors. Therefore $p$ is a sum of squares (SOS) if and only if $p = H(Q)$ for a positive semidefinite $Q$ (any such $Q$ is the sum of matrices of the form $qq^T$). Now, if a polynomial $p$ is SOS then it is automatically nonnegative everywhere. Conversely, one can show than any univariate nonnegative polynomial is a sum of squares. This gives us an exact characterization of nonnegative polynomials in terms of positive semidefinite matrices. Thinking of the matrix $Q$ as a new decision variable, we can write the constraint "$p$ is nonnegative" in a semidefinite program; the solver will find a $Q$ which certifies this.

Similarly, a polynomial $p$ is nonnegative on an interval $[w,\infty)$ if and only if it is of the form $p(x) = SOS_1(x) + (x-w)\cdot SOS_2(x)$ for some SOS polynomials $SOS_i$. Again, one direction is obvious and the other requires a little effort (write $p$ in factored form and group the factors cleverly). Therefore we can also write nonnegativity on an interval in terms of positive semidefinite matrices, and hence use it as a constraint in an SDP.

Now note that when we use $H$ to define the constraint for a polynomial $p$ to be SOS, we are writing linear equality constraints between the coefficients of $p$ and some linear functionals of the matrix inside $H$. Similarly for when we write the constraint that $p$ is nonnegative on an interval: it is a linear equality between some decision variables, plus the constraint that certain matrices (the ones defining the SOS polynomials) are positive semidefinite.

Until now we've been thinking of $p$ as a constant polynomial. But because arbitrary linear equalities between decision variables are allowed in an SDP, we can just as easily write the constraint "$\sum_{i=1}^N x_i f_i(y)$ is nonnegative for $y\in[2,\infty)$" using this method, because our polynomial in question has coefficients which are linear in the decision variables.

Putting this all together gives a semidefinite program which would express exactly what you want in the case that the $f_i$ are polynomials. One could then either find a feasible $x_i$, prove that none exists, or optimize a linear functional of the $x_i$ all in polynomial time. Unfortunately because of the way we have used the algebraic structure of the problem, this is unlikely to extend to non-polynomial $f_i$.

Finally, I should note that if you're interested in this sort of thing I highly recommend checking out my advisor Pablo Parrilo's course notes on MIT OpenCourseWare. You can find the link on his website.

Answer 2

Hi again. I've been poking around, and it looks like you may be able to solve your problem as a semi-infinite program (SIP).

There is a solver caled nsips that handles SIPs. It accepts input in the AMPL modeling language. http://neos.mcs.anl.gov/neos/solvers/sio:nsips/AMPL.html

In the documentation -- http://plato.la.asu.edu/ftp/sipampl.pdf -- it says it accepts problems of this form:

$\min_{x\in\mathbb{R}^n} f(x)$

s.t. $g_{i}(x,t) \leq 0$

$h_{i}(x) \leq 0$

$h_{i}(x) = 0$

$\forall t \in T$.

However, I'm not sure if it assumes $T$ is bounded (in your case, it is not). However, even if it does, I'm sure you can define $T$ as the interval $(2,M]$ where $M$ is a very big number standing in for $\infty$.

Disclaimer: I've never worked with SIPs before and I have no idea how well this will work out for you. However, there is a good deal of theory behind it. SIPs are related to semi-definite programs; the latter can be solved very efficiently using interior-point methods. See http://www.stanford.edu/~boyd/papers/pdf/sip.pdf

Noah Stein (comment above) probably knows much much more about this. ;-)