User gilead - MathOverflow

Answer by Gilead for Avoiding epsilon in mixed integer linear and quadratically constrained programs

2013-03-08T00:32:58Z

If you're solving this problem numerically (using floating point numbers), $x > 0$ really means $x \geq \epsilon$ where $\epsilon$ is the next representable floating-point number on the number line for the given precision, in the increasing direction.

If you think about it, floating point numbers are essentially a discrete approximation to real numbers, therefore strict inequalities can only be approximated. Therefore in numerical optimization, one cannot directly formulate constraints with $\gt, \lt$ or $\neq$. There is always an $\epsilon$ involved.

As for $\epsilon$ not being numerically nice, in theory yes. But in practice that rarely happens if one is careful to avoid catastrophic cancellations. Also in many applications $\epsilon$ is chosen to be the numerical tolerance (which is usually many orders of magnitude larger than the machine epsilon).

Answer by Gilead for If then condition on mixed linear integer programming

2012-08-14T20:57:52Z

Yes.

$a_{\text{max}} y \leq a \leq (1-y)(a_{\text{max}} - \epsilon) + a_{\text{max}}y$

$0 \leq b \leq b_{\text{max}} y$

where $y \in \{{0,1\}}$ and $\epsilon$ is a small positive real.

Answer by Gilead for Quadratic problem solving with absolute value constraint

2011-09-26T16:33:18Z

Well, you have a non-smooth unconstrained problem but there exists a standard reformulation of the absolute function into linear constraints.

$ \begin{align} &\max_{x} -\frac{1}{2}x^T A x + b^T x - C \sum_{i} z_{i}\\ s.t.\quad & z_{i} = s_{i}^{+} + s_{i}^{-}, \quad \forall i\\ & x_{i} = s_{i}^{+} - s_{i}^{-}, \quad \forall i\\ & s_{i}^{+}, s_{i}^{-} \geq 0, \quad \forall i\\ \end{align} $

with $A \succ 0$.

Mind you, this formulation may not give you the correct results if you decide to add constraints on $x_{i}$; if you are constraining $x_{i}$, you may need to reformulate this into a Mixed Integer Quadratic Program (MIQP).

Answer by Gilead for Another question about formulating absolute value constraint in optimization (Sorry for the mistake I made in the previous posting)

2011-01-14T06:37:01Z

This is really a question for http://www.or-exchange.com (since the answer requires practical know-how rather than mathematical abillity). However, since linear programming questions are of some interest to some folks here, I'll make an attempt at an answer that is not completely useless.

There are many ways of formulating an absolute function in an LP (some bad, some good). I'm going to discuss the bad ways, in case you are tempted to use them.

The bad ways

You can try these approaches, but bear in mind that they are not rigorous unless the absolute function is the sole term in the objective (I believe the $l_{1}$-norm LP in compressed sensing fulfills this criterion).

As far as I am aware, there is no 100% reliable way of formulating an absolute function that appears in the constraint set in an LP if there exists a competing objective $\Phi$.

The standard LP approach that is often used (but IMHO, is a poor method) is as follows: introduce a dummy variable $z$ (that is, $z=|C - D|$) and nonnegative slack variables $s_{0},s_{1}$. Write the LP as below: $$ \begin{align} &\min \Phi + z\\ s.t.\;\; & z = s_{0} + s_{1}\\ & C - D = s_{0} - s_{1}\\ & s_{0} \geq 0, s_{1} \geq 0\\ & A - B = C - D\\ &A + B = z \end{align} $$ where $\Phi$ is the original objective. In theory, this seems like it will work, but in practice, depending on how the other constraints are posed, and the "downward pressure" of $z$ with respect to $\Phi$, this might not always give you the correct answer. For instance, if $\Phi$ makes the absolute function $z$ tend toward a non-minimum value, it will depend on the weighting between $z$ and $\Phi$ to determine which term "wins".
A similar (but equally flawed) approach is to use the fact that $|C - D| = \max(C-D,D-C)$, and to write this: $$ \begin{align} &\min \Phi + z\\ s.t.\;\; & z \geq D-C \\ & z \geq C-D\\ & A - B = C - D\\ &A + B = z \end{align} $$

The good ways (but your problem will no longer remain an LP)

The only reliable way to formulate an absolute function in the constraint set is to reformulate your LP as MIP (Mixed Integer Program).

If your solver supports indicator constraints, you can write the following: $$ \begin{align} &\min \Phi\\ s.t.\;\; & z \geq D - C\\ & z \geq C - D\\ & z \leq D - C\text{ or }z\leq C - D \\ & A - B = C - D\\ &A + B = z \end{align} $$ where the "or" is handled as an indicator constraint. Most solvers will use a Big-M formulation to convert the problem into an MIP.
The best way is to use a mixed-integer (MIP) formulation. In order to do that, you need to know the upper bound for $C \in [0,C^{U}]$. (Since $D$ is a known, we'll assume it is constant.) If you have no idea what the upper bound is, choose an adequately large value for $C^U$, bearing in mind that very large values of $C^{U}$ can cause conditioning problems. Also, the larger the $C^{U}$, the poorer your LP-relaxation for branching will be, which in turn will adversely impact the performance of the solution process. So choose $C^{U}$ carefully. First, define an upper-bound $U$ as follows, $U = \max(D,C^{U})$. Then write the following constraints:

$$ \begin{align} &\min \Phi\\ s.t.\;\; & 0 \leq C \leq C^{U}\\ & 0 \leq z - (C - D) \leq (2U)\delta_{1}\\ & 0 \leq z - (D - C) \leq (2U)\delta_{2}\\ & \delta_{1} + \delta_{2} = 1\\ & A - B = C - D\\ &A + B = z \end{align} $$ where $\delta_{1},\delta_{2} \in \{0,1\}$ (binary variables).

Answer by Gilead for Sampling arbitrary point from system of linear inequalities

2010-11-15T15:55:10Z

I agree with Boris -- you want to solve an LP. If you would like to solve it numerically, just pose it like this: $$ \begin{align} & \min 0 \\ s.t. \quad & Ax + e\leq b \end{align} $$ where $e \in \mathbb{R}^m$ is a vector containing $\epsilon$'s, which is some numerical tolerance.

Optimal knot placement for fitting piecewise-continuous linear functions to a nonlinear function

2010-11-12T00:56:18Z

I encountered this problem in my research and it is turning out to be a surprisingly difficult one(for me, at least).

Suppose we have a univariate nonlinear function $f(x)$ where $x \in [L,U]$. Our goal is to approximate this nonlinear function with $n$ piecewise-continuous linear functions $g_{i}(x)$ within the given domain. We assume that $n$ is a pre-specified number. We define each line segment as follows: $$ g_{i}(x) = \frac{f(a_{i}) - f(a_{i-1})}{a_{i} - a_{i-1}} (x - a_{i-1}) + f(a_{i-1})\text{ for }a_{i-1} \leq x \leq a_{i} $$ where $a_{i}$ are knot points in $[L,U]$ and $i = 1,\ldots,n$. The first and the last knot points are fixed at the boundaries, that is, $a_{0} = L, a_{n} = U$. Also, the knot points are ordered and unique: $ a_{i} > a_{i-1}$ for $i=1,\ldots,n$.

I want to find the optimal placements for the knot points $a_{1},\ldots,a_{n-1}$, such that the overall squared-approximation error $e$ is minimized. We can pose the objective as follows: $$ \min_{a_{1},\ldots,a_{n-1}} \left\{ e = \int_{L}^{U} [f(x) - g_{i}(x)]^2 dx \right\} $$

This picture illustrates the problem:

The final optimization problem looks like the following (after a simple reformulation into a optimal-control-like form): $$ \begin{align*} &\min_{a_{1},\ldots,a_{n-1}} e(U)\\ s.t.\quad & \frac{de(x)}{dx} = [f(x) - g_{i}(x)]^2, \quad e(L) = 0\\ &g_{i}(x) = \frac{f(a_{i}) - f(a_{i-1})}{a_{i} - a_{i-1}} (x - a_{i-1}) + f(a_{i-1})\text{ for }a_{i-1} \leq x \leq a_{i}\\ & a_{0} = L, a_{n} = U\\ & a_{i} \geq a_{i-1} + \epsilon,\quad i=1,\ldots,n \end{align*} $$ This optimization problem is extremely difficult to solve numerically, owing to its nonsmoothness and nonconvexity.

Question: How do I solve this problem to global optimality? Can anyone provide any attacks (even partial ones)? Any simplifying properties?

Answer by Gilead for One-Variable Optimization Problem

2010-11-08T15:20:20Z

Hmm, you could rewrite it this way: (I'm going to assume that $\pi_0, \pi_1, \xi'', \mu, \lambda$ are pre-defined constants, and $\alpha,W$ are variables)

$$ \begin{align} &\max_{W}\; \pi_{0} (1 - e^{-\mu W}) - \frac{\pi_{1}}{W}z_{1}(W)\\ s.t.\;& \frac{dz_{0}(\alpha)}{d\alpha} = 1 - e^{-\mu \alpha},\quad z_{0}(0) = 0\\ & \frac{dz_{1}(\alpha)}{d\alpha} = 1 - e^{-\lambda \alpha},\quad z_{1}(0) = 0\\ & z_{0}(W) + \epsilon \leq \xi''\\ & W \geq 0 \end{align} $$ where $z_{0},z_{1}$ are auxiliary variables and $\epsilon$ is a numerical tolerance value. This then becomes a DAE (differential-algebraic equation) optimization problem, which can be solved numerically. (though given that your decision variable is also the independent variable in the differential equations, some further bilinear transformations may be required. See http://dx.doi.org/10.1016/j.na.2005.03.066. Some software packages do this automatically.).

Edit: I just realized, if indeed $\xi''$ is a constant as I have assumed, the inequality constraint can easily be converted into a bound. $$ \begin{align} \int_{0}^{W^U} 1-e^{-\mu \alpha}\,d\alpha = \xi'' - \epsilon\\ \frac{e^{-\mu W^{U}}}{\mu} + W^{U} - \frac{1}{\mu} = \xi'' - \epsilon\\ \end{align} $$ Solve for $W^{U}$, and replace the above inequality constraints with: $$ 0 \leq W \leq W^{U} $$ You may be able to solve this using optimal control methods.

Answer by Gilead for How do I approach Optimal Control?

2010-10-29T03:43:13Z

My field is mathematical programming, and I tend to look at optimal control as just optimization with ODEs in the constraint set; that is, it is the optimization of dynamic systems. I would start by studying some optimization theory (not LPs but NLPs) and getting an intuitive feel for the motivations behind stationarity and optimality conditions -- that will lead naturally into optimal control theory.

I should mention there is another facet of optimal control, related to control systems. The systems considered are discrete time (as opposed to continuous in PMP) therefore it's difference equations instead of differential equations. Examples of optimal control laws in this latter sense are Linear Quadratic Regulators (LQRs), Linear Quadratic Gaussian (LQGs), Model Predictive Control (MPC). It is this latter type of optimal control that is actually applied in industry. The Pontryagin principle, while useful for analysis, is generally intractable for real-time application to nontrivial plants.

Answer by Gilead for Expressions for the Square of an Integral

2010-10-20T03:27:02Z

I'm not sure about simplifying, but you can easily write your objective functional in Bolza form like this:

$$ \begin{align} &\min_{u(t) \in \Omega(t)} \, J = z(T)^2 + \int_{0}^{T} s(t)u(t)dt \\ s.t. &\frac{dz(t)}{dt} = r(t)u(t),\quad z(0) = 0 \end{align} $$

Answer by Gilead for Kalman filtering: 1D case

2010-09-22T01:05:48Z

A few remarks on your problem:

You have to assume something for your initial variance (not covariance in this case, since it's univariate). The same applies in the multivariate case -- you have to know something about $P_{0|0}$. You do not calculate the initial variance.
If you really have no idea what to choose for your initial variance, choose a large number. This is equivalent to saying "I don't know what's going on in the system, so I'm going to be conservative and assume the worst." As the Kalman filter iterates, it will generally converge and the variance will tend to decrease.
Given a measurement $z_{0}$, you can do the rest (Kalman gain, prediction etc.). In fact in the linear case, it is proven that the Kalman gain can be calculated off-line (see "Separation Principle" http://en.wikipedia.org/wiki/Separation_principle).
If your filter is having trouble converging (very unlikely in this simple case), you can use something called a Re-iterative Kalman Filter (http://tinyurl.com/2fokknm). This Kalman filters iterates $n$ steps and uses the information collected to correct $x_{0}$. At $n+1$, it uses to the corrected $x_{0}$ and recursively calculates $x_{n+1}$; thereafter the Kalman filter will usually converge rapidly.

Peter D. Joseph (a pioneer in the use of Kalman Filters in the 1960s) wrote a simple tutorial on the subject in which he gives the reader an intuitive understanding of what these filters do -- in it he motivates the subject through the derivation of a 1-D example. Unfortunately the webpage no longer exists; however I managed to find the original document in text format: http://www.humintel.com/hajek/kalman.txt

If you're willing to reformat it into $\LaTeX$, I think you'll find the document helpful.

Answer by Gilead for Applications of Math: Theory vs. Practice

2010-09-11T07:29:25Z

My work draws on various bodies of mathematics. Here's a brief description (by no means exhaustive) of how I use math in my work:

Mathematical programming/optimization: Optimization is used for anything from reconciling actual data to a model's (regression), estimating unknown parameters in a system, to finding the best inputs that will extremize some functional in a dynamic system. The applications are endless. When people think mathematical programming, they think Linear Programming. But convex nonlinear programming is actually very well-established. In fact, large problems in nonconvex optimization are routinely solved (although modeling a nonconvex system can be quite an art).
Real/functional analysis: useful for understanding optimization algorithms. An understanding of convex functions and sets is crucial -- they lead to global solutions (with guarantees) without solving an NP-hard problem, so we exploit convexity properties whenever possible. (Lipschitz) continuity is another important idea, subgradients etc. are important concepts in nonsmooth optimization. Real analysis is not applied directly, but a good understanding of it is required for reading convergence proofs or descriptions of optimization algorithms.
Computational Geometry: ideas like convex hulls, Voronoi diagrams, etc. are useful in optimization. I use them to partition a problem space into convex regions, or to parametrize a space. The region bounded by convex polytopes can be represented by a set of inequality constraints that can be enforced in an optimization problem. Discrete optimization is used to optimally switch between these regions.
ODE/DAE theory: used for modeling dynamic systems. In particular, understanding the notion of index in DAEs can help one develop models that are amenable to reliable numerical solution.
Calculus. Differential calculus is used everywhere (e.g. model sensitivity analysis, automatic differentiation, postoptimality analysis)
Statistics: projection methods like the Karhunen-Loewe transform (related to SVD) are used to reduce the dimensionality of large models constructed from data. They're also the only way to handled correlated/collinear data (in practice, most large datasets in the real world are correlated. The assumption of factor independence built into standard regression techniques often does not hold, so methods like multiple linear regression often have to be modified for instance into principal components regression in order for them to be usable on real world large datasets). Also, tools like time-series analysis are used to construct time series models from data.
Linear algebra: used almost everywhere. They're the basic building blocks for working with nonlinear systems. In particular, efficient numerical solution of sparse structured matrices is crucial to the efficiency of large-scale nonlinear optimization algorithms (the bottleneck is often in the linear algebra solvers, not in the optimization algorithm itself). Tools like SVD are frequently used.
Numerical methods: used everywhere. Understanding concepts like numerical conditioning is crucial; when modeling, one wants to end up with a system with a Jacobian that is well-conditioned with respect to inversion.
Misc: Diophantine equations are used to derive certain control laws. Laplace transforms are used for modeling linear-time-invariant systems because they allow differential equations to be manipulated as algebraic ones. Algebraic Riccati equations are solved in the derivation of the Kalman gain. Fixed-point iteration is used to converge decomposed models.

Answer by Gilead for Converting an ODE system to State space formulation

2010-09-09T16:54:17Z

Ok, let me give this a shot:

Because your system is nonlinear, I'm assuming you want the nonlinear state-space form. You can easily get the linear form by doing a Taylor series expansion on it around some equilibrium point.
The fact that most of the states are not measurable is not a big problem. You can estimate them using your output variables (subject to observability conditions), using an state observer such as a Kalman filter, Moving Horizon Estimator (MHE) or a Luenberger observer. Also note for a nonlinear system, only local observability can be checked.
Because you have conditional statements, I don't believe you'll be able to write the above as a single state-space system. You have 4 conditions, but they can be reduced to 3 disjunctions, so you'll need 3 state-space systems and a conditional switching equation that "activates" the correct state-space system depending on the values of $v$ and $y$. This is known as a hybrid (or switched) system.
As to the treatment of a hybrid system, perhaps you could clarify what the purpose is of getting your model into state-space form. Is it for simulation reasons? Do you need it in order to do analysis (i.e. controllability, observability)? Or do you need to do optimization? If it is the last case, you can write the logic as a disjunctive program, which will allow you generate a very efficient mixed-integer programming (MIP) problem.

Anyway, this is one way of writing your state-space system:

$$ \begin{align} \frac{dy_{i}}{dt} &= g_{i}^{m}(v_{i},y_{i}) + p_{i} c_{i} + q_{i} e_{i},\quad i=0,\ldots,N-1\\ v_{i} &= \sum_{j=0}^{N-1} a_{ij}y_{j},\quad i=0,\ldots,N-1 \end{align} $$ where $p_{i},q_{i} \in \{0,1\}$ = coefficients, $N$ = number of states, and $m \in \{1,2,3\}$ = modes of the system. In addition, you will need a switching function $T(m,v_{i},y_{i}) = 0$ to select the appropriate modes based on the current states of your system. This can be done programmatically through IF-THEN-ELSE clauses (or via integer variables in an optimization problem).

For $m = 1$ (where $v_{i} < 0, y_{i} < \varepsilon$):

$$g_{i}^{1}(v_{i},y_{i}) = v_{i} \left[1-\exp\left(\frac{{v_i}^2S_i}{{v_i}(\varepsilon -y_{i})^2}\right)\right]$$

For $m = 2$ (where $v_{i} \leq 0, y_{i} \geq \varepsilon$):

$$g_{i}^{2}(v_{i},y_{i}) = v_{i}$$

For $m = 3$ (where $v_{i} > 0, y_{i} \in \mathbb{R}$):

$$g_{i}^{3}(v_{i},y_{i}) = \frac{{v_i}}{1+\frac{{v_i}}{S_i}(1-\exp(-{v_i}/S_i))}$$

Answer by Gilead for Continuous optimization

2010-08-19T15:26:24Z

This is how I would develop the formulation (conceptually).

Investment costs, $Q$: $$Q = \sum_{i} m_{i} c_{i}$$
Investment returns for $i$ at time $t$: we know that the return for $i$ is $m_{i} r_{i}$ if $t \geq t_{i}$, otherwise it is 0. To model this if logical condition: $$ R = \sum_{i} \delta_{i} m_{i} r_{i} $$ $$ \delta_{i} = 1, \text{ if } t \geq t_{i} $$ $$ \delta_{i} = 0, \text{ if } t < t_{i} $$ where $R$ = overall returns. In this case, $t$ and $t_{i}$ are parameters, therefore $\delta_{i}$ are parameters too, not binary variables. They can be pre-calculated for a specified $t$.

In summary, your problem can be represented as follows:

$$ \max_{m_{i}} (R - Q)$$ s.t. $$ Q = \sum_{i} m_{i} c_{i}$$ $$ R = \sum_{i} \delta_{i} m_{i} r_{i}$$ $$ Q \leq C $$ $$ m_{i} \geq 0\quad \forall i$$

where $C, c_{i}, r_{i}, \delta_{i}, t, t_{i}$ are parameters.

So, if $m_{i} \in \mathbb{R}$, this would be a linear program. If $m_{i} \in \mathbb{N}$, this becomes a mixed-integer linear program. (Note: for tractability in the integer case, you may need to specify a reasonably small upper-bound for $m_{i}$ or use partial integer variables)

For $t \rightarrow \infty$ (i.e. all the investments have reached maturity), simply set all $\delta_{i} = 1$.

Answer by Gilead for Efficient computation of AB^-1 for matrices

2010-08-01T21:19:45Z

What are the sizes of $\mathbf{A}$ and $\mathbf{B}$? This information is important.

Let me assume you mean you want the efficient numerical computation of the matrix $ \mathbf{A} \mathbf{B}^{-1}$.

The general strategy would be to do this: Let $\mathbf{J} = \mathbf{A} \mathbf{B}^{-1}$; therefore $\mathbf{J}\mathbf{B} = \mathbf{A}$. You must then rewrite this into $\mathbf{P}\mathbf{x} = \mathbf{Q}$ form (this depends on the dimensions of your matrices).

For instance, for $\mathbf{A} \in \mathbb{R}^{m \times n}$, $\mathbf{B} \in \mathbb{R}^{n \times n}$ and $\mathbf{J} \in \mathbb{R}^{m \times n}$, you can write a system of equations:

$J_{i,*} \cdot B_{*,j} = A_{i,j}$

(Notation: for a matrix $\mathbf{X}$, we define $X_{i,*}$ as the $i$-th row vector and $X_{*,j}$ as the $j$-th column vector, and $X_{i,j}$ as the element in the $i$-th row and $j$-th column).

From here, you can use a fast linear solver to solve the resulting linear equation system -- you will get a solution for the elements of $\mathbf{J}$.

By solving a linear system of equations and not taking the inverse directly, you're not only cutting down on the no. of operations required, you can exploit also properties like sparsity, inertia, etc. and have capabilities like pre-conditioning at your disposal.

Answer by Gilead for Continuous Linear Programming: Estimating a Solution

2010-07-25T07:12:05Z

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. ;-)

Answer by Gilead for Minimizing a function containing an integral

2010-07-23T20:06:18Z

It looks like a mixed-integer dynamic optimization problem. Your problem can be rewritten as follows: (notice the transformation of the integral into a differential equation? It's a standard trick. Also, note that you need an initial condition for $Y$)

$\min_{x(t)} L(T)$

s.t. $ \frac{dL(t)}{dt} = AR(t)-x(t)$, with $L(0) = 0$

$ \frac{dR(t)}{dt} = ax(t)R(t)Y(t) - bR(t)$, with $R(0) = R_{0}$

$\frac{dY(t)}{dt}=−x(t)R(t)Y(t)$, with $Y(0) = Y_{0}$

$x(t) = \delta(t) x_{min} + (1 - \delta(t)) x_{max}$ where $\delta(t) \in {0,1}$

To solve this problem numerically, simply discretize the differential equations using backward Euler (easy), or implicit Runge Kutta (harder, but more accurate). Pose this as a Mixed Integer Nonlinear Program (MINLP) and use one of these solvers to find the solution.

Bonmin https://projects.coin-or.org/Bonmin

Couenne https://projects.coin-or.org/Couenne

These solvers will traverse the branch-and-bound tree more intelligently and efficiently than your method of enumerating every single case, which will grow with the no. of discretization grid points you have. (e.g. let's say you discretize over 20 points; the no. of cases you have to search is $2^{20} = 1048576$. Not nice.)

With a branch-and-bound|cut|reduce MINLP solver (and a bit of luck), on average you are unlikely to hit the worst case scenario where every single case is enumerated.

There are other ways of solving this problem -- multiple-shooting, sequential dynamic optimization, etc. In my opinion, optimal control methods (Pontryagin's maximum principle) are typically intractable on problems like this.

Answer by Gilead for Minimization under non-linear constraints

2010-07-23T19:32:40Z

As I understand, your problem looks like this:

$\min_{x,y} \Phi=a_{1} x + a_{2}y$

s.t. $f(x,y) = 0$

where $f(x,y)$ looks something like this: $f(x,y) = a_{3} \exp{(a_{4}x)} + a_{5} \exp{(a_{6}y)}$

This looks like a nonconvex NLP can be trivially solved using any NLP solver.

Or are you looking for a closed form solution?

Normally, the first thing I would try is to see if I can substitute constraints into the objective to transform the problem into an unconstrained one, but it looks like it's not possible here.

As mentioned above, you can write the first order optimality conditions for the above system and solve the resulting nonlinear system of equations.

$\nabla L(x,y,\lambda) = \nabla\Phi + \lambda \nabla f(x,y) = 0$

In your case, it would be:

$a_{1} + \lambda a_{3} a_{4} \exp{(a_{4}x)} = 0$

$a_{2} + \lambda a_{5} a_{6} \exp{(a_{6}x)} = 0$

$a_{3} \exp{(a_{4}x)} + a_{5} \exp{(a_{6}y)} = 0$

Solve the above system for $x,y,\lambda$. And bam! You're done.

Answer by Gilead for Is all non-convex optimization heuristic?

2010-07-23T00:00:22Z

Hi there, I'm coming to this from a practitioner's point of view. Your question as to whether non-convex optimization is always heuristically driven can be answered as follows:

No.

There are many gradient-based techniques for nonconvex global optimization out there that DO NOT rely on any heuristics at all. They are usually based on partitioning the solution space, and performing some sort of branch and bound search using tight convex relaxations (the tightest relaxations obtainable for nonlinear functions are McCormick relaxations). As mentioned, these algorithms have a worst case exponential complexity, but they are rigorous (non-heuristic) and are able to give you a provably global solution.

Global nonconvex optimization is an active area of research: http://www.mat.univie.ac.at/~neum/glopt/techniques.html#branch

The well-known BARON software for instance, can rigorously find the global optimum of a nonlinear nonconvex problem.

Other software/algorithms include:

LaGO https://projects.coin-or.org/LaGO

Couenne https://projects.coin-or.org/Couenne

Papers that provide mathematical details for all the above solvers can be found in the open literature.

Practitioners in this area have come to realize that procedures for finding the global solution to a general nonconvex problem are usually NP-hard (so far no exceptions have been found).

A special case of this can be seen in polynomial programming, where a nonconvex polynomial optimization problem can be solved by decomposing its KKT (optimality) conditions into its Groebner basis. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.6266

On the surface, this looks attractive because it would seem that any nonconvex optimization problem can then be approximated as a polynomial programming problem by taking its Taylor expansion. However, the computation of Groebner bases is NP-hard.

I hope that gives you a few leads.

Comment by Gilead

2013-05-01T17:09:27Z

p.s. KKT conditions are undefined for open-sets. However if $\Omega + \lambda I \succeq 0$, it would seem that the infimum is 0, attained at $r = 0$.

Comment by Gilead

2013-05-01T17:00:58Z

Is $\Omega + \lambda I$ positive semi-definite? Are you interested in an approximate numerical solution, or a closed-form expression?

Comment by Gilead

2013-03-31T19:28:29Z

It looks like a special-case of a symmetric matrix, but I'm not sure if there's a specific name for it. You might want to go through the list here: en.wikipedia.org/wiki/List_of_matrices

Comment by Gilead

2013-03-10T14:26:18Z

Does VNODE do what you want? cas.mcmaster.ca/~nedialk/Software/VNODE/… Or VSPODE? www3.nd.edu/~markst/…

Comment by Gilead

2013-03-04T16:39:21Z

It looks like a nonlinear fractional program (en.wikipedia.org/wiki/Fractional_programming). You can solve it as a general nonlinear program, but there may be specific properties you can exploit if you treat it as a fractional program.

Comment by Gilead

2013-01-18T00:18:12Z

This sounds like a homework problem.

Comment by Gilead

2012-11-10T21:16:12Z

Do you mean "else"? What does "or" mean in this case? Also, try this forum: or-exchange.com for formulation questions instead of MathOverflow.

Comment by Gilead

2012-08-15T05:30:26Z

In general, you cannot write strict inequality constraints (< or >) in numerical math programs. The $\epsilon$ is used as an approximation, and is usually set to your solver's numerical tolerance.

Comment by Gilead

2012-07-18T20:45:14Z

Why do you want to avoid linear programming?

Comment by Gilead

2012-07-17T21:31:28Z

You can either solve this to a local optima using nonlinear programming (fast), solve to global optima using global optimization code (potentially computationally expensive), or solve a separable programming problem using piecewise-linear approximation of the [possibly] nonconvex objective function using SOS2 constraints--which is designed for problems like this--or simply using MIP constraints. Your method of minimizing $x$'s one-by-one (essentially solving $n$ univariate problems) may not give you a global solution unless a global solver is used in each instance.

Comment by Gilead

2012-06-22T14:56:32Z

See here: amsterdamoptimization.com/pdf/lineq.pdf

Comment by Gilead

2012-06-22T14:56:18Z

As it turns out, using an LP solver to solve linear equations isn't a terrible idea in practice. This may sound counter-intuitive, because it sounds overkill. However, because the efficiency of an LP solution depends so much on the linear algebra, good commercial LP solvers like CPLEX or Gurobi tend to incorporate highly optimized linear algebra routines (with excellent sparse matrix handling). There are other advantages to using an LP solvers: 1. matrix does not have to be square; 2. detection of inconsistencies through slacks; 3. inequality constrained solutions for underdefined systems.

Comment by Gilead

2012-06-22T14:51:51Z

Is $m > n$ or vice-versa?

Comment by Gilead

2012-04-29T20:07:35Z

Since yours is an operations research questions, you may want to ask your question here: or-exchange.com.

Comment by Gilead

2012-03-12T22:50:09Z

@Brian, in a comment, he mentioned that the objective $f$ was a sum of hinge losses. en.wikipedia.org/wiki/Hinge_loss