suncaiyun
Y=(3 0 0 0) ',X是4阶矩阵，X=(x1, x2, x3,x4),x1=(1 20 60 0)',
x2=(0 1 0 0)',x3=(0 0 1 0)',x4=(0 10 9 0)'，Y与X1-X4都是（4*1）的列向量;
(1) 计算矩阵X'YY'X 的特征值为0的(k-1)个正交特征向量b(1),b(2),...,b(k-1)。
上面得到的正交特征向量记为b(1),b(2),...,b(k-1)，把这(k-1)个特征向量作为列向量
而成一个k*(k-1) 阶的矩阵B=(b(1),b(2),...,b(k-1))。
(2) 计算矩阵 B'X'XB的(k-1)个特征值p(1),p(2),...,P(k-1)及其相应的标准特征向量
a(1),a(2),...,a(k-1)。然后选择s个最大特征值，使得它们的和占全部特征值和的
绝大多数（>90%）。这s个最大特征值对应的特征向量被记为a(1),a(2),...,a(s)。
(3) 令T=XBA，其中矩阵A被定义为A=(a(1),a(2),...,a(s))。
（4）计算X(I-BA(T'T)（-1）T'X)=XD，其中，I是同阶单位矩阵。
(说明：前面(T'T)后面的（-1）表示(T'T)的逆矩阵，由于公式不好输入，只好用此法代替，见笑）
邮箱：gmscy@126.com
da_da_mao
in SAS, you can use
proc IML;
y=shape(0,4,1);
y[1,1]=3;
x=shape(0,4,4);
x[1,1]=1;
x[1,2]=20;
...
H=X*Y*t(Y)*t(X);
B=eigen9H);
eiB=t(B)*t(X)*X*B;
...
quit;
just like this, you can work out the result you want.
Hope this helps.
