Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 07-23-2019 03:52 PM
(636 views)

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

Paige Miller

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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.

Build your skills. Make connections. Enjoy creative freedom. Maybe change the world. **Registration is now open through August 30th**. Visit the SAS Hackathon homepage.

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.