Quadratic Programming Algorithm
/******************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: QUADPRG */
/* TITLE: Quadratic Programming Algorithm */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: MATH MATRIX */
/* PROCS: IML */
/* DATA: */
/* */
/* SUPPORT: RHD UPDATE: */
/* REF: METHOD OF THEIL & VAN DE PANNE */
/* MISC: CONVERTED FROM MATRIX TO IML USING MATIML */
/******************************************************************/
*---------------- QUADRATIC PROGRAMMING IN IML ------------------*
| THE FOLLOWING MATRIX PROGRAM IMPLEMENTS THE METHOD OF THEIL & |
| VAN DE PANNE (MANAGEMENT SCIENCE V.7,PP1-20 (1960)) FOR |
| MAXIMIZING A QUADRATIC FUNCTION SUBJECT TO LINEAR CONSTRAINTS. |
| THE METHOD STARTS WITH THE OPTIMAL UNCONSTRINED VALUES AND THEN|
| SEARCHES VARIOUS CONFIGURATIONS OF THE CONSTRAINTS FOR THE |
| OPTIMAL FEASIBLE VALUES. |
| |
| THE OBVIOUS APPLICATION OF THIS METHOD IS IN LEAST SQUARES |
| REGRESSION WITH INEQUALITY CONSTRAINTS ON THE PARAMETERS. |
*---------------------------------------------------JPS 23OCT79--*;
PROC IML;
RESET AUTONAME LINESIZE=132 PAGESIZE=60;
START MAIN;
*---DATA FOR THIS PROBLEM---;
XPY={18,16,22,20};
XPX={ 6 1 8 0,
1 10 1 4,
8 1 17 3,
0 4 3 11};
R= {-1 0 0 0,
0 -1 0 0,
0 0 -1 0,
0 0 0 -1,
1 1 1 1,
5 0 10 0,
0 4 0 5};
C={0 0 0 0 1.66666667 2 3}`;
*---FIND DIMENSIONS OF PROBLEM---;
NX= NROW(XPX); * NUMBER OF VARIABLES;
NR= NROW(R); * NUMBER OF INEQ. CONSTRAINTS;
NT=NX+NR+1; * TOTAL COLUMNS IN Z;
IR=NX+1:NX+NR; * INDEX FOR CONSTRAINTS;
Z1= J(NR,NR,0);
*---SET UP CROSSPRODUCTS---;
Z= ( XPX || R` || XPY )//
( R || Z1 || C )//
(XPY` || C` || 0 );
*---SOLVE FOR UNCONSTRAINED MAXIMUM---;
Z= SWEEP(Z,1:NX);
ZTRY=Z;
I=1;
NS=1;
*---LOOK FOR UNSATISFIED CONSTRAINTS---;
STABLE= Z[IR,NT]<0;
CTABLE= J(NR,1,0);
IF ( SUM(STABLE)=0) THEN GOTO DONE;
LOOP: C= CTABLE[,I];
S= STABLE[,I];
NLOC= LOC(S);
NNEG= NCOL(NLOC);
DO J=1 TO NNEG;
IJ= NLOC[,J];
CTRY=C;
CTRY[IJ,]=1;
* FORM NEW COMBINATION;
DO K=I+1 TO NS;
* SEE IF ALREADY IN STABLE;
IF ( CTRY= CTABLE[,K]) THEN GOTO DUP;
END;
CTABLE=CTABLE||CTRY;
NS=NS+1;
LINK SWP;
STABLE=STABLE||STRY;
IF ( SUM(STRY)=0) THEN LINK TEST;
DUP: END;
I=I+1;
IF ( I<=NS) THEN GOTO LOOP;
PRINT " -----NO FEASIBLE SOLUTION-----";
STOP;
*---ROUTINE TO SWEEP IN CONSTRAINTS---;
SWP: ITRY= LOC(CTRY)+NX;
ZTRY= SWEEP(Z,ITRY);
STRY= ZTRY[IR,NT]<0;
STRY=STRY#^CTRY;
RETURN;
*---ROUTINE TO TEST IF OPTIMAL---;
TEST: CTLOC= LOC(CTRY);
NTLOC= NCOL(CTLOC);
DO III=1 TO NTLOC;
IIJ= CTLOC[,III];
CT=CTRY;
CT[IIJ,]=0;
LOCCT= LOC(CT);
IF ( NCOL(LOCCT)>0) THEN ZCT= SWEEP(Z,LOCCT+NX);
ELSE ZCT=Z;
CJ=( ZCT[IR,NT] <0)#^CT;
IF ( ^CJ[IIJ,]) THEN RETURN;
END;
*---PRINT OUT RESULTS---;
DONE: PRINT CTRY[L="CONSTRAINT CONFIGURATION FOR FEASIBLE OPTIMUM"];
BETA= ZTRY[1:NX,NT];
PRINT BETA[L="PARAMETER ESTIMATES FOR OPTIMAL"];
LAMBDA= ZTRY[IR,NT];
PRINT LAMBDA[L="CONSTRAINT SLACKS AND LAGRANGIANS"];
PRINT CTABLE[F=4. L="CONFIGURATIONS TRIED"];
PRINT STABLE[F=4. L="PATTERN OF CONSTRAINT VIOLATIONS FOR CONFIGURATIONS"];
STOP;
FINISH MAIN;
RUN MAIN;