Loading...
Searching...
No Matches
glassofast-real.f90
Go to the documentation of this file.
! Copyright 2012 Matyas M. Sustik <msustik@gmail.com>
! 2012 B. Calderhead
! 2013 L. Kajan <lkajan@debian.org>
!
! License:
! URL: http://www.cs.utexas.edu/users/sustik/glassofast/
!
! GLASSOFAST software
! The program
! GLASSOFAST is a fast implementation of GLASSO (latest release 1.0). It solves
! the l1 regularized Gaussian maximum likelihood estimation of the inverse of a
! covariance matrix.
! Download
! We implemented the algorithm in Fortran and we provide a MEX package released
! under the GNU General Public License version 3 or later (GPLv3).
!
! Laszlo Kajan <lkajan@rostlab.org> 2012 - patched glassofast to work with
! either single or double precision. Patch was sent upstream.
!
subroutine glassofast(n,S,L,thr,maxIt,msg,warm,X,W,Wd,WXj,info,brk)
!
! .. Scalar Arguments ..
integer n, warm, msg, maxIt
! ..
! .. Array Arguments ..
real S(n,n), L(n,n), X(n,n), W(n,n), Wd(n), WXj(n)
! ..
!
! Purpose
! =======
!
! This subroutine computes the L1 regularized covariance matrix estimate
! using the algorithm described in the paper:
! J. Friedman, T. Hastie, R. Tibshirani:
! Sparse inverse covariance estimation with the graphical lasso
! Biostatistics, 9(3):432-441, July 2008.
! This implementation is documented in the technical report:
! M. A. Sustik B. Calderhead:
! GLASSOFAST: An efficient GLASSO implementation
! Technical Report TR-12-29, University of Texas at Austin
!
! Arguments
! =========
!
! n (input) integer
! The dimension of the input matrix s
!
! S (input) double precision array, dimension n x n
! The empirical covariance matrix
!
! L (input) double precision array, dimension n x n
! Regularization matrix (symmetric)
!
! thr (input) double precision
! Convergence threshold
!
! maxIt (input/output) integer
! Maximum number of whole matrix sweeps on input, actual number
! of sweeps on output
!
! msg (input) integer
! Controls amount of messages printed
!
! warm (input) integer
! flag indicating cold versus warm start, see also X, W
!
! X (input/output) double precision array, dimension n x n
! Inverse covariance matrix estimate
!
! W (input/output) double precision array, dimension n x n
! Covariance matrix estimate
!
! Wd (input) double precision array, dimension n
! Internally used workspace
!
! WXj (input) double precision array, dimension n
! Internally used workspace
!
! info (output) integer
! Indicates errors
! 256 = brk was activated
!
! brk (input) integer
! dlx loop control: return when brk .ne. 0
!
! Wd, WDXj are not allocatable arrays, because gfortran/gdb has trouble
! displaying those.
integer iter
real EPS
parameter(eps = epsilon(eps))
info = 0
shr = sum(abs(s))
do i = 1,n
shr = shr - abs(s(i, i))
enddo
if (shr .eq. 0.0) then
! S is diagonal.
w = 0.0
x = 0.0
do i = 1,n
w(i,i) = w(i,i) + l(i,i)
enddo
x = 0.0
do i = 1,n
x(i,i) = 1.0/max(w(i,i),eps)
enddo
return
endif
shr = thr*shr/(n-1)
thrlasso = shr/n
if (thrlasso .lt. 2*eps) then
thrlasso = 2*eps
end if
if (warm .eq. 0) then
w = s
x = 0.0
else
do i = 1,n
x(1:n,i) = -x(1:n,i)/x(i,i)
x(i,i) = 0
end do
end if
do i = 1,n
wd(i) = s(i,i) + l(i,i)
w(i,i) = wd(i)
end do
do iter = 1,maxit
if (msg .ne. 0) write(0,"(A,I4)",advance='no') "␛[1K␛[1G iteration =", iter
dw = 0.0
do j = 1,n
wxj(1:n) = 0.0
! We exploit sparsity of X when computing column j of W*X*D:
do i = 1,n
if (x(i,j) .ne. 0.0) then
wxj = wxj + w(:,i)*x(i,j)
endif
enddo
do
dlx = 0.0
do i = 1,n
if (i .ne. j) then
a = s(i,j) - wxj(i) + wd(i)*x(i,j)
b = abs(a) - l(i,j)
if (b .gt. 0.0) then
c = sign(b, a)/wd(i)
else
c = 0.0
endif
delta = c - x(i,j)
if (delta .ne. 0.0) then
x(i,j) = c
wxj(1:n) = wxj(1:n) + w(:,i)*delta
dlx = max(dlx, abs(delta))
endif
endif
enddo
if (dlx .lt. thrlasso) then
exit
endif
if (brk .ne. 0) then
info = 256; return
endif
enddo
wxj(j) = wd(j)
dw = max(dw, sum(abs(wxj(1:n) - w(:,j))))
w(:,j) = wxj(1:n)
w(j,:) = wxj(1:n)
enddo
! write(0,*) " dw =", dw
if (dw .le. shr) then
exit
endif
enddo
if (msg .ne. 0) write(0,*)
do i = 1,n
tmp = 1/(wd(i) - sum(x(:,i)*w(:,i)))
x(1:n,i) = -tmp*x(1:n,i)
x(i,i) = tmp
enddo
do i = 1,n-1
x(i+1:n,i) = (x(i+1:n,i) + x(i,i+1:n))/2;
x(i,i+1:n) = x(i+1:n,i)
enddo
maxit = iter
return
end subroutine glassofast
glassofast
subroutine glassofast(n, s, l, thr, maxit, msg, warm, x, w, wd, wxj, info, brk
Definition glassofast-real.f90:21