ERROR: Overflow error in DET.

richard_hu2003 · Posted 07-15-2011 11:52 AM

I got the error message below. This is from IML which is within a big macro "LARGE". I have a do loop outside "LARGE".

1. What does the overflow error mean? Does it mean that the determinant of the matrix is too large?

2. In the current code, the IML will abort if singularity occurs (or determinant is zero). If I want the IML to goto exit whenever it encounters ANY types of errors, how should I code it in this situation? I want the IML to goto exit and continue to execute for next iteration WITHOUT setting OBS=0.

I also include code for your convenience. Thanks.

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Error message:

NOTE: IML Ready
NOTE: Module MODULE1 defined.
NOTE: Module MODULE2 defined.
ERROR: Overflow error in DET.
ERROR: Execution error as noted previously. (rc=100)

operation : DET at line 253640 column 4
operands : ym2Sum
ym2Sum 145 rows 145 cols (numeric)

statement : ASSIGN at line 253638 column 247
traceback : module MODULE1 at line 253638 column 247

NOTE: Paused in module MODULE1.
NOTE: The data set OUT.ALPHA has 145 observations and 1 variables.
NOTE: The data set OUT.YM2SUM has 145 observations and 145 variables.
NOTE: The data set OUT.YMPMSUM has 145 observations and 1 variables.
NOTE: The data set OUT.INV_YM2SUM has 145 observations and 145 variables.
NOTE: The data set OUT.C17201_1_IML_IDX has 145 observations and 1 variables.
NOTE: Exiting IML.
NOTE: 32 workspace compresses.
NOTE: The SAS System stopped processing this step because of errors.
NOTE: SAS set option OBS=0 and will continue to check statements. This may cause NOTE: No observations in data set.
NOTE: PROCEDURE IML used (Total process time):
real time 0.13 seconds
cpu time 0.10 seconds

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Code:

/* IML procedure */

%LET IsError = 0;

proc iml;

start module1; /* start module1 */

..................................

ym2Sum = y2Sum - nyBar2Sum + n1yTilde2Sum;

   ympmSum = ypSum - nypBarSum + n1ypTildeSum;

      singular=det(ym2Sum);

      if singular = 0 then do;

       CALL SYMPUT("IsError", "1");
       print singular;
       abort;

     end;

   finish;                                   /* finish module1 */

start module2; /* start module2 */
inv_ym2Sum = inv(ym2Sum);

..................

finish; /* finish module2 */

run module1;

   run module2;

quit;

%put IsError is &IsError;

%if &IsError = 1 %then %do;
    %put Matrix is singular.;
    %goto exit;
%end;

Rick_SAS · Posted 07-15-2011 01:49 PM

1) Yes, the overflow means that the determinant is so large (in magnitude) that it cannot be stored as a double precision number. This happens often for large matrices. Your matrix dimensions are 145x145. If you take a matrix that large and fill it with uniform random numbers between 0 and 10, the determinant might be 1E200. For example: d = det(10*uniform( j(145) ));

2) See Chap 18 of the IML doc, and especially http://support.sas.com/documentation/cdl/en/imlug/63541/HTML/default/viewer.htm#imlug_genstmts_sect0...

richard_hu2003 · Posted 07-15-2011 02:58 PM

Thanks, Rick. I use determinant = 0 to detect singularity. Is there a better way to check singularity before inverse so that I can avoid overflow error of DET?

Rick_SAS · Posted 07-15-2011 04:04 PM

DET is not cheap. If the point of your computation is to solve a linear system, why not use SOLVE or INV? You're planning to abort if the matrix is nonsingular, anyway, so just try to solve the system.

For a difference between SOLVE and INV, see p. 50-51 of Statistical Programming with SAS/IML Software, which is available as a free download at http://support.sas.com/publishing/authors/wicklin.html (Also, Ch 15 discusses performance considerations between the two functions.)

Your code, as written, probably won't detect most matrices that are numerically singluar. For example, if DET returns 1E-100, your code won't set the error condition.

Rick Wicklin

SAS/IML Blog: http://blogs.sas.com/iml

Hailemichael · Posted 12-17-2011 05:03 AM

Dear Rick,

Thanks for contributing such a wonderful book on SAS IML.

I also bout the book and expecting it within coming days.

I need your support to sove the following errors from SAS IML Studio.

=========================================================================

»ERROR: NLPNRR call: Module GRAD must compute a 6 element row vector.

»ERROR: Execution error as noted previously. (rc=100)

»

» operation : NLPNRR at line 212 column 1

» operands : *LIT1085, beta0, optn, , , , , *LIT1086, *LIT1087

»

»*LIT1085 1 row 1 col (character, size 2)

»

» LL

»

»beta0 6 rows 1 col (numeric)

»

» 0

»

»optn 1 row 2 cols (numeric)

»

» 0 2

»

»*LIT1086 1 row 1 col (character, size 4)

»

» GRAD

»

»*LIT1087 1 row 1 col (character, size 4)

»

» HESS

»

» statement : CALL at line 212 column 1

==============================================================================

Please find below my SAS IML programs:

==============================================================================

SUBMIT;

%macro IPC(dat=_last_, y=, x=, beta=0);

PROC IML;

RESET noname;

use &dat;

read all var {&x} into X;

read all var {&y} into Y;

/* Design matrix with intercept */

X = j(nrow(X),1,1)||X;

P = NCOL(X);

START FMat(Y);

N1 = NROW(Y);

G1 = NCOL(UNIQUE(Y));

F = j(N1, G1, 0);

DO i = 1 TO N1;

g1 = Y;

F[i, g1] = 1;

END;

RETURN(F);

FINISH FMat;

/* Caling transformaion module */

F = FMat(Y);

N = NROW(F);

G = NCOL(F);

D = G - 1;

gamma = {0 0 0, 1 0 0, 0 1 0, 0 0 1};

* Define log likelihood function;

START LL(beta) global(F, X, gamma, N, P, G, D);

loglik = 0;

eta = j(N, D, 0);

/* Ideal points */

DO i=1 TO N;

DO j=1 TO D;

tmp = 0;

DO k=1 TO P;

tmpjk = X[i, k] * beta[(j-1)*P + k];

tmp = tmp + tmpjk;

END;

eta[i, j] = tmp;

END;

/* Unfolding/Eucliean Distance */

delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1)

+ j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma)))

- 2 * eta * t(gamma);

/* Category probability */

pi = j(N, G, 0);

DO i = 1 to N;

DO j=1 TO G;

pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]);

END;

DO i = 1 TO N;

DO j = 1 TO G;

ijloglik = F[i, j] * log(pi[i, j]);

loglik = loglik + ijloglik;

END;

RETURN(loglik);

FINISH LL;

* Define gradient;

START GRAD(beta) global(F, X, gamma, N, P, G, D);

dl = j(P*D, 1, 0);

eta = j(N, D, 0);

/* Ideal points */

DO i=1 TO N;

DO j=1 TO D;

tmp = 0;

DO k=1 TO P;

tmpjk = X[i, k] * beta[(j-1)*P + k];

tmp = tmp + tmpjk;

END;

eta[i, j] = tmp;

END;

/* Unfolding/Eucliean Distance */

delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1)

+ j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma)))

- 2 * eta * t(gamma);

/* Category probability */

pi = j(N, G, 0);

DO i = 1 to N;

DO j=1 TO G;

pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]);

END;

DO j = 1 to D;

DO k = 1 TO P;

dl[(j-1)*P + k, 1] = (F[, j] - pi[, j])` * X[, k];

END;

RETURN(dl);

FINISH GRAD;

* Define Hessian matrix;

START HESS(beta) global(F, X, gamma, N, P, G, D);

d2l = j(P*D, P*D, 0);

eta = j(N, D, 0);

/* Ideal points */

DO i=1 TO N;

DO j=1 TO D;

tmp = 0;

DO k=1 TO P;

tmpjk = X[i, k] * beta[(j-1)*P + k];

tmp = tmp + tmpjk;

END;

eta[i, j] = tmp;

END;

/* Unfolding/Eucliean Distance */

delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1)

+ j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma)))

- 2 * eta * t(gamma);

/* Category probability */

pi = j(N, G, 0);

DO i = 1 to N;

DO j=1 TO G;

pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]);

END;

jj = 1; /* jj maintains a running counter of the current row of beta_tmp and xtwx_tmp */

DO j = 1 TO D;

w1 = pi[ ,j] * (1 - pi[ ,j])`;

DO k = 1 TO P;

kk = jj - 1; /* kk maintains a running counter of the current column index xtwx_tmp */

/* second derivative when jprime=j */

DO kprime = k TO P;

kk = kk + 1;

d2tempj = d2l[jj, kk];

d2l[jj, kk] = d2tempj + (X[ ,k])` * w1 * X[ ,kprime];

d2l[kk,jj] = d2l[jj,kk];

END;

/* second derivative when jprime^=j */

DO jprime = j+1 TO D;

w2 = - pi[ ,j] * (pi[ ,jprime])`;

DO kprime = 1 TO P;

kk = kk + 1;

d2tempdprime = d2l[jj, kk];

d2l[jj, kk] = d2tempdprime + (X[ ,k])` * w2 * X[ ,kprime];

d2l[kk,jj] = d2l[jj,kk];

END;

jj = jj + 1;

END;

RETURN(d2l);

FINISH HESS;

optn = {0 2};

beta0 = β

if beta0=0 then beta0 = j(P*D, 1, 0);

else print 'User-supplied starting values for beta = ' beta0;

CALL NLPNRR(rc, beta, 'LL', beta0, optn) grd='GRAD' hes='HESS';

CALL NLPFDD(f, g, h, 'LL', beta);

var = inv(-h);

sd = sqrt(vecdiag(var));

print 'Hessian = ' h ' Covariance matrix = ' var 'Standard errors = ' sd;

run;

%mend;

%IPC(dat=nesda2, y=resp, x=sex);

ENDSUBMIT;

==============================================================

Thnk you very much for your help.

Best,

Haile

NL

Rick_SAS · Posted 12-19-2011 09:16 AM

1) For debugging, I suggest that you get rid of the macro and the SUBMIT/ENDSUBMIT block. After you have a program that works, you can wrap it in a macro if you want. For your program, use %LET stmts to define the data-specific parameters:

%let dat=nesda2; %let y=resp; %let x=sex; %let beta = 0;

2) When your program works, it will be very slow because you are assigning and manipulating matrices element-by-element. Use matrix and vector computations to speed up the computations. For example, the double loop under the comment "category probability" can be replaced by the single expression

pi = exp(-0.5 * delta) / exp(-0.5 * delta[, +]);

3) To interpret the ERROR msg, see p. 31 of my book. In your example, the error message says "Module GRAD must compute a 6 element row vector." Consequently, you should change the GRAD module to return a row vector instead of a column vector.

For an example of using SAS/IML for log-likelihood estimation, see

http://blogs.sas.com/content/iml/2011/10/12/maximum-likelihood-estimation-in-sasiml/

For two tips on defining the GRAD and HESS functions, see

http://blogs.sas.com/content/iml/2011/10/14/hints-for-derivatives/

Hailemichael · Posted 12-19-2011 10:31 AM

Hi Rick,

Thank you very much for your immediate response.

It seems now working. Your additional comments are really appreciated.

All Best,

Halie