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
|
# The covariance selection example at the end of chapter 7 (Sparse linear
# equations).
from cvxopt import matrix, spmatrix, log, mul, blas, lapack, amd, cholmod
from pickle import load
def covsel(Y):
"""
Returns the solution of
minimize -log det K + tr(KY)
subject to K_ij = 0 if (i,j) not in zip(I, J).
Y is a symmetric sparse matrix with nonzero diagonal elements.
I = Y.I, J = Y.J.
"""
cholmod.options['supernodal'] = 2
I, J = Y.I, Y.J
n, m = Y.size[0], len(I)
# non-zero positions for one-argument indexing
N = I + J*n
# position of diagonal elements
D = [ k for k in range(m) if I[k]==J[k] ]
# starting point: symmetric identity with nonzero pattern I,J
K = spmatrix(0.0, I, J)
K[::n+1] = 1.0
# Kn is used in the line search
Kn = spmatrix(0.0, I, J)
# symbolic factorization of K
F = cholmod.symbolic(K)
# Kinv will be the inverse of K
Kinv = matrix(0.0, (n,n))
for iters in range(100):
# numeric factorization of K
cholmod.numeric(K, F)
d = cholmod.diag(F)
# compute Kinv by solving K*X = I
Kinv[:] = 0.0
Kinv[::n+1] = 1.0
cholmod.solve(F, Kinv)
# solve Newton system
grad = 2 * (Y.V - Kinv[N])
hess = 2 * ( mul(Kinv[I,J], Kinv[J,I]) +
mul(Kinv[I,I], Kinv[J,J]) )
v = -grad
lapack.posv(hess,v)
# stopping criterion
sqntdecr = -blas.dot(grad,v)
print("Newton decrement squared:%- 7.5e" %sqntdecr)
if (sqntdecr < 1e-12):
print("number of iterations: %d" %(iters+1))
break
# line search
dx = +v
dx[D] *= 2
f = -2.0*sum(log(d)) # f = -log det K
s = 1
for lsiter in range(50):
Kn.V = K.V + s*dx
try:
cholmod.numeric(Kn, F)
except ArithmeticError:
s *= 0.5
else:
d = cholmod.diag(F)
fn = -2.0 * sum(log(d)) + 2*s*blas.dot(v,Y.V)
if (fn < f - 0.01*s*sqntdecr): break
else: s *= 0.5
K.V = Kn.V
return K
Y = load(open("covsel.bin","rb"))
print("%d rows/columns, %d nonzeros\n" %(Y.size[0], len(Y)))
covsel(Y)
|
