Skip to content
Permalink
master
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
/*********************************************************************************/
/* This is the Gill/Murray cholesky routine. Reference: */
/* */
/* Gill, Jeff and Gary King. ``What to do When Your Hessian is Not Invertible: */
/* Alternatives to Model Respecification in Nonlinear Estimation,'' Sociological */
/* Methods and Research, Vol. 32, No. 1 (2004): Pp. 54--87. */
/* */
/* Abstract */
/* */
/* What should a researcher do when statistical analysis software terminates */
/* before completion with a message that the Hessian is not invertable? The */
/* standard textbook advice is to respecify the model, but this is another way */
/* of saying that the researcher should change the question being asked. */
/* Obviously, however, computer programs should not be in the business of */
/* deciding what questions are worthy of study. Although noninvertable */
/* Hessians are sometimes signals of poorly posed questions, nonsensical */
/* models, or inappropriate estimators, they also frequently occur when */
/* information about the quantities of interest exists in the data, through */
/* the likelihood function. We explain the problem in some detail and lay out */
/* two preliminary proposals for ways of dealing with noninvertable Hessians */
/* without changing the question asked. */
/* */
/* Also available is the software to implement the procedure described in this */
/* paper in R format. */
/* */
/* This procedure produces: */
/* */
/* y = chol(A+E), where E is a diagonal matrix with each element as small */
/* as possible, and A and E are the same size. E diagonal values are */
/* constrained by iteravely updated Gerschgorin bounds. */
/* */
/* Robert B. Schnabel and Elizabeth Eskow. 1990. "A New Modified Cholesky */
/* Factorization," SIAM Journal of Scientific Statistical Computating, */
/* 11, 6: 1136-58. */
/*********************************************************************************/
proc sechol(A);
local i,j,k,m,n,gamm,tau,delta,deltaprev,sum1,sum2,eigvals,dlist,dmax,
P,Ptemp,Pprod,L,norm_A,normj,g,gmax;
n = rows(A);
m = cols(A);
L = zeros(n,n);
deltaprev=0;
gamm = maxc(abs(diag(A)));
tau = __macheps^(1/3);
if minc(eig(A)') > 0;
/* print("not pivoting"); */
tau = -1000000;
endif;
if m ne n;
print("sechol: input matrix must be square");
retp(A);
endif;
norm_A = maxc(sumc(abs(A)));
gamm = maxc(abs(diag(A)));
delta = maxc(maxc(__macheps*norm_A~__macheps));
Pprod = eye(n);
if n > 2;
for k (1,(n-2),1);
trap 1,1;
if minc((diag(A[(k+1):n,(k+1):n])' - A[k,(k+1):n]^2/A[k,k])') < tau*gamm
and minc(eig(A[(k+1):n,(k+1):n])) < 0;
dmax = maxindc(diag(A[k:n,k:n]));
if (A[(k+dmax-1),(k+dmax-1)] > A[k,k]);
/* print("pivot using dmax:"); print(dmax); */
P = eye(n);
Ptemp = P[k,.]; P[k,.] = P[(k+dmax-1),.]; P[(k+dmax-1),.] = Ptemp;
A = P*A*P;
L = P*L*P;
Pprod = P*Pprod;
endif;
g = zeros(n-(k-1),1);
for i ((k), (n), 1);
if i == 1;
sum1 = 0;
else;
sum1 = sumc(abs(A[i,k:(i-1)])');
endif;
if i == n;
sum2 = 0;
else;
sum2 = sumc(abs(A[(i+1):n,i]));
endif;
g[i-(k-1)] = A[i,i] - sum1 - sum2;
endfor;
gmax = maxindc(g);
if gmax /= k;
/* print("gerschgorin pivot on cycle:"); print(k); */
P = eye(n);
Ptemp = P[k,.]; P[k,.] = P[(k+dmax-1),.]; P[(k+dmax-1),.] = Ptemp;
A = P*A*P;
L = P*L*P;
Pprod = P*Pprod;
endif;
normj = sumc(abs(A[(k+1):n,k]));
delta = maxc((0~deltaprev~-A[k,k]+(maxc(normj~tau*gamm))')');
if delta > 0;
A[k,k] = A[k,k] + delta;
deltaprev = delta;
endif;
endif;
A[k,k] = sqrt(A[k,k]);
L[k,k] = A[k,k];
for i ((k+1), (n), 1);
if L[k,k] > __macheps; A[i,k] = A[i,k]/L[k,k]; endif;
L[i,k] = A[i,k];
A[i,(k+1):i] = A[i,(k+1):i] - L[i,k]*L[(k+1):i,k]';
if A[i,i] < 0; A[i,i] = 0; endif;
endfor;
endfor;
endif;
A[(n-1),n] = A[n,(n-1)];
eigvals = eig(A[(n-1):n,(n-1):n]);
dlist = ( (0|deltaprev|-minc(eigvals)+tau*maxc( (1/(1-tau))*(maxc(eigvals)-minc(eigvals))|gamm)) );
if dlist[1] > dlist[2];
delta = dlist[1];
else;
delta = dlist[2];
endif;
if delta < dlist[3];
delta = dlist[3];
endif;
if delta > 0;
A[(n-1),(n-1)] = A[(n-1),(n-1)] + delta;
A[n,n] = A[n,n] + delta;
deltaprev = delta;
endif;
A[(n-1),(n-1)] = sqrt(A[(n-1),(n-1)]);
L[(n-1),(n-1)] = A[(n-1),(n-1)];
A[n,(n-1)] = A[n,(n-1)]/L[(n-1),(n-1)];
L[n,(n-1)] = A[n,(n-1)];
A[n,n] = sqrt(A[n,n] - L[n,(n-1)]^2);
L[n,n] = A[n,n];
retp(Pprod'*L'*Pprod');
endp;