Recently, I am working on coding in SAS for a set of regularized regressions and need to compute trace of the projection matrix:
$$ S=X(X'X + \lambda I)^{-1}X' $$.
Wikipedia has a well written introduction to Trace @
here.
To obtain the inverse of matrix (X'X + \lambda I) in SAS/STAT, there are multiple ways:
1. Build the SSCP matrix first, then inverse it using the following methods:
1.1. Using SVD based method, shown @
here but not demonstrated in this post;
1.2. Using PROC REG, shown @
here , and also demonstrated in the code below;
1.3. Using SWEEP operator, shown @
here under Method 0 and Method 1;
2. Since (X'X + \lambda I) is the SSCP for a Ridge Regression with ridge parameter=\lambda, it is handy to directly use the ODS OUTPUT InvXPX= statement to obtain the inversed matrix, when X is appended by a diagonal matrix of \lambda*I_{p, p}. SAS code is demonstrated below;
Trace of the project matrix S is a key concept in modern regression analysis. For example, the effective degree of freedom of the model in a regularized linear regression is trace(S)/N, see [1] for details. For another example, in approximating leave-one-out-cross-validation using GCV, trace(S)/N is a key component (formula 7.52 of ~[1]).
Check the reference and the cited works therein for more information.
Ref:
1. Hastie et al., The Elements of Statistical Learning, 2nd Ed.
/* Use the SocEcon example data created in
Example A.1: A TYPE=CORR Data Set Produced by PROC CORR
on page 8153 of SAS/STAT 9.22 User Doc
*/
data SocEcon;
input Pop School Employ Services House;
datalines;
5700 12.8 2500 270 25000
1000 10.9 600 10 10000
3400 8.8 1000 10 9000
3800 13.6 1700 140 25000
4000 12.8 1600 140 25000
8200 8.3 2600 60 12000
1200 11.4 400 10 16000
9100 11.5 3300 60 14000
9900 12.5 3400 180 18000
9600 13.7 3600 390 25000
9600 9.6 3300 80 12000
9400 11.4 4000 100 13000
;
run;
%let depvar=HOUSE;
%let covars=pop school employ services;
%let lambda=1;
/* Below is the way to obtain trace(S), where S is the project matrix in a (regularized) linar regression.
For further information, check pp.68, pp.153 of Elements of Statistical Learning,2nd Ed.
*/
/* For details about TYPE=SSCP special SAS data, consult:
Appendix A: Special SAS Data Sets, pp.8159 of SAS/STAT 9.22 User's Guide
*/
proc corr data=SocEcon sscp out=xtx(where=(_TYPE_='SSCP')) noprint;
var &covars;
run;
data xtx2;
set xtx;
array _n{*} _numeric_;
array _i{*} i1-i5 (5*0);
do j=1 to 5;
if j=_n_ then _i[_n_]=λ
else _i[j]=0;
end;
_n[_n_]=_n[_n_]+1;
drop j _TYPE_ _NAME_;
run;
/* Obtain the inverse of (XTX+\lambda*I)
Note that we explicitly specified Intercept term in the
covariate list and fit a model without implicit intercept
term in the model.
*/
proc reg data=xtx2
outest=S0(type=SSCP
drop=i1-i5 _MODEL_ _DEPVAR_ _RMSE_)
singular=1E-17;
model i1-i5 = Intercept &covars / noint noprint;
run;quit;
data S0;
set S0;
length _NAME_ $8;
_NAME_=cats('X_', _n_);
run;
proc score data=SocEcon score=S0 out=XS0(keep=X_:) type=parms;
var &covars;
run;
data XS0X;
merge XS0 X;
array _x1{*} X_:;
array _x0{*} intercept pop school employ services;
do i=1 to dim(_X1);
_x1[i]=_x1[i]*_x0[i];
end;
rowSum=sum(of _x1[*]);
keep rowSum;
run;
proc means data=XS0X noprint;
var rowSum;
output out=trace sum(rowSum)=Trace;
run;
Verify the result using R:
> socecon<-read.csv('c:/socecon.csv', header=T)
> x<-as.matrix(cbind(1, socecon[,-5]))
> xtx<-t(x)%*%x
> phi<-xtx+diag(rep(1, 5))
>
> # method 1.
> S<-x%*%solve(phi)*x
> sum(S)
[1] 4.077865
> # method 2.
> S<-(x%*%solve(phi))%*%t(x)
> sum(diag(S))
[1] 4.077865
>
Of course, except for method 1.2 shown above, we can also use Method 2 mentioned above, and obtain the same inverse matrix:
data SocEcon2 /view=SocEcon2;
set x end=eof;
array x{5} Intercept &covars (5*0);
Intercept=1;
output;
if eof then do;
do j=1 to dim(x);
x[j]=λ
output;
drop j;
x[j]=0;
end;
end;
run;
ods select none;
ods output InvXPX=S1;
proc reg data=SocEcon2 singular=1E-17;
model y = Intercept &covars /noint i;
run; quit;
ods select all;
As a side point, it is also of interests to compare the computational performance of both illustrated methods. SVD-based approach and SWEEP operator based approach are omitted.
Using a SAS data set of 1001 covariates (including Intercept term) and 1E6 observations, total 7.6GB on a windows PC, the test was done on two dramatically different machines: one WindowsXP laptop with mediocre HDD and a 2-core T7300 CPU; the other one is a high-end Server running Linux64 with Disk Array and fast CPUs totaled 16 cores.
On the PC,:
-> Using method 1.2, the time used decomposition is listed below:
PROC CORR: real time: 25:47.72; CPU time: 15:26.15
DATA step on XTX: real time: 2.28 sec; CPU time: 0.07 sec
PROC REG : real time: 15.45 sec; CPU time: 27.31 sec;
Total real time: 26:05.45
-> Using method 2, the time used is:
PROC REG: real time: 48:28.44; CPU time: 1:16:46.41
On the server:
-> Using method 1.2, the time used decomposition is listed below:
PROC CORR: real time: 5:40.61; CPU time: 5:40.58
DATA step on XTX: real time: 0.05 sec; CPU time: 0.05 sec
PROC REG: real time 1.71 sec; CPU time: 5.27 sec
Total real time: 5:42.37
-> Using method 2, the time used is:
PROC REG: real time: 6:01.46; CPU time: 19.13.49
The performance is summarized below:
