/* This code is stored as matrix.sas It uses the IML procedure in SAS. The start command gives a name to a set (or module) of commands. The finish command specifies the end of the module. This module can be run anywhere in your SAS code by applying the run command to the name of the module.*/ proc iml; /* Nmae the module matrix1 */ start matrix1; /* Add and subtract matices */ a = {3 6 , 2 1}; b = {7 -4, -3 2}; c = a+b; d = a-b; print a, b, c, d; /* Compute the sum of all elements in a matrix */ suma=sum(a); print a suma; /* Compute row sums */ rowsuma = a[ ,+]; print 'row sums' rowsuma; /* Compute column sums */ colsuma = a[+, ]; print 'columns sums' colsuma; /* Select the first row of a */ row1 = a[1, ]; print 'row1' row1; /* Select the second column of a */ col2 = a[ ,2]; print 'col2' col2; /* Select the element in the first row and the second column of a */ a12 = a[1,2]; print a12; /* Multiplication by a scalar */ c = {2 -1 3, 0 4 -2 }; d = 2*c; print c, d; /* Transpose of a matrix */ ct = t(c); print c, ct; /* Matrix multiplication */ a = {3 0 -2, 1 -1 4}; b = {1 1, 1 2, 1 3}; c = a*b; print,,, a, b, c; /* Here are two ways to compute the inner product of two vectors. The second method is more efficient */ x = {1, 7, -6, 4}; y = {2, -2, 1, 5}; xy1 = t(x)*y; xy2 = sum(x#y); print,,,, x, y, xy1, xy2; /* Length of a vector */ ynorm = sqrt(sum(y#y)); print ynorm; /* Number of rows in the column vector vector y */ ny = nrow(y); print y ny; /* Elementwise multiplication */ a = { 3 6, 2 1}; b = { 7 -4, -3 2}; ab = a#b; print,,, a, b, ab; /* Kronecker Product */ a = {2 4, 0 -2, 3 -1}; b = {5 3, 2 1}; ab = a@b; print,,, a, b, ab; /* Create an identity matrix */ a = i(4); print a; a = i(3); print a; /* Create a matrix where every element is the same. The first argument is the number of rows, the second argument is the number of columns, the last argument is the value of each element */ a = j(2,3,4.2); print a; /* Two methods of computing the trace of a matrix */ w = { 1 2 3, 4 5 6, 7 8 10}; trw1 = trace(w); trw2 = sum(diag(w)); print,,, w trw1 trw2; /* Inverse of a matrix */ winv = inv(w); prod = w*winv; print,,, w, winv, prod; /* Determinant of a matrix */ detw = det(w); print,,, w detw; x1 = { 1 2 3, 4 5 6, 7 8 9}; detx1 = det(x1); print,,, x1 detx1; /* IML does not have a function that computes the rank of a matrix. You can use the following which makes use of the generalized inverse function */ a = {1 1 1, 2 5 -1, 0 1 -1}; ranka=round(trace(ginv(a)*a)); print a ranka; /* You can also compute the QR decomposition of the matrix. Here a = q*r and the number of linearly dependent columns of a is stored in lindep */ CALL QR( q, r, piv, lindep, a); ranka = ncol(a)-lindep; print a, ranka, q, r; /* Other examples */ c = {1 1 1, 2 5 -1, 0 1 1}; rankc=round(trace(ginv(c)*c)); print c rankc; B = {1 1 0 0, 1 1 0 0, 1 0 1 0, 1 0 1 0, 1 0 0 1, 1 0 0 1}; rankB=round(trace(ginv(B)*B)); print B rankB; /* Note the rank function in IML is related to sorting. It computes the ranks of numbers in a matrix */ c = {1.2 5.1, 3.5 9.8}; rc = rank(c); print,,,,, c, rc; /* Eigenvalues and eigenvectors: In call eigen(e, u, a) the eigenvalues of a are stored in the vector e and the eigenvectors of a are stored in u. */ a = {1.96 0.72, 0.72 1.54}; call eigen( e, u, a); print,,, a,,, 'Eigenvalues and eigenvectors', e u; /* Singular value decomposition: The singular values are stored in q in the following example and a = u*diag(q)*t(v) */ a = { 2 0 1 1, 0 2 1 1, 1 1 1 1}; call svd(u, q, v, a); ar = u*diag(q)*t(v); print a, 'Singular value decomposition', u, q, v, ar; /* An example where the singular values are the eigenvalues */ a = {1.96 0.72, 0.72 1.54}; call svd( u, q, v, a); print a, 'Singular value decomposition', u q v; /* Eigenvalues of a square symmetric positive definite matrix are all positive. They match the singular values. */ a = {4 2 -1, 2 6 -4, -1 -4 9}; call eigen( e, u, a); print,,, a,,, 'Eigenvalues and eigenvectors', e, u; call svd( u, q, v, a); print 'Singular value decomposition', u, q, v; /* Eigenvalues of a square symmetric matrix that is not positive definite. How do they match the singular values? */ a = {4 2 -1, 2 6 -4, -1 -4 -9}; call eigen( e, u, a); print,,, a,,, 'Eigenvalues and eigenvectors', e, u; call svd( u, q, v, a); print 'Singular value decomposition', u, q, v; /* Inverse of a non-singular matrix */ a = {1.96 0.72, 0.72 1.54}; ainv = inv(a); ident = a*ainv; print,, a, ainv, ident; /* Use the spectral decomposition to compute an inverse square root matrix */ a = {1.96 0.72, 0.72 1.54}; call eigen( e, u, a); print,,, a,,, 'Eigenvalues and eigenvectors', e, u; ainvm2 = u*inv(diag(sqrt(e)))*t(u); ainv = ainvm2*ainvm2; ident = ainvm2*a*ainvm2; print ainvm2, ainv, ident; /* Solve a set of linear equations Ab=x These two methods provide the same results here, but the solve function is more efficient and more accurate for larger sets of equations */ A = {1.96 0.72, 0.72 1.54}; x = {1, 1}; b1 = inv(A)*x; b2 = solve(A,x); print,, A, x, 'Solutions', b1, b2; finish; /* Execute the set of commands in the module named matrix1. The results requested by the print commands appear in the output window */ run matrix1;