/
sechol.g.txt
138 lines (134 loc) · 5.22 KB
/
sechol.g.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
/*********************************************************************************/
/* 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;