title 'IMLRANK: Matrix rank and linear independence';
proc iml;
   *-- Define a function module to find the rank of a matrix;
start r(A);
    reset noprint;
   *-- rank of a matrix (cant call it rank, since that name
       is used for sorting);
   *-- rank = number of nonzero rows/cols in echelon form;
   e = echelon(A);
   do i=1 to nrow(e);   *-- Find rows which are not all zero;
      if any( e[i,] <> 0 ) then rows = rows || i;
      end;
   e = e[rows,];        *-- Keep only non-zero rows;
   do i=1 to ncol(e);   *-- Find cols which are not all zero;
      if any( e[,i] <> 0 ) then cols = cols || i;
      end;
   e = e[,cols];        *-- Keep only non-zero cols;
   reset print;
   return( min( nrow(e), ncol(e)) );
   finish;
 
   reset print log fuzz;
* 1. MATRIX RANK AND LINEAR INDEPENDENCE;
*    A set of vectors is linearly DEPENDENT if there are scalars;
*    (not all =0) which give a linear combination = zero vector.;
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 are linear combinations of X1, X2;
print (t(X3)) '=' (t(-X1)) '+' (t(7#X2));
print (t(X4)) '=' (t(2#X1)) '+' (t(-4#X2));
* only two of x1, x2, x3, x4 are linearly independent;
X = t(X1) || t(X2) || t(X3)  || t(X4);
r = r(X);
* rank = number of nonzero rows in echelon form;
r = echelon(X);
* change last element so X4 not dependent;
X[3,4] = 4;
r = r(X);
r = echelon(X);
skip 3;
* -----------------------------------------------------;
* 2. MATRIX RANK BY ELEMENTARY ROW OPERATIONS;
*    Use elementary row ops to reduce matrix to echelon
*    form. Zero rows/cols do not contribute to rank.;
* -----------------------------------------------------;
A = {4 1 8, 5 2 7,  5 1 11,  8 3 12};
SAVE = A;
* ROW 4 - 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;