Timer unit: 1e-06 s Total time: 68.7961 s File: /user/cso/NoBackup/Magni/magni/cs/reconstruction/gamp/_algorithm.py Function: run at line 68 Line # Hits Time Per Hit % Time Line Contents ============================================================== 68 def run(y, A, A_asq=None): 69 """ 70 Run the GAMP reconstruction algorithm. 71 72 Parameters 73 ---------- 74 y : ndarray 75 The m x 1 measurement vector. 76 A : ndarray or magni.utils.matrices.{Matrix, MatrixCollection} 77 The m x n matrix which is the product of the measurement matrix and the 78 dictionary matrix. 79 A_asq : ndarray or magni.utils.matrices.{Matrix, MatrixCollection} or None 80 The m x n matrix which is the entrywise absolute value squared product 81 of the measurement matrix and the dictionary matrix (the default is 82 None, which implies that a sum approximation is used). 83 84 Returns 85 ------- 86 alpha : ndarray 87 The n x 1 reconstructed coefficient vector. 88 history : dict, optional 89 The dictionary of various measures tracked in the GAMP iterations. 90 91 See Also 92 -------- 93 magni.cs.reconstruction.gamp._config : Configuration options. 94 magni.cs.reconstruction.gamp.input_channel : Input channels. 95 magni.cs.reconstruction.gamp.output_channel : Output channels. 96 magni.cs.reconstruction.gamp.stop_criterion : Stop criteria. 97 98 Notes 99 ----- 100 Optionally, the algorithm may be configured to save and return the 101 iteration history along with the reconstruction result. The returned 102 history contains the following: 103 104 * alpha_bar : Mean coefficient estimates (the reconstruction coefficients). 105 * alpha_tilde : Variance coefficient estimates. 106 * MSE : solution Mean squared error (if the true solution is known). 107 * input_channel_parameters : The state of the input channel. 108 * output_channel_parameters : The state of the output channel. 109 * stop_criterion : The currently used stop criterion. 110 * stop_criterion_value : The value of the stop criterion. 111 * stop_iteration : The iteration at which the algorithm stopped. 112 * stop_reason : The reason for termination of the algorithm. 113 114 Examples 115 -------- 116 For example, recovering a vector from AWGN noisy measurements using GAMP 117 118 >>> import numpy as np 119 >>> from magni.cs.reconstruction.gamp import run, config 120 >>> np.random.seed(seed=6028) 121 >>> k, m, n = 10, 200, 400 122 >>> tau = float(k) / n 123 >>> A = 1 / np.sqrt(m) * np.random.randn(m, n) 124 >>> A_asq = np.abs(A)**2 125 >>> alpha = np.zeros((n, 1)) 126 >>> alpha[:k] = np.random.normal(scale=2, size=(k, 1)) 127 >>> np_printoptions = np.get_printoptions() 128 >>> np.set_printoptions(suppress=True, threshold=k+2) 129 >>> alpha[:k + 2] 130 array([[ 1.92709461], 131 [ 0.74378508], 132 [-3.2418159 ], 133 [-1.32277347], 134 [ 0.90118 ], 135 [-0.19157262], 136 [ 0.82855712], 137 [ 0.24817994], 138 [-1.43034777], 139 [-0.21232344], 140 [ 0. ], 141 [ 0. ]]) 142 >>> sigma = 0.15 143 >>> y = A.dot(alpha) + np.random.normal(scale=sigma, size=(A.shape[0], 1)) 144 >>> input_channel_params = {'tau': tau, 'theta_bar': 0, 'theta_tilde': 4, 145 ... 'use_em': False} 146 >>> config['input_channel_parameters'] = input_channel_params 147 >>> output_channel_params = {'sigma_sq': sigma**2, 148 ... 'noise_level_estimation': 'fixed'} 149 >>> config['output_channel_parameters'] = output_channel_params 150 >>> alpha_hat = run(y, A, A_asq) 151 >>> alpha_hat[:k + 2] 152 array([[ 1.93810961], 153 [ 0.6955502 ], 154 [-3.39759349], 155 [-1.35533562], 156 [ 1.10524227], 157 [-0.00594848], 158 [ 0.79274671], 159 [ 0.04895264], 160 [-1.08726071], 161 [-0.00142911], 162 [ 0.00022861], 163 [-0.00004272]]) 164 >>> np.sum(np.abs(alpha - alpha_hat) > sigma * 3) 165 0 166 167 or recover the same vector returning a history comparing the pr. iteration 168 solution to the true vector and printing the A_asq details 169 170 >>> config['report_A_asq_setup'] = True 171 >>> config['report_history'] = True 172 >>> config['true_solution'] = alpha 173 >>> alpha_hat, history = run(y, A, A_asq) # doctest: +NORMALIZE_WHITESPACE 174 GAMP is using the A_asq: [[ 0.024 ..., 0.002] 175 ..., 176 [ 0. ..., 0.014]] 177 The sum approximation method is: None 178 >>> alpha_hat[:k + 2] 179 array([[ 1.93810961], 180 [ 0.6955502 ], 181 [-3.39759349], 182 [-1.35533562], 183 [ 1.10524227], 184 [-0.00594848], 185 [ 0.79274671], 186 [ 0.04895264], 187 [-1.08726071], 188 [-0.00142911], 189 [ 0.00022861], 190 [-0.00004272]]) 191 >>> np.array(history['MSE']).reshape(-1, 1)[1:11] 192 array([[ 0.04562729], 193 [ 0.01328304], 194 [ 0.00112098], 195 [ 0.00074968], 196 [ 0.00080175], 197 [ 0.00076615], 198 [ 0.00077043]]) 199 200 or recover the same vector using sample variance AWGN noise level 201 estimation 202 203 >>> config['report_A_asq_setup'] = False 204 >>> config['report_history'] = False 205 >>> output_channel_params['noise_level_estimation'] = 'sample_variance' 206 >>> config['output_channel_parameters'] = output_channel_params 207 >>> alpha_hat = run(y, A, A_asq) 208 >>> alpha_hat[:k + 2] 209 array([[ 1.94820622], 210 [ 0.72162206], 211 [-3.39978431], 212 [-1.35357001], 213 [ 1.10701779], 214 [-0.00834467], 215 [ 0.79790879], 216 [ 0.08441384], 217 [-1.08946306], 218 [-0.0015894 ], 219 [ 0.00020561], 220 [-0.00003623]]) 221 >>> np.sum(np.abs(alpha - alpha_hat) > sigma * 3) 222 0 223 224 or recover the same vector using median AWGN noise level estimation 225 226 >>> output_channel_params['noise_level_estimation'] = 'median' 227 >>> config['output_channel_parameters'] = output_channel_params 228 >>> alpha_hat = run(y, A, A_asq) 229 >>> alpha_hat[:k + 2] 230 array([[ 1.93356483], 231 [ 0.65232347], 232 [-3.39440429], 233 [-1.35437724], 234 [ 1.10312573], 235 [-0.0050555 ], 236 [ 0.78743162], 237 [ 0.03616397], 238 [-1.08589927], 239 [-0.00136802], 240 [ 0.00024121], 241 [-0.00004498]]) 242 >>> np.sum(np.abs(alpha - alpha_hat) > sigma * 3) 243 0 244 245 or recover the same vector learning the AWGN noise level using expectation 246 maximization (EM) 247 248 >>> output_channel_params['noise_level_estimation'] = 'em' 249 >>> config['output_channel_parameters'] = output_channel_params 250 >>> alpha_hat = run(y, A, A_asq) 251 >>> alpha_hat[:k + 2] 252 array([[ 1.94118089], 253 [ 0.71553983], 254 [-3.40076165], 255 [-1.35662005], 256 [ 1.1099417 ], 257 [-0.00688125], 258 [ 0.79442879], 259 [ 0.06258856], 260 [-1.08792606], 261 [-0.00148811], 262 [ 0.00022266], 263 [-0.00003785]]) 264 >>> np.sum(np.abs(alpha - alpha_hat) > sigma * 3) 265 0 266 267 >>> np.set_printoptions(**np_printoptions) 268 269 """ 270 271 10 169 16.9 0.0 @_decorate_validation 272 def validate_input(): 273 _numeric('y', ('integer', 'floating', 'complex'), shape=(-1, 1)) 274 _numeric('A', ('integer', 'floating', 'complex'), shape=( 275 y.shape[0], 276 _conf['true_solution'].shape[0] 277 if _conf['true_solution'] is not None else -1)) 278 279 if isinstance(A_asq, _MatrixBase): 280 # It is not possible to validate the range of an implicit matrix 281 # Thus allow all possible values 282 range_ = '[-inf;inf]' 283 else: 284 range_ = '[0;inf)' 285 _numeric('A_asq', ('integer', 'floating', 'complex'), range_=range_, 286 shape=A.shape, ignore_none=True) 287 288 10 155 15.5 0.0 @_decorate_validation 289 def validate_output(): 290 # complex128 is two float64 (real and imaginary part) each taking 8*8 291 # bits. Thus, in total 2*8*8=128 bits. However, we only consider it to 292 # be "64 bit precision" since that is what each part is. 293 bits_pr_nbytes = 4 if np.iscomplexobj(convert(0)) else 8 294 _numeric('alpha', ('integer', 'floating', 'complex'), 295 shape=(A.shape[1], 1), 296 precision=convert(0).nbytes * bits_pr_nbytes) 297 _generic('history', 'mapping', 298 keys_in=('alpha_bar', 'alpha_tilde', 'MSE', 299 'input_channel_parameters', 300 'output_channel_parameters', 'stop_criterion', 301 'stop_criterion_value', 'stop_iteration', 302 'stop_reason')) 303 304 10 35 3.5 0.0 validate_input() 305 306 # Initialisation 307 10 28254 2825.4 0.0 init = _get_gamp_initialisation(y, A, A_asq) 308 10 38 3.8 0.0 AH = init['AH'] 309 10 28 2.8 0.0 A_asq = init['A_asq'] 310 10 29 2.9 0.0 AT_asq = init['AT_asq'] 311 10 28 2.8 0.0 o = init['o'] 312 10 32 3.2 0.0 q = init['q'] 313 10 29 2.9 0.0 s = init['s'] 314 10 27 2.7 0.0 r = init['r'] 315 10 29 2.9 0.0 m = init['m'] 316 10 28 2.8 0.0 n = init['n'] 317 10 28 2.8 0.0 alpha_bar = init['alpha_bar'] 318 10 29 2.9 0.0 alpha_tilde = init['alpha_tilde'] 319 10 27 2.7 0.0 alpha_breve = init['alpha_breve'] 320 10 26 2.6 0.0 A_dot_alpha_bar = init['A_dot_alpha_bar'] 321 10 29 2.9 0.0 z_bar = init['z_bar'] 322 10 29 2.9 0.0 damping = init['damping'] 323 10 29 2.9 0.0 sum_approximation_method = init['sum_approximation_method'] 324 10 28 2.8 0.0 A_frob_sq = init['A_frob_sq'] 325 10 27 2.7 0.0 output_channel = init['output_channel'] 326 10 26 2.6 0.0 output_channel_parameters = init['output_channel_parameters'] 327 10 30 3.0 0.0 input_channel = init['input_channel'] 328 10 29 2.9 0.0 input_channel_parameters = init['input_channel_parameters'] 329 10 29 2.9 0.0 stop_criterion = init['stop_criterion'] 330 10 28 2.8 0.0 stop_criterion_name = init['stop_criterion_name'] 331 10 26 2.6 0.0 iterations = init['iterations'] 332 10 28 2.8 0.0 tolerance = init['tolerance'] 333 10 28 2.8 0.0 convert = init['convert'] 334 10 28 2.8 0.0 report_history = init['report_history'] 335 10 29 2.9 0.0 history = init['history'] 336 10 29 2.9 0.0 true_solution = init['true_solution'] 337 338 # GAMP iterations 339 5010 13798 2.8 0.0 for it in range(iterations): 340 # Save previous state 341 5000 12907 2.6 0.0 alpha_bar_prev = alpha_bar # Used in stop criterion 342 343 # GAMP state updates 344 5000 14069 2.8 0.0 if sum_approximation_method == 'rangan': 345 # Rangan's scalar variance sum approximation. 346 347 # Factor side updates 348 v = 1.0 / m * A_frob_sq * alpha_breve 349 o = A_dot_alpha_bar - v * q 350 z_bar, z_tilde = output_channel.compute(locals()) 351 q = (z_bar - o) / v 352 u = 1.0 / m * np.sum((v - z_tilde) / v**2) 353 354 # Variable side updates (including damping) 355 s_full = 1 / (1.0 / n * A_frob_sq * u) 356 s = (1 - damping) * s_full + damping * s 357 r = alpha_bar + s * AH.dot(q) 358 alpha_bar_full, alpha_tilde = input_channel.compute(locals()) 359 alpha_bar = (1 - damping) * alpha_bar_full + damping * alpha_bar 360 alpha_breve = 1.0 / n * np.sum(alpha_tilde) 361 362 else: 363 # Either "full" GAMP or Krzakala's sum approximation. 364 # For "full" GAMP, A_asq is the entrywise absolute value squared 365 # matrix. 366 # For Krzakala's sum approximation, A_asq is a 367 # magni.utils.matrices.SumApproximationMatrix that implements the 368 # scaled sum approximation to the full A_asq matrix. 369 370 # Factor side updates 371 5000 15683931 3136.8 22.8 v = A_asq.dot(alpha_tilde) 372 5000 180630 36.1 0.3 o = A_dot_alpha_bar - v * q 373 5000 1739912 348.0 2.5 z_bar, z_tilde = output_channel.compute(locals()) 374 5000 262888 52.6 0.4 q = (z_bar - o) / v 375 5000 319540 63.9 0.5 u = (v - z_tilde) / v**2 376 377 # Variable side updates (including damping) 378 5000 16601338 3320.3 24.1 s_full = 1 / (AT_asq.dot(u)) 379 5000 778120 155.6 1.1 s = (1 - damping) * s_full + damping * s 380 5000 6949608 1389.9 10.1 r = alpha_bar + s * AH.dot(q) 381 5000 18503827 3700.8 26.9 alpha_bar_full, alpha_tilde = input_channel.compute(locals()) 382 5000 743619 148.7 1.1 alpha_bar = (1 - damping) * alpha_bar_full + damping * alpha_bar 383 384 # Stop criterion 385 5000 6163475 1232.7 9.0 A_dot_alpha_bar = A.dot(alpha_bar) # Used in residual stop criteria 386 5000 770945 154.2 1.1 stop, stop_criterion_value = stop_criterion.compute(locals()) 387 388 # History reporting 389 5000 13536 2.7 0.0 if report_history: 390 history['alpha_bar'].append(alpha_bar) 391 history['alpha_tilde'].append(alpha_tilde) 392 history['input_channel_parameters'].append( 393 copy.deepcopy(input_channel.__dict__)) 394 history['output_channel_parameters'].append( 395 copy.deepcopy(output_channel.__dict__)) 396 history['stop_criterion_value'].append(stop_criterion_value) 397 history['stop_iteration'] = it 398 if true_solution is not None: 399 history['MSE'].append( 400 1/n * np.linalg.norm(true_solution - alpha_bar)**2) 401 402 5000 14355 2.9 0.0 if stop: 403 history['stop_reason'] = stop_criterion_name.upper() 404 break 405 406 10 26 2.6 0.0 alpha = alpha_bar 407 408 10 37 3.7 0.0 validate_output() 409 410 10 27 2.7 0.0 if report_history: 411 return alpha, history 412 else: 413 10 25 2.5 0.0 return alpha Total time: 18.271 s File: /user/cso/NoBackup/Magni/magni/cs/reconstruction/gamp/input_channel.py Function: compute at line 151 Line # Hits Time Per Hit % Time Line Contents ============================================================== 151 def compute(self, var): 152 """ 153 Compute the IIDBG input channel value. 154 155 Parameters 156 ---------- 157 var : dict 158 The variables used in computing of the input channel value. 159 160 Returns 161 ------- 162 mean : ndarray 163 The computed input channel mean. 164 variance : ndarray 165 The computed input channel variance. 166 167 """ 168 169 5000 105298 21.1 0.6 super(IIDBG, self).compute(var) 170 171 5000 5502 1.1 0.0 s = var['s'] 172 5000 3706 0.7 0.0 r = var['r'] 173 5000 4223 0.8 0.0 tau = self.tau 174 5000 3945 0.8 0.0 theta_bar = self.theta_bar 175 5000 4012 0.8 0.0 theta_tilde = self.theta_tilde 176 177 5000 1131524 226.3 6.2 h_pzc = (theta_tilde * s) / (s + theta_tilde) 178 5000 1172033 234.4 6.4 g_pzc = (theta_bar / theta_tilde + r / s) * h_pzc 179 5000 1604026 320.8 8.8 k_pzc = 1 + (1 - tau) / tau * np.sqrt(theta_tilde / h_pzc) * np.exp( 180 5000 9386440 1877.3 51.4 1/2 * ((r - theta_bar)**2 / (s + theta_tilde) - r**2 / s)) 181 5000 688537 137.7 3.8 l_pzc = 1.0 / k_pzc 182 183 # New values of alpha_bar (mean) and alpha_tilde (variance) 184 5000 288969 57.8 1.6 mean = l_pzc * g_pzc 185 5000 2202044 440.4 12.1 variance = l_pzc * (h_pzc + g_pzc**2) - l_pzc**2 * g_pzc**2 186 187 5000 9818 2.0 0.1 if self.use_em: 188 # EM update of tau 189 # See Eq. (19) in [3] 190 5000 211022 42.2 1.2 tau = 1.0 / self.n * np.sum(l_pzc) 191 5000 17803 3.6 0.1 self.tau = var['convert'](tau) 192 193 # EM update of theta_bar 194 # See Eq. (25) in [3] 195 5000 192234 38.4 1.1 theta_bar = 1.0 / (self.n * tau) * np.sum(mean) 196 5000 15840 3.2 0.1 self.theta_bar = var['convert'](theta_bar) 197 198 # EM update of theta_tilde 199 # See Eq. (32) in [3] 200 5000 8735 1.7 0.0 theta_tilde = 1.0 / (self.n * tau) * np.sum( 201 5000 1194598 238.9 6.5 l_pzc * ((theta_bar - g_pzc)**2 + h_pzc)) 202 5000 16426 3.3 0.1 self.theta_tilde = var['convert'](theta_tilde) 203 204 5000 4244 0.8 0.0 return mean, variance Total time: 1.62602 s File: /user/cso/NoBackup/Magni/magni/cs/reconstruction/gamp/output_channel.py Function: compute at line 153 Line # Hits Time Per Hit % Time Line Contents ============================================================== 153 def compute(self, var): 154 """ 155 Compute the AWGN output channel value. 156 157 Parameters 158 ---------- 159 var : dict 160 The variables used in computing of the output channel value. 161 162 Returns 163 ------- 164 mean : ndarray 165 The computed output channel mean. 166 variance : ndarray 167 The computed output channel variance. 168 169 """ 170 171 5000 104571 20.9 6.4 super(AWGN, self).compute(var) 172 173 5000 5225 1.0 0.3 y = var['y'] 174 175 5000 4302 0.9 0.3 if var['it'] != 0: # Allow for fixed initialisation of sigma_sq 176 # Noise level estimation for current iteration 177 4990 4898 1.0 0.3 if self.noise_level_estimation == 'sample_variance': 178 self.sigma_sq = 1.0/self.m * np.sum( 179 (y - var['A_dot_alpha_bar'])**2) 180 4990 3886 0.8 0.2 elif self.noise_level_estimation == 'median': 181 # Estimate variance based on median 182 # std(x) ~= median(|x|) / stdQ1 for x~N(0, std**2) 183 self.sigma_sq = ( 184 calculate_median(np.abs(y - var['A_dot_alpha_bar'])) / 185 self.stdQ1) ** 2 186 187 5000 77126 15.4 4.7 sigma_sq_plus_v = self.sigma_sq + var['v'] 188 189 5000 395990 79.2 24.4 mean = (var['v'] * y + self.sigma_sq * var['o']) / sigma_sq_plus_v 190 5000 246946 49.4 15.2 variance = (self.sigma_sq * var['v']) / sigma_sq_plus_v 191 192 5000 9141 1.8 0.6 if self.noise_level_estimation == 'em': 193 # EM noise level recursion 194 # See Eq. (77) in [1] 195 self.sigma_sq = self.sigma_sq * np.sum( 196 5000 466451 93.3 28.7 (np.abs(mean - var['o']) / var['v'])**2) / np.sum( 197 5000 291437 58.3 17.9 1.0 / sigma_sq_plus_v) 198 199 5000 12136 2.4 0.7 if len(variance.shape) != 2: 200 # When using a sum approximation of A_asq, the variance becomes 201 # scalar. However, our interface dictates that it must be a column 202 # vector for the sum approximation to be correct. Also, the 203 # denominator sum of the sigma_sq EM update must be scaled 204 # correctly as if sigma_sq_plus_v had been a vector. 205 206 variance = np.ones(mean.shape, dtype=var['convert']) * variance 207 self.sigma_sq = self.sigma_sq / len(y) 208 209 5000 3911 0.8 0.2 return mean, variance Total time: 0.665457 s File: /user/cso/NoBackup/Magni/magni/cs/reconstruction/gamp/stop_criterion.py Function: compute at line 124 Line # Hits Time Per Hit % Time Line Contents ============================================================== 124 def compute(self, var): 125 """ 126 Compute the MSE convergence stop criterion value. 127 128 The GAMP algorithm should converge to a fixed point. This criterion 129 is based on the mean squared error of the difference between the 130 proposed solution in this iteration and the proposed solution in the 131 previous solution. 132 133 Parameters 134 ---------- 135 var : dict 136 Dictionary of variables used in calculating of the stop criterion. 137 138 Returns 139 ------- 140 stop : bool 141 The indicator of whether or not the stop criterion is satisfied. 142 value : float 143 The stop criterion value. 144 145 """ 146 147 5000 100494 20.1 15.1 super(MSEConvergence, self).compute(var) 148 149 5000 8195 1.6 1.2 mse = 1/self.n * np.linalg.norm( 150 5000 545376 109.1 82.0 var['alpha_bar_prev'] - var['alpha_bar'])**2 151 5000 8429 1.7 1.3 stop = mse < self.tolerance 152 153 5000 2963 0.6 0.4 return stop, mse