Quadratic and Geometric Programming

A quadratic program (QP) is the simplest nonlinear program that is still genuinely nonlinear: the objective is quadratic, every constraint is linear. It sits exactly one step beyond the linear programs of Part 2, and it is the local model that Chapter 21's sequential methods solve at every iteration — which is why it deserves a treatment of its own.

Here \(\mathbf{Q}\) is symmetric and \(n \times n\) — any non-symmetric matrix may be replaced by \(\tfrac12(\mathbf{Q} + \mathbf{Q}^{\top})\) without changing the value of the quadratic form, so symmetry costs nothing. The feasible region is a polyhedron, exactly as in linear programming; only the objective has curvature. Setting \(\mathbf{Q} = \mathbf{0}\) recovers an LP.

Everything hinges on the definiteness of \(\mathbf{Q}\), read exactly as in Section 19-5. Since \(\nabla^2 f = \mathbf{Q}\) everywhere, \(f\) is convex iff \(\mathbf{Q}\) is positive semidefinite. A convex QP over a polyhedron is the best-behaved nonlinear program there is: any KKT point is a global minimiser, and it can be found in finitely many pivots. An indefinite \(\mathbf{Q}\) makes the problem NP-hard in general — local minima proliferate, and the methods below only certify a KKT point.

Unlike an LP, a QP's optimum need not sit at a vertex. If the unconstrained minimiser \(-\mathbf{Q}^{-1}\mathbf{c}\) happens to be feasible, it is the answer and no constraint matters. Otherwise the solution is pushed onto the boundary — onto a face, an edge, or a vertex — wherever the elliptical contours first touch the polyhedron.

Apply the KKT conditions of Chapter 21 to the standard QP. Attach multipliers \(\boldsymbol{\lambda} \ge \mathbf{0}\) to \(\mathbf{A}\mathbf{x} \le \mathbf{b}\) and \(\boldsymbol{\mu} \ge \mathbf{0}\) to \(\mathbf{x} \ge \mathbf{0}\), and add slacks \(\mathbf{s} \ge \mathbf{0}\) to the inequalities. The stationarity condition \(\nabla f + \mathbf{A}^{\top}\boldsymbol{\lambda} - \boldsymbol{\mu} = \mathbf{0}\) is where the magic happens: because \(f\) is quadratic, \(\nabla f = \mathbf{Q}\mathbf{x} + \mathbf{c}\) is linear in \(\mathbf{x}\).

Read what has happened. Every condition is a linear equation or a sign restriction, except the two complementarity products. Strip those away and you are left with a system of \(n + m\) linear equations in \(2(n+m)\) non-negative unknowns — precisely the shape of a linear-programming constraint set. All the nonlinearity of the original problem has been squeezed into the combinatorial statement

The pairs are the complementary pairs; a solution in which at most one member of each pair is basic is a complementary basic feasible solution. Since the number of ways to choose one from each pair is finite, so is the search — and for convex QPs, any such solution found is globally optimal. Section 22-4 turns this observation into an algorithm.

Drop the inequalities and the combinatorics vanish. With only equalities, \(\min \tfrac12\mathbf{x}^{\top}\mathbf{Q}\mathbf{x} + \mathbf{c}^{\top}\mathbf{x}\) subject to \(\mathbf{A}\mathbf{x} = \mathbf{b}\), the Lagrangian is \(L = \tfrac12\mathbf{x}^{\top}\mathbf{Q}\mathbf{x} + \mathbf{c}^{\top}\mathbf{x} + \boldsymbol{\lambda}^{\top}(\mathbf{A}\mathbf{x} - \mathbf{b})\), and setting both blocks of partials to zero gives one symmetric linear system.

There is no iteration here at all: one factorisation and the equality-constrained QP is solved exactly, complete with its multipliers. The system is nonsingular whenever \(\mathbf{A}\) has full row rank and \(\mathbf{Q}\) is positive definite on the null space of \(\mathbf{A}\) — that is, the objective curves upward in every direction that stays feasible. Note the matrix is symmetric but indefinite (the zero block guarantees it), so a plain Cholesky factorisation will not do.

This is the engine inside almost every practical QP solver. An active-set method guesses which inequalities are tight, treats that guess (the working set) as equalities, solves the KKT system above, then repairs the guess — dropping a constraint whose multiplier came out negative (it was pushing the wrong way) or adding one that the step would have violated. Each repair is a single pivot, and the process is finite.

Wolfe's method is the classic way to exploit the LCP structure of Section 22-2 with nothing but a simplex code. The KKT equations are linear, so introduce artificial variables \(\mathbf{a} \ge \mathbf{0}\) and run a phase-one LP that drives them to zero. If the artificials reach zero, the remaining variables satisfy the KKT equations — and if we also never let a complementary pair go basic together, they satisfy complementarity too.

The procedure, then, is short:

Two cousins are worth naming. Beale's method partitions variables into basic and non-basic and expresses the objective in the non-basics, choosing the entering variable by the sign of the resulting free derivative — it handles the case where the optimum is interior to a face more naturally. Lemke's complementary pivoting algorithm attacks the LCP directly with an artificial variable \(z_0\) and a clever "almost-complementary path", and it remains the standard textbook LCP solver. Modern codes have largely replaced all three with active-set and interior-point methods, but the logic — find the active set, then solve a linear system — is unchanged.

Geometric programming (GP), developed by Duffin, Peterson and Zener in the 1960s, addresses a completely different structure: objectives built from sums of positive terms with arbitrary real exponents. Such functions describe engineering design almost perfectly — costs of material scale with area, of pumping with velocity cubed, of inventory inversely with order size.

Each \(u_k\) is a monomial (in the GP sense): a positive coefficient times a product of powers, where the exponents \(a_{ik}\) may be any real numbers — negative, fractional, whatever the physics demands. A sum of monomials is a posynomial ("positive polynomial"). Note the two standing restrictions: coefficients must be positive and variables must be strictly positive. Both are natural for a design problem in which \(x_i\) is a length, a flow rate, or a thickness.

with every \(g_m\) a posynomial. The form looks restrictive but is not: any constraint \(h(\mathbf{x}) \le d\) with posynomial \(h\) and \(d > 0\) becomes \(h(\mathbf{x})/d \le 1\), and equality constraints are allowed when both sides are monomials. A posynomial is not convex in \(\mathbf{x}\) — but under the substitution \(x_i = e^{y_i}\) each term becomes \(e^{\mathbf{a}_k^{\top}\mathbf{y} + \ln c_k}\), and a sum of exponentials of affine functions is convex. This is why GP is a convex problem in disguise, and why its local optimum is always global.

The whole of geometric programming rests on one classical inequality: a weighted arithmetic mean is never smaller than the corresponding weighted geometric mean.

Now play the trick. Given a posynomial \(g(\mathbf{x}) = \sum_k u_k(\mathbf{x})\) and any weights \(\delta_k\) summing to one, substitute \(U_k = u_k(\mathbf{x})/\delta_k\). The left side becomes \(\sum_k u_k = g(\mathbf{x})\) exactly, and the inequality reads

Look at the last product: the variables have collected into a single power of \(x_i\) with exponent \(\sum_k a_{ik}\delta_k\). If we insist that every one of those exponents be zero, the bound stops depending on \(\mathbf{x}\) altogether — it becomes a number, valid for every feasible design. Those two demands on the weights are the dual constraints, and maximising the resulting number is the dual problem.

Stare at what has been achieved. The primal is a nonlinear minimisation over \(n\) design variables; the dual constraints are \(n + 1\) linear equations in \(N\) weights. The nastiness has migrated from the objective into a linear system — and the dual objective involves the design variables not at all.

The dual has \(N\) unknowns tied down by \(n + 1\) linear equations. The gap between them measures how much work is left after the linear algebra is done, and it earns its own name.

The zero-degree-of-difficulty case is the showpiece of the subject: the optimal cost is obtained from a linear system, before a single design variable is known. And when the constraints are consistent but \(D < 0\), the dual system is over-determined and the formulation should be re-examined.

Recovering the design is the last step, and it follows from the condition under which AM–GM holds with equality — namely when all the \(U_k\) are equal. Tracing that back gives the interpretation that makes GP so useful to engineers:

Each \(\delta_k^\star\) is the fraction of the total optimal cost contributed by term \(k\). That is a design insight in its own right — it says where the money goes in an optimal design, and it is often invariant to the coefficients. Setting each monomial equal to its share \(\delta_k^\star f^\star\) yields \(N\) equations in the \(n\) unknowns \(x_i\); taking logarithms makes them linear, and any \(n\) independent ones give \(\mathbf{x}^\star\).

Constrained GPs extend the same idea. Each constraint \(g_m(\mathbf{x}) \le 1\) contributes its own block of weights \(\delta_{mk}\), the normality condition applies only to the objective's block, orthogonality runs over all blocks, and each constraint acquires a multiplier \(\lambda_m = \sum_k \delta_{mk}\) appearing in the dual objective as \(\prod_m \prod_k (c_{mk}\lambda_m/\delta_{mk})^{\delta_{mk}}\). A constraint with \(\lambda_m^\star = 0\) is inactive at the optimum — the same complementary-slackness message as Section 22-2, arriving by a different road.

For the QPs, form \(\mathbf{Q}\), test definiteness, and write the KKT conditions before solving. For the GPs, always count \(D\) first. Difficulty rises down the list.