BookmarkSubscribeRSS Feed
Calcite | Level 5

I am a mathematical statistics masters student at the University of Pretoria, South-Africa.


The SWEEP function in SAS IML works very well to calculate the G1 generalized inverse of a large matrix. I need to calculate the G1 generalized inverse (only satisfying the condition AGA = A, not the Moore-Penrose inverse) of a large matrix using the SWEEP operator. I want to program the sweep operator.

The matrices I need to calculate the G1 generalized inverse of is very large, dimension approximately 3000 x 3000.

I have tried to program the sweep operator to calculate the G1 generalized inverse using the following article as a reference
"A Tutorial on the SWEEP Operator" Author: James H. Goodnight.

My program is not generating the correct G1 generalized inverse.


I also want to make the program efficient.

Can you please guide me and give me advice on how to program the sweep operator?

Is it wise to use a loop in the program?

ginv(A) = SWEEP(A,1:nrow(A))
The sweep operator applied to all the rows/columns of the matrix = G1 generalized inverse
Is there some trick that can help me make my program more efficient? Since I am working with very large matrices.

proc iml;
/*SWEEP operator to calculate G1 - first attempt*/
start sweep_ginv(A);
	n = nrow(A);
	do k = 1 to n;
		B = J(n,n,0);
		print A;
		p = A[k,k];
		B[setdif(1:n,k),k] = -A[setdif(1:n,k),k]/p;
		B[k,setdif(1:n,k)] = A[k,setdif(1:n,k)]/p;
		B[setdif(1:n,k),setdif(1:n,k)] = A[setdif(1:n,k),setdif(1:n,k)]-(1/p)*A[setdif(1:n,k),k]*A[k,setdif(1:n,k)];
		B[k,k] = 1/p;
		print B;
		A = B;
		print A;
finish sweep_ginv;

/*Second attempt*/
start sweep_ginv2(F);
	n = nrow(F);
	range = 1:n;
	do k = 1 to n;
		B = J(n,n,0);
		p = F[k,k];
		B[loc(range ^= k),k] = -F[loc(range ^= k),k]/p;
		B[k,loc(range ^= k)] = F[k,loc(range ^= k)]/p;
		B[loc(range ^= k),loc(range ^= k)] = F[loc(range ^= k),loc(range ^= k)]-(1/p)*F[loc(range ^= k),k]*F[k,loc(range ^= k)];
		B[k,k] = 1/p;
		F = B;
finish sweep_ginv2;

/*Test for small matrices first*/
seed = 1234;
n = 40;
c = J(n, n, seed); 
A = normal(c);

/*Test for generalized inverse
the sum should be zero*/
T = sweep_ginv2(A);
sum = sum(A*T*A - A);
print sum;
Diamond | Level 26

See if this document from @Rick_SAS helps:


In general, when you are trying to do some common mathematical algorithm in PROC IML, Rick has already explained how to do it, so his blog is a good place to start looking for help.

Paige Miller

I'm on vacation, but the SWEEP operator generates an inverse when applied to a SYMMETRIC matrix. Change your example to use a symmetric matrix.


In general, a good debugging strategy is to use a small matrix that you can print during each iteration of an algorithm. Use a 4x4 or 5x5 symmetric develop and debug your program.




Time is running out to save with the early bird rate. Register by Friday, March 1 for just $695 - $100 off the standard rate.


Check out the agenda and get ready for a jam-packed event featuring workshops, super demos, breakout sessions, roundtables, inspiring keynotes and incredible networking events. 


Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 2 replies
  • 3 in conversation