sas/iml矩阵算术

这样的资源网络上是有的，不过大都是英文，这里姑且进行一些拾牙慧，用汉语进行标注。

1.行列式

*行列式特征;
proc iml;  reset print log;
   A = {3 1, 2 4};
   r = det(A);     *行列式的值;
 
   r = det(A[{2 1},]);*交换行列式的行行列式变号;
   r = det(A[,{2 1}]);*交换行列式的列行列式变号;
 
   r = det( t(A) );*行列式与它的转置行列式相等;
 
 *行列式的某行/列所有元素乘以K等于K乘以该行列式;
   r = diag({3 1}) * A;
   r = det( diag({3 1})  * A);
 
 *矩阵乘以K的行列式，等于K的n方乘以该矩阵行列式;
   r = det(3 # A);
 
 *方阵乘积的行列式等于两方阵行列式的乘积;
   B={ 4 2, 3 5};
   r = det(B);
   r = det(A * B);
   r = det(A) # det(B);
 
 *行列式中行列成比例，则行列式值为0;
   C={1 5, 2 6, 4 4};
   C=C || C[,1];
   r = det(C);
 
 *把行列式某一列/行的各元素乘以同一数后加到另一列/行
   对应元素上去，行列式不变;
   A[2,] = A[2,] - 2#A[1,];
   r = det(A);
quit;

*行列式求值;
proc iml;
   reset print log;
  *计算行列式的余子式;
start cofactor(A,i,j);
   reset noprint;
   size = nrow(A);                *矩阵维度;
   rows = loc( (1:size) ^= i );   *删除第i行，loc函数用来生成由非零元素位置组成的行向量;
   cols = loc( (1:size) ^= j );   *删除第j列;
   minor = det( A[rows,cols] );
   reset print;
   return ( ((-1)##(i+j)) # minor );
finish;
 
 *行列式按余子式展开;
    A ={4 2 1, 5 6 7, 1 0 3};
    r = det(A);
    *按照第一行的余子式展开cofactors of row 1 elements;
    print (cofactor(A,1,1)) (A[{2 3},{2 3}]);
    print (cofactor(A,1,2)) (A[{2 3},{1 3}]);
    print (cofactor(A,1,3)) (A[{2 3},{1 2}]);
    *行列式等于行元素与其余子式的乘积和
	det = product of row with cofactors;
    r = A[1,] * {18, -8, -6};
 
    *行列式按照高斯消元法求值;
    M = {2 3 1 2,  4 2 3 4,  1 4 2 2,  3 1 0 1};
    d = det(M);
 
    *以M[1,1]为主元素进行消元，行列式等于各主元素的乘积;
    D=M[1,1];
    *首先消掉第一行,第一列同时被消掉;
    M[1,] = M[1,] / M[1,1];
    M = M - M[,1] *  M[1,];
    *剩下2、3、4行列组成新的行列式;
    M = M[{2 3 4},{2 3 4}];
    *主元素进行累积;
    D = D # M[1,1];
 
    *重复上述消元过程;
    M[1,] = M[1,] / M[1,1];
    M = M - M[,1] *  M[1,];
    M = M[{2 3},{2 3}];
    D = D # M[1,1];
 
    *再一次消元，等到行列式值d = det(m);
    M[1,] = M[1,] / M[1,1];
    M = M - M[,1] *  M[1,];
    M = M[2,2];
    D = D # M[1,1];
quit;

 

2.矩阵

a.基本的矩阵

*基本的矩阵;
options nodate;
proc iml;
   reset print log;
   *一般矩阵;
   X = {1 2 3, 4 0 6, 9 1 4, 6 2 4};
 
   print (nrow(X)) (ncol(X));
   rowvec = { 1 2 3};
   colvec = { 1, 4, 6};
   *矩阵转置;
   xt = t(x);
   xt = t(x);
   *方阵;
   A = {5 -1 3, -1 4 7, 3 3 -4};
   *对角线元素组成的列向量:vecdiag函数;
   d = vecdiag(A);
   *对角阵;
   d = diag({1 4 9});
   *矩阵的迹:对角元素之和;
   t = trace(A);
   *对称阵;
   sym = ( A = t(A) );
   sym = all( A = t(A) );
   *方阵性质;
   B = ( A + t(A) ) / 2;
   sym = all( B = t(B));
   *上三角阵;
   upper = {-5 2 7, 0 -2 4, 0 0 3};
   *下三角阵;
   lower = t(upper);
   *单位阵;
   I = I(3);
   *纯量矩阵;
   S = diag({10 10 10});
   *单位元素1组成的矩阵;
   J = J(2,5);
   J = J(3,1);
   *零元素组成的矩阵;
   Z = J(2,5,0);
quit;

b.矩阵的基本运算

options nodate;
proc iml;
   reset print log;
   *加法和减法;
   a = {1 2 3, 4 5 6};
   b = {5 1 0, 3 -2 1};
   c = {2 2 2, 1 1 1};
   *对应元素相加;
   d = a + b;
   d = a + {1 3, 3 10};
   *对应元素相减-- matrices are subtracted elementwise;
   d = a - b;
   *负数看成减法;
   d = -a;
   *矩阵加减的性质;
   *   a. 交换性: A + B = B + A;
   print ( A + B ) ( B + A );
   equal = all( (A + B) = (B + A) );
   *   b. 联合性: A + (B+C) = (A+B) + C;
   print ( (A+B) + C ) ( A + (B+C) );
  
   *矩阵与常数乘积;
   d = 3 # A;
   * 乘法分配律;
   print ( 3 # (A + B) )  ( (3#A) + (3#B) );
   *矩阵对应元素乘除法;
   d = A # B;
   d = A / 2;
quit;

c.逆矩阵性质

proc iml; reset print log fuzz;
   A = {5 1 0,  3 -1 2,  4 0 -1};
   r = det(A);
 
   *行列式不等于0，矩阵逆存在;
   AI = inv(A);
    *矩阵预期逆乘积为单位阵;
   r = AI * A;
    *矩阵逆的逆矩阵等于该矩阵;
   r = inv(AI);
    *对称阵的逆矩阵是对称的;
   r = inv( t(A) );
   B = {4 2 2, 2 3 1, 2 1 3};
   r = inv(B);
   r = inv(t(B));
   *对角阵的逆矩阵等于其各元素倒数组成的对角阵;
   D =diag({ 1 2 4 });
   r = inv(D);
   r = diag(1 / {1 2 4});
   *转置矩阵的逆等于逆的转置;
   r = inv( t(A) );
   r = t( inv(A) );
   *常量与矩阵乘积的逆等于常量倒数乘以矩阵逆 ;
   r = inv(5#A);
   r = (1/5)#inv(A);
   *矩阵乘积的逆等于逆的乘积：inv(A * B) = inv(B) * inv(A);
   B = {1 2 3, 1 3 2, 2 4 1};
   r = A * B;
   r = inv(A * B);
   r = (inv(B)) * (inv(A));
quit;

*行列式为0的矩阵的逆的情况;
proc iml;
   reset print log fuzz fw=6;
*构造一个降序阵;
   A = {4 4 -2, 4 4 -2, -2 -2 10};
   r = det(A);    *矩阵的行列式;
   r = echelon(A);*初等行运算进行降为，A矩阵秩为2，没有一般的逆矩阵;
   r = inv(A);
   AI = ginv(A);*ginv函数产生广义逆矩阵;
*广义逆矩阵的性质;
   r = A * AI * A;
   r = AI * A * AI;
*A*AI和AI*A 都是对称的,但二者均不等于单位阵;
   r = A * AI;
   r = AI * A;
*对于矩形矩阵，如果(A'A)是(A'A)的广义逆矩阵，那么inv(A'A)*A'是矩阵A的广义逆矩阵;
   A = J(4,1) || {1 0, 1 0, 0 1, 0 1};
   AA= t(A) *  A;
   AAI= ginv(AA);
*矩阵A的广义逆矩阵是AAI * A';
   AI = AAI  *  t(A);
   r = A * AI * A;
   r = AI * A * AI;
quit;

d.矩阵的秩与向量相关性

*矩阵的秩;
proc iml;
   *1.定义函数R，用来计算矩阵A的秩;
start r(A);
    reset noprint;
   *矩阵秩行初等变换后非零行的数目;
   e = echelon(A);
   do i=1 to nrow(e);   *找到不是所有元素都为0的行;
      if any( e[i,] <> 0 ) then rows = rows || i;
      end;
   e = e[rows,];        *保留非0行;
   do i=1 to ncol(e);   *找到非零列;
      if any( e[,i] <> 0 ) then cols = cols || i;
      end;
   e = e[,cols];        *保留非零列;
   reset print;
   return( min( nrow(e), ncol(e)) );
   finish;
 
   reset print log fuzz;
* 秩与线性相关;
X1 = {1 3 5};
X2 = {0 1 2};
X3 = {-1 4 9};
X4 = {2 2 2};
comb = (2#X1) + (-4#X2) + (0#X3) + -1#X4;
* X3, X4 是 X1, X2的线性组合;
print (t(X3)) '=' (t(-X1)) '+' (t(7#X2));
print (t(X4)) '=' (t(2#X1)) '+' (t(-4#X2));
*  x1, x2, x3, x4中只有两个是独立的;
X = t(X1) || t(X2) || t(X3)  || t(X4);*X是X1,X2，X3，X4的联合矩阵;
r = r(X);                            *调用定义函数，计算X的秩;

r = echelon(X); *X进行行初等变换;

X[3,4] = 4;*变换X的三行四列的元素;
r = r(X);
r = echelon(X);
skip 3;   *log中空3行;
* -----------------------------------------------------;
* 2. 初等行运算计算矩阵秩;
A = {4 1 8, 5 2 7,  5 1 11,  8 3 12};
SAVE = A;
* ROW4 - 2#ROW 1;
A[4,] = A[4,] - 2#A[1,];
*  Each elementary row operation is equivalent to premultiplication
   by an identity matrix with the same operation applied to its rows;
T = I(4);
T[4,] = T[4,] - 2#T[1,];
A = T * SAVE;
   * ROW 3 - ROW 2;
A[3,] = A[3,] - A[2,];
   * .25 # ROW 1;
A[1,] = .25 # A[1,];
   * ROW 2 - 5#ROW 1;
A[2,] = A[2,] - 5#A[1,];
   * ROW 4 + 1#ROW 3;
A[4,] = A[4,] + A[3,];
   * 1.33#ROW 2;
A[2,] = (4/3) # A[2,];
   * ROW 2 + ROW 3;
A[2,] = A[2,] + A[3,];
   *RANK = number of nonzero rows/cols;
r = r(A);
quit;

 

3.解方程组，函数R是求秩定义的函数(见上例)

*解方程组;
proc iml;
   reset print log fuzz fw=5;
*1.三个等式三个未知数（齐次方程组）;
A= {1 1 -4, 1 -2 1, 1 1 1};
b= {2, 1, 0};
xx = t('X1' : 'X3');
print A '*' xx '=' b;
print (r(A)) (r(A||b));* 齐次有：r(A) = r(A b);
*利用inv() or solve()函数解方程;
x = inv(A) * b;
x = solve(A,b);
print A ' * ' x '=' (A * x) '=' b;
r = echelon(A||b);*(A || b)初等行变换Echelon显示解放方程过程;

*2.四个等式三个未知数的方程组（齐次）;
A= {1 1 -4, 1 -2 1, 1 1 1, 2 -1 2};
b= {2, 1, 0, 1};
print (r(A)) (r(A||b));*检验齐次：r(A) = r(A b);
*方程比未知数多,可以考虑广义逆矩阵ginv;
x = ginv(A) * b;
print A ' * ' x '=' (A * x) '=' b;
r = echelon(A||b);*初等行变换Echelon显示解放方程过程;

*3.三个等式三个未知数方程组（非齐次）;
A={1 3 1, 1 -2 -2, 2 1 -1};
b={2,3,6};
print (r(A)) (r(A || b));*检验非齐次： r(A) < r(A b);
r = echelon(A||b);*初等行变换过程;
x = inv(A) * b;*由于A是奇异阵所有函数inv()失效;

x = ginv(A) * b;*考虑用广义逆矩阵函数解方程;
print A '*' x ' = ' (A * x);

*给定b向量，解方程组;
b = {2,3,5};
r = r(A||b);
r = echelon(A||b);

x = ginv(A) * b;
print A '*' x ' = ' (A * x) '=' b;
*等价解法2;
A11 = A[{1 2},{1 2}];
A12 = A[{1 2}, 3];
b1  = b[{1 2},];
x2  = -1.1;

x1 = inv(A11) * b1 - inv(A11) * A12 * x2;* 两个独立的未知变量的解;
quit;

 

4.特征值与特征向量

proc iml;
   reset print log fuzz fw=5;
A={13 -4 2, -4 13 -2, 2 -2 10};
call eigen(values, vectors, A); *call eigen语句得到特征向量vectors，特征值value;
*特征值性质;
r = sum(values);*矩阵的迹等于所有特征值和;
r = sum(vecdiag(A));
*矩阵对应元素平方和等于特征值的平方和ssq(A) = ssq(eigenvalues);
print (ssq(A)) '=' (ssq(values));
print (det(A)) '=' (values[#]);*矩阵行列式等于特征值乘积;
*矩阵的秩等于非零特征值的个数;
r = echelon(A);
rank = (((echelon(A)^=0)[,+]) ^=0)[+,];
rank = sum(values ^= 0);
*逆矩阵的特征值等于矩阵特征值的倒数,但有相同的特征向量;
AI = inv(A);
call eigen(val, vec, AI);
* similar relation for other powers: roots(A##n) = (roots(A) ## n);
call eigen(val, vec, A*A);
call eigen(val, vec, A*A*A);
page;

*特征根与特征向量;
A={13 -4 2, -4 13 -2, 2 -2 10};
call eigen(L, V, A);
*1.向量与特征根的乘积等于特征根与其对应的特征向量的乘积;
r = A * V[,1];
r = L[1] # V[,1];
print A '*' (V[,1]) '=' (A*V[,1]) '=' (L[1]) '#' (V[,1]);
r = A * V;
r = V  *  diag(L);
print A '*' V  '='  V '*' (diag(L)) ;
*利用特征值对矩阵进行对角化: V'AV = diagonal;
r = t(V) * A * V;
*2.特征向量系正交的，且其与转置阵的乘积为单位阵V'V=I;
r = t(V) * V;
*3.矩阵的谱分析: A = sum of rank 1 products;
A1 = L[1] # V[,1] * t(V[,1]);
A2 = L[2] # V[,2] * t(V[,2]);
A3 = L[3] # V[,3] * t(V[,3]);
r = A1 + A2 + A3;
r = ssq(A);
r = ssq(A1) || ssq(A2) || ssq(A3);
*-- Each root is sum of squares accounted for;
r = L##2;
*--The first i roots and vectors give a rank i approximation;
r = echelon(A1+A2);
r = ssq(A1+A2);
quit;

 参考：《线性代数》同济大学第四版

代码摘自http://www.yorku.ca/health/psyc/