Statistical programming, matrix languages, and more

Matrix is singular

Accepted Solution Solved
Reply
Contributor
Posts: 39
Accepted Solution

Matrix is singular

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;


Accepted Solutions
Solution
‎02-03-2014 09:30 AM
SAS Super FREQ
Posts: 3,234

Re: Matrix is singular

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.

View solution in original post


All Replies
Trusted Advisor
Posts: 1,294

Re: Matrix is singular

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.

Solution
‎02-03-2014 09:30 AM
SAS Super FREQ
Posts: 3,234

Re: Matrix is singular

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.

Contributor
Posts: 39

Re: Matrix is singular

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

SAS Super FREQ
Posts: 3,234

Re: Matrix is singular

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!

Contributor
Posts: 39

Re: Matrix is singular

Really interesting!

Thank you Rick!

Post a Question
Discussion Stats
  • 5 replies
  • 870 views
  • 6 likes
  • 3 in conversation