Applications of QR, RUPDT, and RDODT Subroutines
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: QR */
/* TITLE: Applications of QR, RUPDT, and RDODT Subroutines */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: MATRIX */
/* PROCS: IML */
/* DATA: */
/* */
/* SUPPORT: MZC UPDATE: */
/* REF: */
/* MISC: */
/* */
/****************************************************************/
proc iml;
a={ -74 80 18 -11 -4 -8,
14 -69 21 28 0 7,
66 -72 -5 7 1 4,
-12 66 -30 -23 3 -3,
3 8 -7 -4 1 0,
4 -12 4 4 0 1};
/* Second column of b is orthogonal to the columns of a */
/* third column of b is equal the sum of the first two */
b= { 51 -56 -5,
-61 52 -9,
-56 764 708,
69 4096 4165,
10 -13276 -13266,
-12 8421 8409};
aa = a[,1:3];
m= nrow(aa);
n= ncol(aa);
p= ncol(b);
call QR(q,r,piv,lindep,aa);
r0 = j(n,n,0.);
qtb = j(n,p,0.);
ssq = j(1,p,0.);
bb=b;
CALL RUPDT(rup,bup,sup,r0,aa,qtb,bb,ssq);
/* r and rup should be identical(maybe in different sign) */
print r rup;
/* check should be close to bup (maybe in different sign */
check=(q`*b)[1:3,];
print check bup;
/* get ssq from from the last 3 rows of q`*b */
check=(q`*b)[4:6,];
check=sqrt(diag(check`*check));
/* diagonal elements of check should be close to sup */
print check sup;
/* downdateing the last row of aa */
z=aa[6,]; bb=b[6,];
CALL Rdodt(def,rup1,bup1,sup1,rup,z,bup,bb,sup);
/* def should be 0 */
print def;
aa = a[1:5,1:3];
m= nrow(aa);
n= ncol(aa);
p= ncol(b);
/* updateing the first 5 rows of aa recursively */
r = j(n,n,0.);
qtb = j(n,p,0.);
ssq = j(1,p,0.);
do i = 1 to m;
z = aa[i,1:n]; bb = b[i,1:p];
CALL RUPDT(rup,bup,sup,r,z,qtb,bb,ssq);
r = rup; qtb = bup; ssq = sup;
end;
/* updateing the first 5 rows of aa at one time */
r = j(n,n,0.);
qtb = j(n,p,0.);
ssq = j(1,p,0.);
z = aa[1:5,1:n]; bb = b[1:5,1:p];
CALL RUPDT(rup2,bup2,sup2,r,z,qtb,bb,ssq);
/* each group of matrices should be identical */
print sup sup1 sup2, rup rup1 rup2, bup bup1 bup2;
proc iml;
a={ -74 80 18 -11 -4 -8,
14 -69 21 28 0 7,
66 -72 -5 7 1 4,
-12 66 -30 -23 3 -3,
3 8 -7 -4 1 0,
4 -12 4 4 0 1};
/* Second column of b is orthogonal to the columns of a */
/* third column of b is equal the sum of the first two */
b= { 51 -56 -5,
-61 52 -9,
-56 764 708,
69 4096 4165,
10 -13276 -13266,
-12 8421 8409};
aa = a[1:2,1:3]//a[5:6,1:3];
m= nrow(aa);
n= ncol(aa);
p= ncol(b);
call QR(q,rr,piv,lindep,aa);
check=q`*(b[1:2,]//b[5:6,]);
checkssq=sqrt(check[4,]##2);
r0 = j(n,n,0.);
qtb = j(n,p,0.);
ssq = j(1,p,0.);
bb=b[1:2,]//b[5:6,];
CALL RUPDT(rup,bup,sup,r0,aa,qtb,bb,ssq);
/* downdateing the middle 2 rows of a[,1:3] */
CALL RUPDT(rupp,bup,supp,r0,a[,1:3],qtb,b,ssq);
z=a[3:4,1:3]; bb=b[3:4,];
CALL Rdodt(def,rup1,bup1,sup1,rupp,z,bup,bb,supp);
/* rr, rup and rup1 should be identical(maybe in different */
/* sign), */
/* checkssq, sup and sup1 should be same, def should be 0 */
print rr rup rup1, checkssq sup sup1, def;