/******************* Quadratic Programming ************************
Author: Rob Dittmar dittmar@stls.frb.org
Research Department
Federal Reserve Bank of St. Louis
Date: 20 Sep 1996
This code is provided free, without performance guarantees,
for public, non-commercial use.
To solve a quadratic programming problem, you usually write the first order
conditions in the form of a what is referred to as a linear complementarity
problem and use a variant of the simplex algorithm to solve it.
If the quadratic programming problem to be solved is of the form
Min c'x + x'Qx
s. t. A x >= b
x >= 0
then the first order conditions take the form
v = 2 Q x - A' mu + c
s = A x - b
x, mu, v, s >= 0
x'v + s'mu = 0
where v are the multipliers on the non-negativity constraints, the mus are
multipliers on the inequality constraints, and the s are slack variables for
the inequality contraints.
This problem is in the linear-complementarity problem form which is to find
vectors w and z to solve
w = M z + q
w, z, >= 0
w'z = 0.
I have attached a short procedure that will solve the linear-complementarity
problem. Write the FOC's for your quadratic problem in the above form and
solve the associated linear complementarity problem.
****************************************************************************/
/* Smallest allowable pivot element */
_pivtol = 1e-8;
/* Procedure to perform partial pivoting on a matrix. The
procedure takes a matrix argument, tableau, and two scalar
arguments, r and c. The procedure performs a row reduction
using the [r,c]-element of the matrix as a pivot. The returned
matrix will have a 1 as its [r,c]-element, and zeros everywhere
else in column c. */
PROC (1) = rowred(tableau,r,c);
LOCAL reduced,pivot,i;
reduced = tableau;
pivot = tableau[r,.]/tableau[r,c];
i = 1;
DO WHILE i <= ROWS(tableau);
reduced[i,.] = tableau[i,.] - tableau[i,c]*pivot;
i = i + 1;
ENDO;
reduced[r,.] = pivot;
RETP(reduced);
ENDP;
/* Procedure to solve the linear complementarity problem:
w = M z + q
w and z >= 0
w'z = 0
The procedure takes the matrix M and vector q as arguments. The
procedure has three returns. The first and second returns are
the final values of the vectors w and z found by complementary
pivoting. The third return is a 2 by 1 vector. Its first
component is a 1 if the algorithm was successful, and a 2 if a
ray termination resulted. The second component is the value of
the artificial varible upon termination of the algorithm. */
PROC (3) = LCPSolve(M,q);
LOCAL dimen,tableau,basis,locat,cand,rayTerm,eMs,quots,
oldVar,vars,w,z,retcode;
/* Dimension of LCP */
dimen = ROWS(M);
/* Create initial tableau */
tableau = EYE(dimen)~(-M)~(-ONES(dimen,1))~q;
/* Let artificial variable enter the basis */
basis = SEQA(1,1,dimen);
locat = MININDC(tableau[.,2*dimen+2]);
basis[locat] = 2*dimen + 1;
cand = locat + dimen;
tableau = rowred(tableau,locat,2*dimen+1);
/* Perform complementary pivoting */
rayTerm = 0;
DO WHILE MAXC(basis) == 2*dimen+1 AND NOT(rayTerm);
eMs = tableau[.,cand].*MISS(tableau[.,cand] .> 0, 0);
quots = tableau[.,2*dimen+2]./eMs;
locat = MININDC(quots);
IF abs(eMs[locat]) > _pivtol AND NOT ISMISS(PACKR(quots));
@ Reduce tableau @
tableau = rowred(tableau,locat,cand);
oldVar = basis[locat];
@ New variable enters the basis @
basis[locat] = cand;
@ Select next candidate for entering the basis @
IF oldVar > dimen;
cand = oldVar - dimen;
ELSE;
cand = oldVar + dimen;
ENDIF;
ELSE;
rayTerm = 1;
ENDIF;
ENDO;
@ Return solution to LCP @
vars = ZEROS(2*dimen+1,1);
vars[basis] = tableau[.,2*dimen+2];
w = vars[1:dimen];
z = vars[(dimen+1):2*dimen];
retcode = vars[(2*dimen+1)];
IF rayTerm;
retcode = 2|retcode;
ELSE;
retcode = 1|retcode;
ENDIF;
RETP(w,z,retcode);
ENDP;