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;