linear_operators.iterative.dli

1 """ 2 Implements Double loop inference algorithms. 3 4 Reference 5 --------- 6 7 Bayesian Inference and Optimal Design for the Sparse Linear Model, 8 Matthias W. Seeger 9 10 http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.165.8284&rep=rep1&type=pdf 11 12 """ 13 from copy import copy 14 import numpy as np 15 import numexpr as ne 16 from utils import eigendecomposition 17 from algorithms import Algorithm, default_callback, StopCondition 18 from ..operators import diag, identity 19 from optimize import * 20 from norms import norm2, dnorm2, Norm2 21 22 default_stop = StopCondition(maxiter=5) 23 24 # lanczos algorithm 25

26 -class LanczosAlgorithm(Algorithm):

27 """ 28 Tridiagonalization Lanczos step and eigendecomposition at exit. 29 30 http://en.wikipedia.org/wiki/Lanczos_algorithm 31 """

32 - def __init__(self, A, **kwargs):

33 self.A = A 34 self.n = self.A.shape[0] 35 self.kwargs = kwargs 36 # extract appropriate kwargs for stop condition 37 maxiter = getattr(kwargs, "maxiter", 300) 38 tol = getattr(kwargs, "tol", None) 39 gtol = getattr(kwargs, "gtol", None) 40 self.stop_condition = StopCondition(maxiter=maxiter, tol=tol, gtol=gtol) 41 # maxiter default to matrix size if not given. 42 self.maxiter = getattr(self.stop_condition, "maxiter", self.n) 43 Algorithm.__init__(self)

44

45 - def initialize(self):

46 Algorithm.initialize(self) 47 # to store results 48 # tridiagonal matrix coefficients 49 self.alpha = np.zeros(self.maxiter + 1) 50 self.beta = np.zeros(self.maxiter) 51 self.vectors = np.zeros((self.maxiter + 1, self.n)) 52 self.w = np.zeros(self.n) 53 # starting point 54 self.vectors[0] = np.random.randn(self.n) 55 self.vectors[0] /= np.sqrt(norm2(self.vectors[0]))

56

57 - def iterate(self):

58 self.orthogonalization() 59 self.update_alpha() 60 self.update_w() 61 self.update_beta() 62 self.update_vectors() 63 self.iter_ += 1

64

65 - def orthogonalization(self):

66 self.w = self.A * self.vectors[self.iter_] 67 if self.iter_ > 0: 68 self.w -= self.beta[self.iter_ - 1] * self.vectors[self.iter_ - 1]

69

70 - def update_alpha(self):

71 self.alpha[self.iter_] = np.dot(self.w, self.vectors[self.iter_])

72

73 - def update_w(self):

74 self.w -= self.alpha[self.iter_] * self.vectors[self.iter_]

75

76 - def update_beta(self):

77 self.beta[self.iter_] = np.sqrt(norm2(self.w))

78

79 - def update_vectors(self):

80 self.vectors[self.iter_ + 1] = self.w / self.beta[self.iter_]

81

82 - def at_exit(self):

83 """ 84 Convert alpha and beta to a TridiagonalOperator and perform 85 eigendecomposition. 86 """ 87 from ..operators import SymmetricTridiagonal 88 from ..operators import EigendecompositionOperator 89 self.T = SymmetricTridiagonal(2 * (self.maxiter + 1,), 90 self.alpha, self.beta) 91 # use band matrix eigendecomposition as tridiagonal lapack 92 # routine to available 93 self.B = self.T.toband() 94 #select_range = [self.n - 1 - self.maxiter, self.n - 1] 95 #self.E = self.B.eigen(select="i", select_range=select_range) 96 self.E = self.B.eigen() 97 # multiply T eigenvectors with lanczos vectors 98 w = self.E.eigenvalues 99 v = np.zeros((self.n, self.maxiter + 1)) 100 for i in xrange(self.E.eigenvectors.shape[1]): 101 v[:, i] = np.dot(self.vectors.T, self.E.eigenvectors[:, i]) 102 # remove the last eigenpair with negative eigenvalue XXX 103 w = w[1:] 104 v = v[:, 1:] 105 self.current_solution = EigendecompositionOperator(v=v, w=w)

106

107 -class Criterion(object):

108 - def __init__(self, algo):

109 self.algo = algo 110 self.n_variables = self.algo.model.shape[1] 111 # likelihood norm 112 self.norm = Norm2(C=algo.noise_covariance) 113 # storing 114 self.last_u = None 115 self.Xu = None 116 self.Bu = None

117 - def islast(self, u):

118 return np.all(u == self.last_u)

119 - def load_last(self):

120 return self.Xu, self.Bu

121 - def get_projections(self, u):

122 if self.islast(u): 123 return self.load_last() 124 else: 125 self.last_u = copy(u) 126 X = self.algo.model 127 B = self.algo.prior 128 self.Xu = X * u 129 self.Bu = B * u 130 return self.Xu, self.Bu

131 - def likelihood(self, u):

132 sigma = self.algo.sigma 133 X = self.algo.model 134 y = self.algo.data 135 Xu, Bu = self.get_projections(u) 136 return sigma ** (-2) * self.norm(Xu - y)

137 - def dlike(self, u):

138 sigma = self.algo.sigma 139 X = self.algo.model 140 y = self.algo.data 141 Xu, Bu = self.get_projections(u) 142 return sigma ** (-2) * X.T * self.norm.diff(Xu - y)

143 - def d2like(self, u):

144 sigma = self.algo.sigma 145 X = self.algo.model 146 I = identity(2 * (X.shape[0],)) 147 N = getattr(self.algo, "noise_covariance", I) 148 if N is None: 149 N = I 150 return sigma ** (-2) * X.T * N * X

151 - def d2lik_p(self, u, p):

152 return self.d2like(u) * p

153 - def penalization(self, u):

154 sigma = self.algo.sigma 155 B = self.algo.prior 156 t = self.algo.tau 157 z = self.algo.z 158 #e = t * np.sqrt(z + (np.abs(B * u) / sigma) ** 2) 159 Xu, Bu = self.get_projections(u) 160 e = ne.evaluate("2 * t * sqrt(z + (abs(Bu) / sigma) ** 2)") 161 return e.sum()

162 - def dpen(self, u):

163 sigma = self.algo.sigma 164 B = self.algo.prior 165 t = self.algo.tau 166 z = self.algo.z 167 Xu, Bu = self.get_projections(u) 168 #e = 2 * (t * Bu) / np.sqrt(z + (Bu / sigma) ** 2) 169 e = ne.evaluate("2 * (t * Bu) / sqrt(z + (Bu / sigma) ** 2)") 170 return (B.T * e) / (sigma ** 2)

171 - def d2pen(self, u):

172 sigma = self.algo.sigma 173 B = self.algo.prior 174 t = self.algo.tau 175 z = self.algo.z 176 Xu, Bu = self.get_projections(u) 177 rho = ne.evaluate("(t * z) / ((z + (Bu / sigma) ** 2) ** (1.5) * sigma ** 2)") 178 return B.T * diag(rho) * B

179 - def d2pen_p(self, u, p):

180 return self.d2pen(u) * p

181 - def __call__(self, u):

182 return self.likelihood(u) + self.penalization(u)

183 - def gradient(self, u):

184 return self.dlike(u) + self.dpen(u)

185 - def hessian(self, u):

186 return self.d2like(u) + self.d2pen(u)

187 - def hessian_p(self, u, p):

188 return self.hessian(u) * p

189

190 -class DoubleLoopAlgorithm(Algorithm):

191 """ 192 A class implementing the double loop algorithm. 193 194 Parameters 195 ---------- 196 197 model : LinearOperator 198 Linear model linking data and unknowns. 199 data : ndarray 200 Data. 201 prior : LinearOperator 202 Prior. 203 tau : ndarray (optional) 204 Parameters of the Laplace potential on priors coefficients. 205 sigma : float (optional) 206 Likelihood standard deviation. 207 lanczos : dict 208 Keyword arguments of the Lanczos decomposition. 209 fmin_args : dict 210 Keyword arguments of the function minimization. 211 """

212 - def __init__(self, model, data, prior, noise_covariance=None, 213 tau=None, sigma=1., optimizer=FminSLSQP, 214 lanczos={}, fmin_args={}, 215 callback=default_callback, 216 stop_condition=default_stop, 217 ):

218 from ..interface import aslinearoperator 219 220 self.model = aslinearoperator(model) 221 self.data = data 222 self.prior = aslinearoperator(prior) 223 if noise_covariance is not None: 224 noise_covariance = aslinearoperator(noise_covariance) 225 self.noise_covariance = noise_covariance 226 if tau is None: 227 self.tau = np.ones(prior.shape[0]) 228 else: 229 if tau.size != prior.shape[0]: 230 raise ValueError("Incorrect shape for tau.") 231 self.tau = tau 232 self.sigma = sigma 233 self.optimizer = optimizer 234 self.lanczos = lanczos 235 self.fmin_args = fmin_args 236 # 237 self.callback = callback 238 self.stop_condition = stop_condition 239 # to store internal variables 240 self.z = None 241 self.gamma = None 242 self.inv_gamma = None 243 self.g_star = None 244 self.current_solution = None 245 self.last_solution = None 246 self.inv_cov = None 247 self.inv_cov_approx = None 248 self.criterion = None

249 - def initialize(self):

250 """ 251 Set parameters to initial values. 252 """ 253 self.z = 0.05 * np.ones(self.model.shape[1]) 254 self.g_star = 0. 255 self.current_solution = np.zeros(self.model.shape[1]) 256 self.iter_ = 0 257 self.gamma = np.ones(self.prior.shape[0]) 258 self.update_inv_gamma()

259 - def iterate(self):

260 print "outer" 261 self.outer() 262 print "inner" 263 self.inner() 264 return Algorithm.iterate(self)

265 # outer loop

266 - def outer(self):

267 """ 268 Outer loop : Lanczos approximation. 269 """ 270 self.update_inv_cov() 271 self.update_inv_cov_approx() 272 self.update_z() 273 self.update_g_star()

274 - def update_inv_cov(self):

275 D = diag(self.gamma ** (-1), dtype=self.prior.dtypeout) 276 X = self.model 277 B = self.prior 278 N = self.noise_covariance 279 if N is None: 280 self.inv_cov = X.T * X + B.T * D * B 281 else: 282 self.inv_cov = X.T * N * X + B.T * D * B

283 - def update_inv_cov_approx(self):

284 self.lanczos_algorithm = LanczosAlgorithm(self.inv_cov, **self.lanczos) 285 self.inv_cov_approx = self.lanczos_algorithm()

286 - def update_z(self):

287 # get eigenvalues, eigenvectors 288 e = self.inv_cov_approx.eigenvalues 289 v = self.inv_cov_approx.eigenvectors 290 B = self.prior 291 self.z = sum([ei * (B * vi) ** 2 for ei, vi in zip(e, v.T)])

292 - def update_g_star(self):

293 self.g_star = np.dot(self.z.T, self.inv_gamma) 294 self.g_star -= self.inv_cov_approx.logdet()

295 # inner loop

296 - def inner(self):

297 """ 298 Inner loop : Penalized minimization. 299 """ 300 self.update_current_solution() 301 self.update_gamma() 302 self.update_inv_gamma()

303 - def update_current_solution(self):

304 self.inner_criterion = Criterion(self) 305 self.last_solution = copy(self.current_solution) 306 self.inner_algo = self.optimizer(self.inner_criterion, 307 self.last_solution, 308 **self.fmin_args) 309 self.current_solution = self.inner_algo()

310 - def update_gamma(self):

311 s = np.abs(self.prior * self.current_solution) 312 sn2 = (s / self.sigma) ** 2 313 self.gamma = np.sqrt(self.z + sn2) / self.tau

314 - def update_inv_gamma(self):

315 self.inv_gamma = self.gamma ** (-1)

316

Source Code for Module linear_operators.iterative.dli