topic Re: Matrix is singular in SAS/IML Software and Matrix Computations

Matrix is singular

oriti — Mon, 03 Feb 2014 09:49:18 GMT

Dear all,

I have a simulation which in some point it stopped because I get singular matrix.

I read Rick wonderful post about the rareness of singularity ( What is the chance that a random matrix is singular? - The DO Loop)

so it really puzzled me that I get singular matrix.

After discussion in the forum

(checking Rick assumption about the var-covar matrix in the program which has one small eigenvalue)

and thinking again about the model I conclude that my problem is technical.

So..

How can I force SAS to inverse a matrix? The chance that the det=0 is zero..

How can I continue with the simulation without that singular matrix? and how can I print the singular matrix or X data?

Thank you all in advance,

I enclose the code

proc iml;

n=10000;

b=j(6,n);

do j=1 to n;

xtl=j(240, 1);

xsl=j(240, 1);

y=j(240, 1);

xtl[1,1]=0.000844;

xsl[1,1] =2.78896;

y[1,1] =-0.064516129;

a=0.01095 ;*first regression*;

k=0.41450;

ka=0.01349;

kb=0.05634;

kc=-0.00228 ;

kd=-40.63877;

c=-0.00678; *second regression*;

e=0.07520;

d=0.94798;

h=-1.73592;

m=-0.00013407; *third regression*;

l=0.00144;

q=-0.00046335;

r=0.99109;

/** specify the mean and covariance of the population errors **/

Mean = {0, 0, 0};

Cov = {0.0047067768 0.0008647583 0.0000199767,

0.0008647583 0.0081750097 0.0000391468,

0.0000199767 0.0000391468 0.0000022146};

NumSamples = 240;

call randseed(0); /** set seed for the RandNormal module **/

U = RandNormal(NumSamples, Mean, Cov);

do i=2 to 240;

Y[i,1]=a+k*y[(i-1),1]+ka*xsl[(i-1),1]+kb*xtl[(i-1),1]+kc*xsl[(i-1),1]##2+kd*xtl[(i-1),1]##2+u[i,1];

xsl[i,1]=c+e*y[(i-1),1]+d*xsl[(i-1),1]+h*xtl[(i-1),1]+u[i,2];

xtl[i,1]=m+l*y[(i-1),1]+q*xsl[(i-1),1]+r*xtl[(i-1),1]+u[i,3];

end;

ylag=lag(y);

xsl2=xsl##2;

xtl2=xtl##2;

/* Step 1: Compute X`X and X`y */

x = j(nrow(y), 1, 1) || ylag || xsl|| xtl || xsl2||xtl2; /* add intercept column */

v = {1}; /* specify rows to exclude */

idx = setdif(1:nrow(x), v); /* start with 1:240, exclude values in v */

xnew = x[idx, ]; /* extract submatrix */

xpx = xnew`* xnew; /* cross-products */

v = {1}; /* specify rows to exclude */

idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */

ynew = y[idx, ]; /* extract submatrix*/

xpy = xnew` * ynew;

xpxi = inv(xpx); /* form inverse crossproducts */

b[,j] = xpxi * xpy; /* solve for parameter estimates */

*b[,j]= solve(xpx, xpy);/* another way which is more efficient*/;

end;

bt=b`;

varNames = "b0":"b5";

experiment="Sim1":"Sim10000";

print bt[colname=varNames

rowname=experiment

label=" Model"];

quit;

Re: Matrix is singular

PaigeMiller — Mon, 03 Feb 2014 14:17:54 GMT

Sure would be helpful to know which step of this long program produced the error message that a matrix was singular.

Your covariance matrix has a determinant that is close to zero, that may be part of the problem.

Re: Matrix is singular

Rick_SAS — Mon, 03 Feb 2014 14:30:49 GMT

You are misreading my blog post. In my post I generated matrices with elements drawn from a uniform or normal distribution, and showed that the (mathematical) probability of getting a singular matrix is zero (although in finite-precision arithmetic the probability in nonzero).

That is not what you are doing here. You are simulating a timne series model and using a (nearly singular) covariance matrix to generate correlated random errors for the sinulation.

Your problem is blowing up because your time series model is diverging to infinity. IFirst, use a nonzero seed so that you can reproduce the results consistently:

call randseed(1);

Next, insert the following code before you compute the inverse:

if abs(det(xpx))<1e-6 then do;
   z = ynew || xnew;
   create SINGULAR from z[c={Y Int Ylag xsl xtl xsl2 xtl2}];
   append from z;
   close SINGULAR;

   create U from U[c={U1 U2 U3}];
   append from U;
   close U;
end;

If you look at the SINGULAR data set that is created, you will see that the Y, YLag, XSL, and XTL variables are all marching off to infinity.

I don't know how you want to deal with this problem in your model. If you want to ignore the trajectories that diverge, you can use IF(abs(DET(xpx)))>0 to protect the code that computes the parameter estimates. You will have missing values for the bad trajectories.

Re: Matrix is singular

oriti — Tue, 04 Feb 2014 10:03:30 GMT

Thank you Rick for your detailed answer!

I really appreciate your help!

It's really interesting and helpful what you wrote.

When I looked at 'SINGULAR' with your help..I saw that the last 20 simulated observations were marching off to infinity,

So..when I simulated 220 instead of 240 obs everything went well..I got 10000 results as I wanted.

Well, I have to think why and how just the last 20 out of 240 simulated results caused the problem?

Is it technical problem?

I'll appreciate your opinion.

Thank you,

Orit

Re: Matrix is singular

Rick_SAS — Tue, 04 Feb 2014 15:16:05 GMT

It's not a technical problem, it is a mathematical problem. An iterated dynamical system will have a deterministic long-term behavior (bounded or diverging). When you add a stochastic component, you introduce a probability that the system will diverge, even if the deterministic system remains bounded. You can read the main ideas here: Convergence or divergence? A simple iteration with a random component - The DO Loop

A simpler model is a random walk game. Let's play a game where you flip a coin 20 times. You move 1 meter forward whenever the coin is heads; you move 1 meter backwards whenever the coin is tails. Now suppose you begin each game 10 meters aways from the edge of a cliff which has a deadly 1000-foot drop. You will probably be able to play the game for a long time, but the more you play the more you risk encountering an unlucky string of heads that plunge you off the cliff!

Re: Matrix is singular

oriti — Tue, 04 Feb 2014 19:37:58 GMT

Really interesting!

Thank you Rick!