Chapter Contents |
Previous |
Next |
Language Reference |
proc iml; a = {1 1 0 0 1 0 0, 1 1 0 0 1 0 0, 1 1 0 0 0 1 0, 1 1 0 0 0 0 1, 1 0 1 0 1 0 0, 1 0 1 0 0 1 0, 1 0 1 0 0 1 0, 1 0 1 0 0 0 1, 1 0 0 1 1 0 0, 1 0 0 1 0 1 0, 1 0 0 1 0 0 1, 1 0 0 1 0 0 1}; a = a || uniform(j(12,1,1)); aa = a` * a; m = nrow(a); n = ncol(a);Applying the ROOT function to the coefficient matrix A' A of the normal equations,
r1 = root(aa); ss1 = ssq(aa - r1` * r1); print ss1 r1 [format=best6.];generates an upper triangular matrix R1 where linearly dependent rows are zeroed out, and you can verify that A' A = R'1 R1.
ord = j(n,1,0); call qr(q,r2,pivqr,lindqr,a,ord); ss2 = ssq(aa[pivqr,pivqr] - r2` * r2); print ss2 r2 [format=best6.];
r3 = shape(0,n,n); call rupdt(rup,bup,sup,r3,a); r3 = rup; ss3 = ssq(aa - r3` * r3); print ss3 r3 [format=best6.];
call rzlind(lind,r4,bup,r3); ss4 = ssq(aa - r4` * r4); print ss4 r4 [format=best6.];
b = uniform(j(12,1,1)); ab = a` * b; print b a [format=best6.];
call svd(u,d,v,a); t = 1e-12 * d[1]; do i=1 to n; if d[i] < t then d[i] = 0.; else d[i] = 1. / d[i]; end; x1 = v * diag(d) * u` * b; len1 = x1` * x1; ss1 = ssq(a * x1 - b); x1 = x1`; print ss1 len1, x1 [format=best6.];The solution obtained by singular value decomposition, , is of minimum Euclidean length.
ord = j(n,1,0); call qr(qtb,r2,pivqr,lindqr,a,ord,b); nr = n - lindqr; r = r2[1:nr,1:nr]; x2 = shape(0,n,1); x2[pivqr] = trisolv(1,r,qtb[1:nr]) // j(lindqr,1,0.); len2 = x2` * x2; ss2 = ssq(a * x2 - b); x2 = x2`; print ss2 len2, x2 [format=best6.];Note that the residual sum of squares is minimal, but the solution is not of minimum Euclidean length.
r1 = root(aa); nr = n - lind; piv = shape(0,n,1); j1 = 1; j2 = nr + 1; do i=1 to n; if r1[i,i] ^= 0 then do; piv[j1] = i; j1 = j1 + 1; end; else do; piv[j2] = i; j2 = j2 + 1; end; end;Now compute by solving the equation .
r = r1[piv[1:nr],piv[1:nr]]; x = trisolv(2,r,ab[piv[1:nr]]); x = trisolv(1,r,x); x3 = shape(0,n,1); x3[piv] = x // j(lind,1,0.); len3 = x3` * x3; ss3 = ssq(a * x3 - b); x3 = x3`; print ss3 len3, x3 [format=best6.];Note that the residual sum of squares is minimal, but the solution is not of minimum Euclidean length.
r3 = shape(0,n,n); qtb = shape(0,n,1); call rupdt(rup,bup,sup,r3,a,qtb,b); r3 = rup; qtb = bup; call rzlind(lind,r4,bup,r3,,qtb); qtb = bup[piv[1:nr]]; x = trisolv(1,r4[piv[1:nr],piv[1:nr]],qtb); x4 = shape(0,n,1); x4[piv] = x // j(lind,1,0.); len4 = x4` * x4; ss4 = ssq(a * x4 - b); x4 = x4`; print ss4 len4, x4 [format=best6.];Since the matrices R4 and R1 are the same (except for the signs of rows), the solution is the same as .
r = r4[piv[1:nr],]`; call qr(q,r5,piv2,lin2,r); y = trisolv(2,r5,qtb); x5 = q * (y // j(lind,1,0.)); len5 = x5` * x5; ss5 = ssq(a * x5 - b); x5 = x5`; print ss5 len5, x5 [format=best6.];
ord = j(n,1,0); call qr(qtb,r2,pivqr,lindqr,a,ord,b); nr = n - lindqr; r = r2[1:nr,]`; call qr(q,r5,piv2,lin2,r); y = trisolv(2,r5,qtb[1:nr]); x6 = shape(0,n,1); x6[pivqr] = q * (y // j(lindqr,1,0.)); len6 = x6` * x6; ss6 = ssq(a * x6 - b); x6 = x6`; print ss6 len6, x6 [format=best6.];The solution obtained by complete QR decomposition has minimum Euclidean length.
ord = j(n,1,0); call qr(qtb,r2,pivqr,lindqr,a,ord,b); nr = n - lindqr; r = r2[1:nr,1:nr]`; z = r2[1:nr,nr+1:n]`; call lupdt(lup,bup,sup,r,z); rd = trisolv(3,lup,r2[1:nr,]); rd = trisolv(4,lup,rd); x7 = shape(0,n,1); x7[pivqr] = rd` * qtb[1:nr,]; len7 = x7` * x7; ss7 = ssq(a * x7 - b); x7 = x7`; print ss7 len7, x7 [format=best6.];The solution obtained by complete QR decomposition has minimum Euclidean length.
r3 = shape(0,n,n); qtb = shape(0,n,1); call rupdt(rup,bup,sup,r3,a,qtb,b); r3 = rup; qtb = bup; call rzlind(lind,r4,bup,r3,,qtb); nr = n - lind; qtb = bup; r = r4[piv[1:nr],piv[1:nr]]`; z = r4[piv[1:nr],piv[nr+1:n]]`; call lupdt(lup,bup,sup,r,z); rd = trisolv(3,lup,r4[piv[1:nr],]); rd = trisolv(4,lup,rd); x8 = shape(0,n,1); x8 = rd` * qtb[piv[1:nr],]; len8 = x8` * x8; ss8 = ssq(a * x8 - b); x8 = x8`; print ss8 len8, x8 [format=best6.];The solution obtained by complete QR decomposition has minimum Euclidean length. The same result can be obtained with the APPCORT or COMPORT call.
ga = ginv(a); t1 = a * ga; t2 = t1`; t3 = ga * a; t4 = t3`; ss1 = ssq(t1 - t2) + ssq(t3 - t4) + ssq(t1 * a - a) + ssq(t3 * ga - ga); print ss1, ga [format=best6.];
call svd(u,d,v,a); do i=1 to n; if d[i] <= 1e-10 * d[1] then d[i] = 0.; else d[i] = 1. / d[i]; end; ga = v * diag(d) * u`; t1 = a * ga; t2 = t1`; t3 = ga * a; t4 = t3`; ss2 = ssq(t1 - t2) + ssq(t3 - t4) + ssq(t1 * a - a) + ssq(t3 * ga - ga); print ss2;
ord = j(n,1,0); call qr(q1,r2,pivqr,lindqr,a,ord); nr = n - lindqr; q1 = q1[,1:nr]; r = r2[1:nr,]`; call qr(q2,r5,piv2,lin2,r); tt = trisolv(4,r5`,q1`); ga = shape(0,n,m); ga[pivqr,] = q2 * (tt // shape(0,n-nr,m)); t1 = a * ga; t2 = t1`; t3 = ga * a; t4 = t3`; ss3 = ssq(t1 - t2) + ssq(t3 - t4) + ssq(t1 * a - a) + ssq(t3 * ga - ga); print ss3;
ord = j(n,1,0); call qr(q,r2,pivqr,lindqr,a,ord); nr = n - lindqr; r = r2[1:nr,1:nr]`; z = r2[1:nr,nr+1:n]`; call lupdt(lup,bup,sup,r,z); rd = trisolv(3,lup,r2[1:nr,]); rd = trisolv(4,lup,rd); ga = shape(0,n,m); ga[pivqr,] = rd` * q[,1:nr]`; t1 = a * ga; t2 = t1`; t3 = ga * a; t4 = t3`; ss4 = ssq(t1 - t2) + ssq(t3 - t4) + ssq(t1 * a - a) + ssq(t3 * ga - ga); print ss4;
r3 = shape(0,n,n); y = i(m); qtb = shape(0,n,m); call rupdt(rup,bup,sup,r3,a,qtb,y); r3 = rup; qtb = bup; call rzlind(lind,r4,bup,r3,,qtb); nr = n - lind; qtb = bup; r = r4[piv[1:nr],piv[1:nr]]`; z = r4[piv[1:nr],piv[nr+1:n]]`; call lupdt(lup,bup,sup,r,z); rd = trisolv(3,lup,r4[piv[1:nr],]); rd = trisolv(4,lup,rd); ga = shape(0,n,m); ga = rd` * qtb[piv[1:nr],]; t1 = a * ga; t2 = t1`; t3 = ga * a; t4 = t3`; ss5 = ssq(t1 - t2) + ssq(t3 - t4) + ssq(t1 * a - a) + ssq(t3 * ga - ga); print ss5;
Chapter Contents |
Previous |
Next |
Top |
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.