libfreecontact 1.0.21
Loading...
Searching...
No Matches
glassofast-real.f90
Go to the documentation of this file.
1! Copyright 2012 Matyas M. Sustik <msustik@gmail.com>
2! 2012 B. Calderhead
3! 2013 L. Kajan <lkajan@debian.org>
4!
5! License:
6! URL: http://www.cs.utexas.edu/users/sustik/glassofast/
7!
8! GLASSOFAST software
9! The program
10! GLASSOFAST is a fast implementation of GLASSO (latest release 1.0). It solves
11! the l1 regularized Gaussian maximum likelihood estimation of the inverse of a
12! covariance matrix.
13! Download
14! We implemented the algorithm in Fortran and we provide a MEX package released
15! under the GNU General Public License version 3 or later (GPLv3).
16!
17! Laszlo Kajan <lkajan@rostlab.org> 2012 - patched glassofast to work with
18! either single or double precision. Patch was sent upstream.
19!
20 subroutine glassofast(n,S,L,thr,maxIt,msg,warm,X,W,Wd,WXj,info,brk)
21!
22! .. Scalar Arguments ..
23 integer n, warm, msg, maxIt
24! ..
25! .. Array Arguments ..
26 real S(n,n), L(n,n), X(n,n), W(n,n), Wd(n), WXj(n)
27! ..
28!
29! Purpose
30! =======
31!
32! This subroutine computes the L1 regularized covariance matrix estimate
33! using the algorithm described in the paper:
34! J. Friedman, T. Hastie, R. Tibshirani:
35! Sparse inverse covariance estimation with the graphical lasso
36! Biostatistics, 9(3):432-441, July 2008.
37! This implementation is documented in the technical report:
38! M. A. Sustik B. Calderhead:
39! GLASSOFAST: An efficient GLASSO implementation
40! Technical Report TR-12-29, University of Texas at Austin
41!
42! Arguments
43! =========
44!
45! n (input) integer
46! The dimension of the input matrix s
47!
48! S (input) double precision array, dimension n x n
49! The empirical covariance matrix
50!
51! L (input) double precision array, dimension n x n
52! Regularization matrix (symmetric)
53!
54! thr (input) double precision
55! Convergence threshold
56!
57! maxIt (input/output) integer
58! Maximum number of whole matrix sweeps on input, actual number
59! of sweeps on output
60!
61! msg (input) integer
62! Controls amount of messages printed
63!
64! warm (input) integer
65! flag indicating cold versus warm start, see also X, W
66!
67! X (input/output) double precision array, dimension n x n
68! Inverse covariance matrix estimate
69!
70! W (input/output) double precision array, dimension n x n
71! Covariance matrix estimate
72!
73! Wd (input) double precision array, dimension n
74! Internally used workspace
75!
76! WXj (input) double precision array, dimension n
77! Internally used workspace
78!
79! info (output) integer
80! Indicates errors
81! 256 = brk was activated
82!
83! brk (input) integer
84! dlx loop control: return when brk .ne. 0
85!
86! Wd, WDXj are not allocatable arrays, because gfortran/gdb has trouble
87! displaying those.
88
89integer iter
90real EPS
91parameter(eps = epsilon(eps))
92info = 0
93shr = sum(abs(s))
94do i = 1,n
95 shr = shr - abs(s(i, i))
96enddo
97if (shr .eq. 0.0) then
98! S is diagonal.
99 w = 0.0
100 x = 0.0
101 do i = 1,n
102 w(i,i) = w(i,i) + l(i,i)
103 enddo
104 x = 0.0
105 do i = 1,n
106 x(i,i) = 1.0/max(w(i,i),eps)
107 enddo
108 return
109endif
110shr = thr*shr/(n-1)
111thrlasso = shr/n
112if (thrlasso .lt. 2*eps) then
113 thrlasso = 2*eps
114end if
115if (warm .eq. 0) then
116 w = s
117 x = 0.0
118else
119 do i = 1,n
120 x(1:n,i) = -x(1:n,i)/x(i,i)
121 x(i,i) = 0
122 end do
123end if
124do i = 1,n
125 wd(i) = s(i,i) + l(i,i)
126 w(i,i) = wd(i)
127end do
128do iter = 1,maxit
129 if (msg .ne. 0) write(0,"(A,I4)",advance='no') "␛[1K␛[1G iteration =", iter
130 dw = 0.0
131 do j = 1,n
132 wxj(1:n) = 0.0
133! We exploit sparsity of X when computing column j of W*X*D:
134 do i = 1,n
135 if (x(i,j) .ne. 0.0) then
136 wxj = wxj + w(:,i)*x(i,j)
137 endif
138 enddo
139 do
140 dlx = 0.0
141 do i = 1,n
142 if (i .ne. j) then
143 a = s(i,j) - wxj(i) + wd(i)*x(i,j)
144 b = abs(a) - l(i,j)
145 if (b .gt. 0.0) then
146 c = sign(b, a)/wd(i)
147 else
148 c = 0.0
149 endif
150 delta = c - x(i,j)
151 if (delta .ne. 0.0) then
152 x(i,j) = c
153 wxj(1:n) = wxj(1:n) + w(:,i)*delta
154 dlx = max(dlx, abs(delta))
155 endif
156 endif
157 enddo
158 if (dlx .lt. thrlasso) then
159 exit
160 endif
161 if (brk .ne. 0) then
162 info = 256; return
163 endif
164 enddo
165 wxj(j) = wd(j)
166 dw = max(dw, sum(abs(wxj(1:n) - w(:,j))))
167 w(:,j) = wxj(1:n)
168 w(j,:) = wxj(1:n)
169 enddo
170! write(0,*) " dw =", dw
171 if (dw .le. shr) then
172 exit
173 endif
174enddo
175if (msg .ne. 0) write(0,*)
176do i = 1,n
177 tmp = 1/(wd(i) - sum(x(:,i)*w(:,i)))
178 x(1:n,i) = -tmp*x(1:n,i)
179 x(i,i) = tmp
180enddo
181do i = 1,n-1
182 x(i+1:n,i) = (x(i+1:n,i) + x(i,i+1:n))/2;
183 x(i,i+1:n) = x(i+1:n,i)
184enddo
185maxit = iter
186return
187end subroutine glassofast
subroutine glassofast(n, s, l, thr, maxit, msg, warm, x, w, wd, wxj, info, brk