magni.cs.reconstruction.gamp.channel_initialisation module

Module providing utility functions for initilisation of in- and output channels in the Generalised Approximate Message Passing (GAMP) algorithm.

Routine listings

get_em_bg_amp_initialisation(problem_params, method=’vila’)
Get initial parameters for EM Bernoulli-Guassian AMP.
rho_se(delta, zeta, resolution=1000)
Return the theoretical noiseless LASSO phase transition.
magni.cs.reconstruction.gamp.channel_initialisation.get_em_bg_amp_initialisation(problem_params, method='vila')[source]

Get initial parameters for EM Bernoulli-Guassian AMP.

If the initialisation method is vila then the scheme from [1] is used. If it is krzakala then the scheme from [2] is used.

Parameters:
  • problem_params (dict) – The problem parameters used to compute the initialisation.
  • method (str) – The initialisation method to use.
Returns:

  • tau (float) – The initial sparsity level.
  • sigma_sq (float) – The initial AWGN noise level.
  • theta_tilde (float) – The initial Gaussian prior variance.

See also

magni.cs.reconstruction.gamp._output_channel.wrap_calculate_using_AWGN()
Related output channel.
magni.cs.reconstruction.gamp._input_channels.wrap_calculate_using_iidsGB()
Related input channel.

Notes

Independently of the choice of method, the problem_params are:

  • y: the measurements
  • A: the system matrix

If method is vila, one must also specify:

  • SNR: the signal-to-noise ratio

References

[1]J. P. Vila and P. Schniter, “Expectation-Maximization Gaussian-Mixture Approximate Message Passing”, IEEE Transactions on Signal Processing, 2013, vol. 61, no. 19, pp. 4658-4672, Oct. 2013.
[2]F. Krzakala, M. Mezard, F. Sausset, Y. Sun, and L. Zdeborova, “Probabilistic reconstruction in compressed sensing: algorithms, phase diagrams, and threshold achieving matrices”, Journal of Statistical Mechanics: Theory and Experiment, vol. P08009, pp. 1-57, Aug. 2012.

Examples

For example, get the “vila” initialisation for a SNR of 100

>>> import numpy as np
>>> from magni.cs.reconstruction.gamp.channel_initialisation import     ... get_em_bg_amp_initialisation as get_init
>>> np.random.seed(6012)
>>> A = np.random.randn(20, 40)
>>> y = np.random.randn(20, 1)
>>> problem_params = {'A': A, 'y': y, 'SNR': 100}
>>> init_tup = get_init(problem_params, method='vila')
>>> [round(float(elem), 3) for elem in init_tup]
[0.193, 0.01, 0.123]

or get the corresponding “krzakala” initialisation

>>> del problem_params['SNR']
>>> init_tup = get_init(problem_params, method='krzakala')
>>> [round(float(elem), 3) for elem in init_tup]
[0.05, 1.0, 0.479]
magni.cs.reconstruction.gamp.channel_initialisation.rho_se(delta, zeta, resolution=1000)[source]

Return the theoretical noiseless LASSO phase transition.

Parameters:
  • delta (float) – The under sampling ratio.
  • zeta ({1, 2}) – The “problem” to get the phase transition for.
  • resolution (int) – The resolution used in the brute force optimisation.
Returns:

rho_se (float) – The phase transition value.

Notes

The theoretical noiseless LASSO phase transition is computed based on eq. 5 in [3]. A simple brute force optimisation with the specified resolution of that expression is used to find the phase transition. The “problems”, for which the phase transition may be computed, are:

  1. Sparse nonnegative vectors
  2. Sparse signed vectors

References

[3]D.L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing”, Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 45, pp. 18914-18919, Nov. 2009.

Examples

For example, find a phase transition value for a sparse signed vector:

>>> from magni.cs.reconstruction.gamp.channel_initialisation import rho_se
>>> round(float(rho_se(0.19, 2)), 3)
0.238

or find the corresponding value for a sparse nonnegative vector

>>> round(float(rho_se(0.19, 1)), 3)
0.318

and a few more examples

>>> round(float(rho_se(0.0, 1)), 3)
0.0
>>> round(float(rho_se(0.0, 2)), 3)
0.0
>>> round(float(rho_se(0.5, 1)), 3)
0.558
>>> round(float(rho_se(0.5, 2)), 3)
0.386
>>> round(float(rho_se(1.0, 1)), 3)
0.95
>>> round(float(rho_se(1.0, 2)), 3)
0.95