Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Quadratic programming is the problem of finding a vector * x* that
minimizes a quadratic function, possibly subject to linear constraints:

$$\underset{x}{\mathrm{min}}\frac{1}{2}{x}^{T}Hx+{c}^{T}x$$

such that * A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u*.

`interior-point-convex`

`quadprog`

AlgorithmThe `interior-point-convex`

algorithm performs the
following steps:

is a sparse matrix.
For details of the sparse data type, see Sparse Matrices (MATLAB).
Generally, the algorithm is faster for large problems that have relatively
few nonzero terms when you specify H as H`sparse` . Similarly, the algorithm is faster
for small or relatively dense problems when you specify as H`full` . |

The algorithm begins by attempting to simplify the problem by removing redundancies and simplifying constraints. The tasks performed during the presolve step include:

Check if any variables have equal upper and lower bounds. If so, check for feasibility, and then fix and remove the variables.

Check if any linear inequality constraint involves just one variable. If so, check for feasibility, and change the linear constraint to a bound.

Check if any linear equality constraint involves just one variable. If so, check for feasibility, and then fix and remove the variable.

Check if any linear constraint matrix has zero rows. If so, check for feasibility, and delete the rows.

Check if the bounds and linear constraints are consistent.

Check if any variables appear only as linear terms in the objective function and do not appear in any linear constraint. If so, check for feasibility and boundedness, and fix the variables at their appropriate bounds.

Change any linear inequality constraints to linear equality constraints by adding slack variables.

If algorithm detects an infeasible or unbounded problem, it halts and issues an appropriate exit message.

The algorithm might arrive at a single feasible point, which represents the solution.

If the algorithm does not detect an infeasible or unbounded problem in the presolve step, it continues, if necessary, with the other steps. At the end, the algorithm reconstructs the original problem, undoing any presolve transformations. This final step is the postsolve step.

For details, see Gould and Toint [63].

The initial point `x0`

for the algorithm is:

Initialize

`x0`

to`ones(n,1)`

, where`n`

is the number of rows in.*H*For components that have both an upper bound

`ub`

and a lower bound`lb`

, if a component of`x0`

is not strictly inside the bounds, the component is set to`(ub + lb)/2`

.For components that have only one bound, modify the component if necessary to lie strictly inside the bound.

Take a predictor step (see Predictor-Corrector), with minor corrections for feasibility, not a full predictor-corrector step. This places the initial point closer to the

*central path*without entailing the overhead of a full predictor-corrector step. For details of the central path, see Nocedal and Wright [7], page 397.

The sparse and full interior-point-convex algorithms differ mainly in the predictor-corrector phase. The algorithms are similar, but differ in some details. For the basic algorithm description, see Mehrotra [47].

**Sparse Predictor-Corrector. **Similar to the `fmincon`

interior-point algorithm,
the sparse `interior-point-convex`

algorithm tries to find
a point where the Karush-Kuhn-Tucker (KKT) conditions
hold. For the quadratic programming problem described in Quadratic Programming Definition, these conditions are:

$$\begin{array}{c}Hx+c-{A}_{eq}^{T}y-{\overline{A}}^{T}z=0\\ \overline{A}x-\overline{b}-s=0\\ {A}_{eq}x-{b}_{eq}=0\\ {s}_{i}{z}_{i}=0,\text{}i=1,2,\mathrm{...},m\\ s\ge 0\text{\hspace{0.17em}}\\ z\ge 0.\end{array}$$

Here

$$\overline{A}$$ is the extended linear inequality matrix that includes bounds written as linear inequalities. $$\overline{b}$$ is the corresponding linear inequality vector, including bounds.

is the vector of slacks that convert inequality constraints to equalities.*s*has length*s*, the number of linear inequalities and bounds.*m*is the vector of Lagrange multipliers corresponding to*z*.*s*is the vector of Lagrange multipliers associated with the equality constraints.*y*

The algorithm first predicts a step from the Newton-Raphson
formula, then computes a corrector step. The corrector attempts to
better enforce the nonlinear constraint * s_{i}z_{i} = 0*.

Definitions for the predictor step:

, the dual residual:*r*_{d}$${r}_{d}=Hx+c-{A}_{eq}^{T}y-{\overline{A}}^{T}z.$$

, the primal equality constraint residual:*r*_{eq}$${r}_{eq}={A}_{eq}x-{b}_{eq}.$$

, the primal inequality constraint residual, which includes bounds and slacks:*r*_{ineq}$${r}_{ineq}=\overline{A}x-\overline{b}-s.$$

, the complementarity residual:*r*_{sz}=*r*_{sz}.*Sz*is the diagonal matrix of slack terms,*S*is the column matrix of Lagrange multipliers.*z*, the average complementarity:*r*_{c}$${r}_{c}=\frac{{s}^{T}z}{m}.$$

In a Newton step, the changes in * x*,

$$\left(\begin{array}{cccc}H& 0& -{A}_{eq}^{T}& -{\overline{A}}^{T}\\ {A}_{eq}& 0& 0& 0\\ \overline{A}& -I& 0& 0\\ 0& Z& 0& S\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta z\end{array}\right)=-\left(\begin{array}{c}{r}_{d}\\ {r}_{eq}\\ {r}_{ineq}\\ {r}_{sz}\end{array}\right).$$ | (9-1) |

However, a full Newton step might be infeasible, because of
the positivity constraints on * s* and

`quadprog`

shortens the step, if necessary,
to maintain positivity. Additionally, to maintain a "centered" position
in the interior, instead of trying to solve * s_{i}z_{i} = 0*, the algorithm
takes a positive parameter

* s_{i}z_{i}* =

`quadprog`

replaces * r_{sz}* in
the Newton step equation with

`quadprog`

reorders the
Newton equations to obtain a symmetric, more numerically stable system
for the predictor step calculation.After calculating the corrected Newton step, the algorithm performs more calculations to get both a longer current step, and to prepare for better subsequent steps. These multiple correction calculations can improve both performance and robustness. For details, see Gondzio [4].

**Full Predictor-Corrector. **The full predictor-corrector algorithm does not combine bounds
into linear constraints, so it has another set of slack variables
corresponding to the bounds. The algorithm shifts lower bounds to
zero. And, if there is only one bound on a variable, the algorithm
turns it into a lower bound of zero, by negating the inequality of
an upper bound.

$$\overline{A}$$ is the extended linear matrix
that includes both linear inequalities and linear equalities. $$\overline{b}$$ is the corresponding linear
equality vector. $$\overline{A}$$ also includes
terms for extending the vector * x* with slack variables

$$\overline{A}x=\left(\begin{array}{cc}{A}_{eq}& 0\\ A& I\end{array}\right)\left(\begin{array}{c}{x}_{0}\\ s\end{array}\right),$$

where *x*_{0} means the
original * x* vector.

The KKT conditions are

$$\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w=0\\ \overline{A}x=\overline{b}\\ x+t=u\\ {v}_{i}{x}_{i}=0,\text{}i=1,2,\mathrm{...},m\\ {w}_{i}{t}_{i}=0,\text{}i=1,2,\mathrm{...},n\\ x,v,w,t\ge 0.\end{array}$$ | (9-2) |

To find the solution * x*, slack variables and
dual variables to Equation 9-2,
the algorithm basically considers a Newton-Raphson step:

$$\left(\begin{array}{ccccc}H& -{\overline{A}}^{T}& 0& -I& I\\ \overline{A}& 0& 0& 0& 0\\ -I& 0& -I& 0& 0\\ V& 0& 0& X& 0\\ 0& 0& W& 0& T\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta t\\ \begin{array}{l}\Delta v\\ \Delta w\end{array}\end{array}\right)=-\left(\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w\\ \overline{A}x-\overline{b}\\ u-x-t\\ \begin{array}{l}VX\\ WT\end{array}\end{array}\right)=-\left(\begin{array}{c}{r}_{d}\\ {r}_{p}\\ {r}_{ub}\\ \begin{array}{l}{r}_{vx}\\ {r}_{wt}\end{array}\end{array}\right),$$ | (9-3) |

where * X*,

, the dual residual*r*_{d}, the primal residual*r*_{p}, the upper bound residual*r*_{ub}, the lower bound complementarity residual*r*_{vx}, the upper bound complementarity residual*r*_{wt}

The algorithm solves Equation 9-3 by first converting it to the symmetric matrix form

$$\left(\begin{array}{cc}-D& {\overline{A}}^{T}\\ \overline{A}& 0\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=-\left(\begin{array}{c}R\\ {r}_{p}\end{array}\right),$$ | (9-4) |

where

$$\begin{array}{l}D=H+{X}^{-1}V+{T}^{-1}W\\ R=-{r}_{d}-{X}^{-1}{r}_{vx}+{T}^{-1}{r}_{wt}+{T}^{-1}W{r}_{ub}.\end{array}$$

All the matrix inverses in the definitions of * D* and

To derive Equation 9-4 from Equation 9-3, notice that the
second row of Equation 9-4 is
the same as the second matrix row of Equation 9-3. The first row of Equation 9-4 comes from solving
the last two rows of Equation 9-3 for Δ* v* and
Δ

To solve Equation 9-4, the algorithm follows the essential elements of Altman and Gondzio [1]. The algorithm solves the symmetric system by an LDL decomposition. As pointed out by authors such as Vanderbei and Carpenter [2], this decomposition is numerically stable without any pivoting, so can be fast.

After calculating the corrected Newton step, the algorithm performs more calculations to get both a longer current step, and to prepare for better subsequent steps. These multiple correction calculations can improve both performance and robustness. For details, see Gondzio [4].

The full `quadprog`

predictor-corrector algorithm
is largely the same as that in the `linprog`

`'interior-point'`

algorithm,
but includes quadratic terms as well. See Predictor-Corrector.

[1] Altman, Anna and J. Gondzio. *Regularized
symmetric indefinite systems in interior point methods for linear
and quadratic optimization*. Optimization Methods and Software,
1999. Available
for download here.

[2] Vanderbei, R. J. and T. J. Carpenter. *Symmetric
indefinite systems for interior point methods*. Mathematical
Programming 58, 1993. pp. 1–32. Available
for download here.

The predictor-corrector algorithm iterates until it reaches a point that is feasible (satisfies the constraints to within tolerances) and where the relative step sizes are small. Specifically, define

$$\rho =\mathrm{max}\left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert {A}_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert {b}_{eq}\Vert \right)$$

The algorithm stops when all of these conditions are satisfied:

$$\begin{array}{c}{\Vert {r}_{p}\Vert}_{1}+{\Vert {r}_{ub}\Vert}_{1}\le \rho \text{TolCon}\\ {\Vert {r}_{d}\Vert}_{\infty}\le \rho \text{TolFun}\\ {r}_{c}\le \text{TolFun,}\end{array}$$

where

$${r}_{c}=\underset{i}{\mathrm{max}}\left(\mathrm{min}\left(\left|{x}_{i}{v}_{i}\right|,\left|{x}_{i}\right|,\left|{v}_{i}\right|\right),\mathrm{min}\left(\left|{t}_{i}{w}_{i}\right|,\left|{t}_{i}\right|,\left|{w}_{i}\right|\right)\right).$$

* r_{c}* essentially measures
the size of the complementarity residuals

`quadprog`

calculates a *merit function* * φ* at
every iteration. The merit function is a measure of feasibility.

`quadprog`

stops
if the merit function grows too large. In this case, `quadprog`

declares
the problem to be infeasible.The merit function is related to the KKT conditions for the problem—see Predictor-Corrector. Use the following definitions:

$$\begin{array}{c}\rho =\mathrm{max}\left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert {A}_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert {b}_{eq}\Vert \right)\\ {r}_{\text{eq}}={A}_{\text{eq}}x-{b}_{\text{eq}}\\ {r}_{\text{ineq}}=\overline{A}x-\overline{b}+s\\ {r}_{\text{d}}=Hx+c+{A}_{\text{eq}}^{T}{\lambda}_{\text{eq}}+{\overline{A}}^{T}{\overline{\lambda}}_{\text{ineq}}\\ g={x}^{T}Hx+{f}^{T}x-{\overline{b}}^{T}{\overline{\lambda}}_{\text{ineq}}-{b}_{\text{eq}}^{T}{\lambda}_{\text{eq}}.\end{array}$$

The notation $$\overline{A}$$ and $$\overline{b}$$ means the linear inequality
coefficients, augmented with terms to represent bounds for the sparse
algorithm. The notation $${\overline{\lambda}}_{\text{ineq}}$$ similarly represents
Lagrange multipliers for the linear inequality constraints, including
bound constraints. This was called * z* in Predictor-Corrector, and $${\lambda}_{\text{eq}}$$ was called

The merit function * φ* is

$$\frac{1}{\rho}\left(\mathrm{max}\left({\Vert {r}_{\text{eq}}\Vert}_{\infty},{\Vert {r}_{\text{ineq}}\Vert}_{\infty},{\Vert {r}_{\text{d}}\Vert}_{\infty}\right)+g\right).$$

If this merit function becomes too large, `quadprog`

declares
the problem to be infeasible and halts with exit flag `-2`

.

`trust-region-reflective`

`quadprog`

AlgorithmMany of the methods used in Optimization Toolbox™ solvers
are based on *trust regions,* a simple yet powerful
concept in optimization.

To understand the trust-region approach to optimization, consider
the unconstrained minimization problem, minimize * f*(

$$\underset{s}{\mathrm{min}}\left\{q(s),\text{}s\in N\right\}.$$ | (9-5) |

The current point is updated to be * x* +

The key questions in defining a specific trust-region approach
to minimizing * f*(

In the standard trust-region method ([48]), the quadratic approximation * q* is
defined by the first two terms of the Taylor approximation to

$$\mathrm{min}\left\{\frac{1}{2}{s}^{T}Hs+{s}^{T}g\text{suchthat}\Vert Ds\Vert \le \Delta \right\},$$ | (9-6) |

where * g* is the gradient of

$$\frac{1}{\Delta}-\frac{1}{\Vert s\Vert}=0.$$

Such algorithms provide an accurate solution to Equation 9-6. However, they
require time proportional to several factorizations of * H*.
Therefore, for large-scale problems a different approach is needed.
Several approximation and heuristic strategies, based on Equation 9-6, have been proposed
in the literature ([42] and [50]). The approximation approach followed
in Optimization Toolbox solvers is to restrict the trust-region
subproblem to a two-dimensional subspace

The two-dimensional subspace * S* is
determined with the aid of a preconditioned
conjugate gradient process described below. The solver defines

$$H\cdot {s}_{2}=-g,$$ | (9-7) |

or a direction of negative curvature,

$${s}_{2}^{T}\cdot H\cdot {s}_{2}<0.$$ | (9-8) |

The philosophy behind this choice of * S* is
to force global convergence (via the steepest descent direction or
negative curvature direction) and achieve fast local convergence (via
the Newton step, when it exists).

A sketch of unconstrained minimization using trust-region ideas is now easy to give:

Formulate the two-dimensional trust-region subproblem.

Solve Equation 9-6 to determine the trial step

.*s*If

, then*f*(*x*+*s*) <*f*(*x*).*x*=*x*+*s*Adjust Δ.

These four steps are repeated until convergence. The trust-region
dimension Δ is adjusted according to standard rules. In particular,
it is decreased if the trial step is not accepted, i.e., * f(x + s)
≥ f(x)*.
See [46] and [49] for a discussion of this aspect.

Optimization Toolbox solvers treat a few important special
cases of * f* with specialized functions: nonlinear
least-squares, quadratic functions, and linear least-squares. However,
the underlying algorithmic ideas are the same as for the general case.
These special cases are discussed in later sections.

The subspace trust-region method is used to determine a search direction. However, instead of restricting the step to (possibly) one reflection step, as in the nonlinear minimization case, a piecewise reflective line search is conducted at each iteration. See [45] for details of the line search.

A popular way to solve large symmetric positive definite systems
of linear equations * Hp = –g* is the
method of Preconditioned Conjugate Gradients (PCG). This iterative
approach requires the ability to calculate matrix-vector products
of the form

In a minimization context, you can assume that the Hessian matrix * H* is
symmetric. However,

Linear constraints complicate the situation described for unconstrained minimization. However, the underlying ideas described previously can be carried through in a clean and efficient way. The trust-region methods in Optimization Toolbox solvers generate strictly feasible iterates.

The general linear equality constrained minimization problem can be written

$$\mathrm{min}\left\{f(x)\text{suchthat}Ax=b\right\},$$ | (9-9) |

where * A* is an

The method used to solve Equation 9-9 differs from the unconstrained approach
in two significant ways. First, an initial feasible point *x*_{0} is
computed, using a sparse least-squares step, so that * Ax_{0} = b*.
Second, Algorithm PCG is replaced with Reduced Preconditioned Conjugate
Gradients (RPCG), see [46], in order
to compute an approximate reduced Newton step (or a direction of negative
curvature in the null space of

$$\left[\begin{array}{cc}C& {\tilde{A}}^{T}\\ \tilde{A}& 0\end{array}\right]\left[\begin{array}{c}s\\ t\end{array}\right]=\left[\begin{array}{c}r\\ 0\end{array}\right],$$ | (9-10) |

where $$\tilde{A}$$ approximates * A* (small
nonzeros of

The box constrained problem is of the form

$$\mathrm{min}\left\{f(x)\text{suchthat}l\le x\le u\right\},$$ | (9-11) |

where * l* is a vector of lower bounds, and

The scaled modified Newton step arises from examining the Kuhn-Tucker necessary conditions for Equation 9-11,

$${\left(D(x)\right)}^{-2}g=0,$$ | (9-12) |

where

$$D(x)=\text{diag}\left({\left|{v}_{k}\right|}^{-1/2}\right),$$

and the vector * v*(

If

and*g*< 0_{i}then*u*< ∞_{i}*v*=_{i}*x*–_{i}*u*_{i}If

and*g*≥ 0_{i}then*l*> –∞_{i}*v*=_{i}*x*–_{i}*l*_{i}If

and*g*< 0_{i}then*u*= ∞_{i}*v*= –1_{i}If

and*g*≥ 0_{i}then*l*= –∞_{i}*v*= 1_{i}

The nonlinear system Equation 9-12 is not differentiable everywhere.
Nondifferentiability occurs when * v_{i} = 0*. You can avoid
such points by maintaining strict feasibility, i.e., restricting

The scaled modified Newton step * s_{k}* for
the nonlinear system of equations given by Equation 9-12 is defined as the solution to the
linear system

$$\widehat{M}D{s}^{N}=-\widehat{g}$$ | (9-13) |

at the * k*th iteration, where

$$\widehat{g}={D}^{-1}g=\text{diag}\left({\left|v\right|}^{1/2}\right)g,$$ | (9-14) |

and

$$\widehat{M}={D}^{-1}H{D}^{-1}+\text{diag}(g){J}^{v}.$$ | (9-15) |

Here * J^{v}* plays
the role of the Jacobian of |

Second, reflections are used to increase
the step size. A (single) reflection step is defined as follows. Given
a step * p* that intersects a bound constraint, consider
the first bound constraint crossed by

Was this topic helpful?