We want to find the solution of:

where and .

We will assume in this chapter that is positive definite. There is thus only one solution. We will first see how to handle the following simpler problem:

Equality constraints

where can be any vector. Indeed, we have the following: . The matrix has linearly independent columns which are inside the null-space of and therefore act as bases vectors (or reduced coordinate directions) for the null space. At any point any

(9.4) |

where are the components (or reduced variables) in each reduced coordinate direction (see Figure 9.1).

Combining equation 9.2 and 9.3, we obtain the reduced quadratic function:

(9.5) |

If is positive definite then a minimizer exists which solves the linear system:

(9.6) |

The solution is computed using Cholesky factorization. Once is known we can compute using equation 9.3 and using the secant equation 13.29: . Let's recall equation 13.19:

Let's now compute : From equation 9.7 we have:

Pre-Multiply 9.8 by and using , we have:

(9.9) |

Depending on the choice of and , a number of methods exist. We will obtain and from a QR factorization of the matrix (see annexe for more information about the QR factorization of a matrix). This can be written:

(9.10) |

where is an orthogonal matrix, is an upper triangular matrix, and . The choices

(9.11) |

have the correct properties. Moreover the vector in equation 9.3 and figure 9.2 is orthogonal to the constraints. The reduced coordinate directions are also mutually orthogonal. is calculated by forward substitution in followed by forming . The multipliers are calculated by backward substitution in .

Let's define, as usual, . If is infeasible then the length of the step is chosen to solve:

If then a new constraint becomes active, defined by the index which achieves the min in 9.13, and this index is added to the active set . Suppose we have the solution of the EP. We will now make a test to determine if a constraint should become inactive. All the constraints which have a negative associated Lagrange multiplier can become inactive. To summarize the algorithm, we have:

- Set
*a feasible point*. - Solve the EP.
- Compute the Lagrange multipliers . If for all constraints then terminate. Remove the constraints which have negative .
- Solve 9.13 and activate a new constraint if necessary.
- set and go to (2).

subject to | (9.14) |

is given by

subject to | (9.15) |

By eliminating from the first constraint equation, we obtain the bound constrained problem:

This problem can be solved by means of the gradient projection method, which normally allows us to identify the active set more rapidly than with classical active set methods. The matrix is semi-definite positive when is positive definite. This is good. Unfortunately, if we have linearly dependent constraints then the matrix becomes singular (this is always the case when ). This lead to difficulties when solving 9.16. There is, for example, no Cholesky factorization possible.

My first algorithm used gradient projection to identify the active set. It then project the matrix into the space of the active constraints (the projection is straight forward and very easy) and attempt a Cholesky factorization of the reduced matrix. This fails very often. When it fails, it uses "steepest descent algorithm" which is sloooooooow and useless. The final implemented QP works in the primal space and use QR factorization to do the projection.

(9.17) |

The solution to the EP (equation 9.12) is computed in one step using a Cholesky factorization of . This is very fast but, for badly scaled problem, this can lead to big rounding errors in the solution. The technique to choose which constraint enters the active set is very primitive (it's based on equation 9.13) and can also lead to big rounding errors. The algorithm which finds an initial feasible point when the given starting point is infeasible could be improved. All the linear algebra operations are performed with dense matrix.

This QP algorithm is very far from the "state-of-the-art". Some numerical results shows that the QP algorithm is really the weak point of all the optimization code. Nevertheless, for most problems, this implementation gives sufficient results (see numerical results).