STELA Algorithm

class fastmat.algorithms.STELA(fmatA, **kwargs)

Bases: Algorithm

Soft-Thresholding with simplified Exact Line search Algorithm (STELA)

The algorithm is presented in [1] with derivation and convergence results.

Definition and Interface: For a given matrix \(A \in \mathbb{C}^{m \times N}\) with \(m \ll N\) and a vector \(b \in \mathbb{C}^m\) we approximately solve

\[\min\limits_{ x \in \mathbb{C}^N}\Vert{ A \cdot x - b}\Vert^2_2 + \lambda \cdot \Vert x \Vert_1,\]

where \(\lambda > 0\) is a regularization parameter to steer the trade-off between data fidelity and sparsity of the solution.

>>> # import the packages
>>> import numpy.linalg as npl
>>> import numpy as np
>>> import fastmat as fm
>>> import fastmat.algorithms as fma
>>> # define the dimensions and the sparsity
>>> n, k = 512, 3
>>> # define the sampling positions
>>> t = np.linspace(0, 20 * np.pi, n)
>>> # construct the convolution matrix
>>> c = np.cos(2 * t) * np.exp(-t ** 2)
>>> C = fm.Circulant(c)
>>> # create the ground truth
>>> x = np.zeros(n)
>>> x[np.random.choice(range(n), k, replace=0)] = 1
>>> b = C * x
>>> # reconstruct it
>>> stela = fma.STELA(C, numLambda=0.005, numMaxSteps=100)
>>> y = stela.process(b)
>>> # test if they are close in the
>>> # domain of C
>>> print(npl.norm(C * y - b))

We solve a sparse deconvolution problem, where the atoms are harmonics windowed by a gaussian envelope. The ground truth \(x\) is build out of three pulses at arbitrary locations.

Note

The proper choice of \(\lambda\) is crucial for good perfomance of this algorithm, but this is not an easy task. Unfortunately we are not in the place here to give you a rule of thumb what to do, since it highly depends on the application at hand. Again, consult [1] for any further considerations of this matter.

1(1,2)

Y. Yang, M. Pesavento, “A Unified Successive Pseudoconvex Approximation Framework”, IEEE Transactions on Signal Processing, vol. 65, no. 13, pp. 3313-3327, Dec 2017

Parameters

fmatAfm.Matrix

the system matrix

arrBnp.ndarray

the measurement vector

numLambdafloat, optional

the thresholding parameter; default is 0.1

numMaxStepsint, optional

maximum number of steps; default is 100

numMaxErrorfloat, optional

maximum error tolerance; default is 1e-6

Returns

np.ndarray

solution array

__init__(fmatA, **kwargs)

softThreshold(arrX, numAlpha)

Do a soft thresholding step.