libfreecontact  1.0.21
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 
89 integer iter
90 real EPS
91 parameter(eps = epsilon(eps))
92 info = 0
93 shr = sum(abs(s))
94 do i = 1,n
95  shr = shr - abs(s(i, i))
96 enddo
97 if (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
109 endif
110 shr = thr*shr/(n-1)
111 thrlasso = shr/n
112 if (thrlasso .lt. 2*eps) then
113  thrlasso = 2*eps
114 end if
115 if (warm .eq. 0) then
116  w = s
117  x = 0.0
118 else
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
123 end if
124 do i = 1,n
125  wd(i) = s(i,i) + l(i,i)
126  w(i,i) = wd(i)
127 end do
128 do iter = 1,maxit
129  if (msg .ne. 0) write(0,"(A,I4)",advance='no') " 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
174 enddo
175 if (msg .ne. 0) write(0,*)
176 do 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
180 enddo
181 do 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)
184 enddo
185 maxit = iter
186 return
187 end subroutine glassofast
glassofast
subroutine glassofast(n, S, L, thr, maxIt, msg, warm, X, W, Wd, WXj, info, brk
Definition: glassofast-real.f90:21