magni.cs.reconstruction.amp._algorithm module¶

Module providing the core Approximate Message Passing (AMP) algorithm.

Routine listings¶

run(y, A)
Run the AMP reconstruction algorithm.

Notes

The default configuration of the AMP algorithm provides the Donoho, Maleki, Montanari (DMM) AMP from [1] using the soft threshold with the “residual” thresholding strategy from [2].

References

 [1] 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, p. 18914-18919, Nov. 2009.
 [2] A. Montanari, “Graphical models concepts in compressed sensing” in Compressed Sensing: Theory and Applications, Y. C. Eldar and G. Kutyniok (Ed.), Cambridge University Press, ch. 9, pp. 394-438, 2012.
magni.cs.reconstruction.amp._algorithm.run(y, A)[source]

Run the AMP reconstruction algorithm.

Parameters: y (ndarray) – The m x 1 measurement vector. A (ndarray or magni.utils.matrices.{Matrix, MatrixCollection}) – The m x n matrix which is the product of the measurement matrix and the dictionary matrix. alpha (ndarray) – The n x 1 reconstructed coefficient vector. history (dict, optional) – The dictionary of various measures tracked in the GAMP iterations.

Notes

Optionally, the algorithm may be configured to save and return the iteration history along with the reconstruction result. The returned history contains the following:

• alpha_bar : Coefficient estimates (the reconstruction coefficients).
• MSE : solution mean squared error (if the true solution is known).
• threshold_parameters : The state of the threshold parameters.
• stop_criterion : The currently used stop criterion.
• stop_criterion_value : The value of the stop criterion.
• stop_iteration : The iteration at which the algorithm stopped.
• stop_reason : The reason for termination of the algorithm.

Examples

For example, recovering a vector from noiseless measurements using AMP with soft thresholding

>>> import numpy as np
>>> from magni.cs.reconstruction.amp import run, config
>>> from magni.cs.reconstruction.amp.util import theta_mm
>>> np.random.seed(seed=6028)
>>> k, m, n = 10, 200, 400
>>> A = 1 / np.sqrt(m) * np.random.randn(m, n)
>>> alpha = np.zeros((n, 1))
>>> alpha[:k] = np.random.normal(scale=2, size=(k, 1))
>>> np_printoptions = np.get_printoptions()
>>> np.set_printoptions(suppress=True, threshold=k+2)
>>> alpha[:k + 2]
array([[ 1.92709461],
[ 0.74378508],
[-3.2418159 ],
[-1.32277347],
[ 0.90118   ],
[-0.19157262],
[ 0.82855712],
[ 0.24817994],
[-1.43034777],
[-0.21232344],
[ 0.        ],
[ 0.        ]])
>>> y = A.dot(alpha)
>>> threshold_params = {'threshold_level_update_method': 'residual',
... 'theta': theta_mm(float(m) / n), 'tau_hat_sq': 1.0}
>>> config['threshold_parameters'] = threshold_params
>>> alpha_hat = run(y, A)
>>> alpha_hat[:k + 2]
array([[ 1.92612707],
[ 0.74185284],
[-3.24179313],
[-1.32200929],
[ 0.90089208],
[-0.19097308],
[ 0.82658038],
[ 0.24515825],
[-1.42980997],
[-0.2111469 ],
[ 0.00002815],
[ 0.00047293]])
>>> np.sum(np.abs(alpha - alpha_hat) > 1e-2)
0

or recover the same vector returning a history comparing the pr. iteration solution to the true vector

>>> config['report_history'] = True
>>> config['true_solution'] = alpha
>>> alpha_hat, history = run(y, A)
>>> alpha_hat[:k + 2]
array([[ 1.92612707],
[ 0.74185284],
[-3.24179313],
[-1.32200929],
[ 0.90089208],
[-0.19097308],
[ 0.82658038],
[ 0.24515825],
[-1.42980997],
[-0.2111469 ],
[ 0.00002815],
[ 0.00047293]])
>>> np.array(history['MSE']).reshape(-1, 1)[1:11]
array([[ 0.01273323],
[ 0.0053925 ],
[ 0.00276334],
[ 0.00093472],
[ 0.0004473 ],
[ 0.00021142],
[ 0.00009618],
[ 0.00004371],
[ 0.00002275],
[ 0.00001053]])

or recover the same vector using the median based update method

>>> config['report_history'] = False
>>> threshold_params['threshold_level_update_method'] = 'median'
>>> config['threshold_parameters'] = threshold_params
>>> alpha_hat = run(y, A)
>>> alpha_hat[:k + 2]
array([[ 1.92577079],
[ 0.7413161 ],
[-3.24131125],
[-1.32044146],
[ 0.89998407],
[-0.18981235],
[ 0.82584016],
[ 0.24445332],
[-1.42934708],
[-0.21043221],
[ 0.        ],
[ 0.        ]])
>>> np.sum(np.abs(alpha - alpha_hat) > 1e-2)
0
>>> np.set_printoptions(**np_printoptions)
magni.cs.reconstruction.amp._algorithm._get_amp_initialisation(y, A)[source]

Return an initialisation of the AMP algorithm.

Parameters: y (ndarray) – The m x 1 measurement vector. A (ndarray or magni.utils.matrices.{Matrix, MatrixCollection}) – The m x n matrix which is the product of the measurement matrix and the dictionary matrix. init (dict) – The initialisation of the the AMP algorithm.