title 'IMLPROJ: Projections'; proc iml; *-- Define a function for a projection of y on X; start proj(y,X) global(P); reset noprint; if nrow(X)=1 then X = t(X); XPX = t(X) * X; P = X * ginv(XPX) * t(X); reset print; return( P * y ); finish; *-- Define a length function (LEN); start len(X); return (sqrt( X[##,] )); finish; reset print log fuzz fw=6; *PROJECTIONS (BOCK: EX.2.6-1); X= {1 1, 1 1, 1 -1, 1 -1}; Y= t(1:4); * 1. Project y on X[,1] and X[,2] separately; r = proj(y,X[,1]); r = proj(y,X[,2]); * 2. Resolve y into 2 orthogonal components, u and v such that y = u + v and u * v = 0; u = proj(y,X); * 3. the projection matrix is P: ; print P; V=(I(4) - P) * Y; * 4. u and v are orthogonal, and sum to Y: ; r = t(u) * v; y = u + v; * 5. the sq-lengths sum to y's; r = len(Y)##2; r = len(U)##2; r = len(V)##2; ; * 6. P and (I-P) are idempotent matrices; * --> P = P * P; * --> det(P)= det(I-P)=0 (not of full rank); * --> nrow(P) = rank(P) + rank(I-P); r = P * P; r = DET(P); r = echelon(P); r = echelon(I(4) - P); page; /*------------------------------------* | The Gram-Schmidt process | *------------------------------------*/ *-- MAKES SUCCESSIVE COLUMNS OF A MATRIX ORTHOGONAL.; X= {1 2 1 6, -1 0 1 2, 1 4 3 3, 0 1 1 0, 1 2 1 0}; *-- get lengths of each column; r = len(X); XPX=t(X) * X; *-- The diagonal elements of X'X are the squared column lengths; r = vecdiag(XPX); *-- but X is singular - det(XPX) = 0; r = det(XPX); ; * the Gram-Schmidt process finds new matrix, Z, which spans same space.; Z =J(nrow(X), ncol(X), 0); Z[,1] = X[,1]; Z[,2] = X[,2] - proj(X[,2], Z[,1]); * the next col is linearly dependent since orthog. proj = 0; Z[,3] = X[,3] - proj(X[,3], Z[,1]) - proj(X[,3], Z[,2]); Z[,4] = X[,4] - proj(X[,4], Z[,1]) - proj(X[,4], Z[,2]); *-- now, all cols of Z are orthogonal; ZPZ = t(Z) * Z; *-- but they differ in length; r = len(Z); *-- delete column 3 (0 vector); Z = Z[,{1 2 4}]; *-- and normalize the columns to make lengths = 1; Z = Z * diag( 1 / len(Z) ) ; *--the new columns of y span the same space, but are orthogonal and = length; r = t(Z) * Z; *-- The steps above are carried out by the GSORTH routine; * Z is orthonormal matrix, same size as X; * T is an upper-triangular matrix, such that X = Z*T; * flag is 0 if columns of X are linearly independent, else 1; * rank deficiency -> 0 columns of Z and 0 rows of T; call gsorth(Z,T,flag,X); *-- Where orthogonal polynomials come from; C = {1 1 1, 1 2 4, 1 3 9, 1 4 16}; call gsorth(Z,T,flag,C); quit;