BookmarkSubscribeRSS Feed
Esta_Stoop
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;
	end;
	return(B);
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;
	end;
return(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;
quit;
2 REPLIES 2
PaigeMiller
Diamond | Level 26

See if this document from @Rick_SAS helps:

https://blogs.sas.com/content/iml/2018/11/21/generalized-inverses-for-matrices.html

 

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
Rick_SAS
SAS Super FREQ

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 matrix.to develop and debug your program.

 

SAS INNOVATE 2024

Innovate_SAS_Blue.png

Registration is open! SAS is returning to Vegas for an AI and analytics experience like no other! Whether you're an executive, manager, end user or SAS partner, SAS Innovate is designed for everyone on your team. Register for just $495 by 12/31/2023.

If you are interested in speaking, there is still time to submit a session idea. More details are posted on the website. 

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
  • 414 views
  • 0 likes
  • 3 in conversation