Multiple measurement vector (MMV) problems occur in many applications, for example, in tomography [12], communication [13], blurred image reconstruction, or superresolution [14]. In the following, we derive the connection to block-sparsity. In an MMV problem, we assume that we derived d ∈ ℕ measurements y^l ∈ ℝ^m, l = 1, …, d from d sparse signal vectors x^l ∈ ℝⁿ sharing the same support supp(xl)={i:|xil|≠0}, which is referred to as joint-sparsity [15, 16]. The MMV problem can then be presented as solving the following d equations

for l = 1, …, d and with K ∈ ℝ^m×n. This can be rewritten in the following matrix equation form

where X = (x¹, …, x^d) ∈ ℝ^n×d, Y = (y¹, …, y^d) ∈ ℝ^m×d. With the vectorizing operator vec(·), stacking each column of a matrix on each other, we can cast Equation (2) into a block-sparse problem. We have

where we used the well-known vectorization property of matrix equations, see for example Schacke [17]. Here, ⊗ is the Kronecker product. The vector x = vec(X^T) is block-sparse with n blocks of length d, if the signals x^l are jointly sparse. With D=(K⊗Id)∈ℝny×nx, where n_y = m · d and n_x = n · d, we obtain the block-sparse setting considered in this work.

In the more general setting, we want to reconstruct an unknown vector x∈ℝnx from a given matrix D∈ℝny×nx and given y∈ℝny

with n_x = nd, n_y = md for some n, m, d ∈ ℕ. We assume that noise ϵ∈ℝny is added to Dx. In applications, we often have this problem is ill-posed, i.e., n_y < n_x or D is not invertible. We assume that x is the concatenation of n “smaller” vectors of length d, called blocks, i.e., x[i] ∈ ℝ^d,

Following the notation in Yonina and Eldar [18], we define

where I(‖x[i]‖₂ > 0) = 1 if ‖x[i]‖₂ > 0 and equal to zero otherwise and the ℓ₂-norm is defined as ‖x‖22=∑i=1n|xi|2. We call x∈ℝnx s-block-sparse if ‖x‖_2,0 ≤ s. Similar to Equation (4), we can construct the matrix D in Equation (3) from n “smaller” matrices D[i]∈ℝny×d, i = 1, …, n

Without loss of generality, we can assume that these blocks are orthonormal, see Yonina and Eldar [18], i.e., D[i]TD[i]=Id, where I_d is the d × d identity matrix. This assumption simplifies the presentation of several statements below.

Definition 1. The block-coherence of a matrix D∈ℝny×nx is defined as

Here, we use ‖A‖2=λmax(ATA), where λ_max denotes the largest eigenvalue of matrix A^TA. Note that the block coherence can also be introduced with a normalization factor 1/(‖D[i]‖₂‖D[j]‖₂), but since we assume orthonormal blocks, we can neglect this. This reduces to the already known coherence μ for d = 1

where D_{:, i} denotes the ith column of D, see Donoho et al. [19]. In Yonina and Eldar [18], it is shown that 0 ≤ μ_b(D) ≤ μ(D) ≤ 1, and it is possible to derive recovery statements for small μ_b similar to μ, for details see Yonina and Eldar [18]. We can consider also the cross block coherence, which compares two matrices and will be important in the following.

Definition 2. For B,D∈ℝny×nx with B[i]TD[i]=Id the cross block coherence is defined as

Similar to ordinary basis pursuit and LASSO [5, 20–22], we use the following ℓ_2,1-LASSO to solve Equation (3)

where ‖x‖2,1=∑i=1n‖x[i]‖2 is the ℓ_2,1-norm of x, which will promote block-sparsity for the solution of Equation (7). This convex programm can be solved by the fixed point iteration known as the block iterative shrinkage thresholding algroithm (Block-ISTA/BISTA):

where η_α is the block-soft-thresholding operator, given as

BISTA is an already well-studied proximal gradient algorithm based on the functional in Equation (7). It is known that the fixed point iteration in Equation (8) converges for γ∈(0,1L), where L=‖D‖22 is the Lipschitz constant of the least squares term in Equation (7) w.r.t to x, to a solution of Equation (7), if at least one solution exists, see for example Byrne [23], Bauschke et al. [24], and Beck [25]. On the other hand, the choice of the regularization parameter α has to be done empirically and is very crucial for a “good” recovery. If α is set too large, this can lead to too much damping, possibly setting blocks to 0 that actually have a non-zero norm. If α is too small, we get reverse effects. In practice, this leads to problems since computing the iterations require high computational effort. Deep unfolding is a way to tackle these problems, i.e., reduce the number of iterations by learning optimal regularization parameters and step-sizes. There are already classical concepts in increasing the convergence speed by using an additional step in updating the current iterate x^(k), by using the previous x^(k−1), resulting in Block Fast ISTA [25, 26]. On the other hand, the choice of optimal parameters is still solved empirically.

We will now formalize the concept of deep unfolding for an arbitrary operator that depends on a certain set of parameters, before we apply this to the previously presented fixed point iteration. To this end, we define an operator

which depends on a set of parameters θ ∈ Θ, such as the stepsize of a gradient descent operator and an input y. For example, for BISTA, this would be θ = (α, γ) with

We assume that Fix(T(·;θ, y)) ≠ ∅ and that we have convergence for the fixed point iteration

for an arbitrary x⁽⁰⁾ ∈ X. Deep unfolding now interprets each iteration step as the layer of a neural network and uses the parameters θ ∈ Θ as trainable variables. In more detail, this means we look at the Kth iteration of Equation (11), i.e., the composition

for some K ∈ ℕ. By unfolding this iterative scheme, we define

for k = 1, …, K and set θ^(k) as the set of trainable variables in this operator, so they can vary in each iteration step. For example, in LBISTA, we will get θ^(k) = (α^(k), γ^(k)) in the later called tied case. With this we get the following composition

and define the operator

where we get the full parameter space Θ~=Θ×⋯×Θ, with trainable variables θ~=⋃i=1Kθ(i). TΘ~ will then be the neural network which will be trained with respect to θ~. So in the end, after training, we have an iterative algorithm with fixed but optimized parameters. It seems thus that deep unfolding can be applied to any iterative algorithm and help us to estimate the best choice of parameters, but we will only present deep unfolding for BISTA and consider deep unfolding for arbitrary operators in future studies.

In this section, we give an overview for the training procedure used in this study. The general idea of supervised learning is to choose model parameters such that the predictions are close, in some sense, to the unknown target, i.e., in our case, the unknown vector x in Equation (3) generating the measurement y. Hence, we aim to minimize an expected loss over an unknown distribution D:

where ℓ(·) is a given loss function, x^=Tθ~(x(0);y) is the output of the model, and x^* is the ground-truth. Here, we will use the squared ℓ₂-loss ℓ(x)=12‖x‖22. The objective functional R(θ)=𝔼x*˜D[ℓ(x^-x*)] is also called risk of the model. Since the underlying distribution D is unknown, we take a batch of S independently drawn samples of input and output data (xj*,yj) for j = 1, …, n_train according to Equation (3) and minimize instead data-driven the empirical risk

Proceeding in this way for all layer at once is sometimes referred as end-to-end learning. Because of the special structure of our deep unfolding models, and inspired by Musa et al. [29], we instead train the network layer-wise, by optimizing only θcase(k-1) for layer k yielding the following training procedure: Let k ∈ {1, …, K},

where x^(k) is the output of the kth layer. We realized this training as follows, we generate a validation set (xi,validation*,yi,validation), used to evaluate the model while training and a training set (xi,train*,yi,train), i = 1, …, n_train, used to calculate (15). This objective is locally minimized by gradient descent methods. As a stopping criteria, we evaluate the normalized mean square error, defined as

depending on the validation set, and stop if the maximum of all evaluated NMSE stays the same for a given number of iterations. See Algorithm 1, where Adam is the ADAM Optimizer [30] depending on an training rate t_r and the functional which should be minimized with respect to given variables, here the loss function ℓ(·) with respect to θ^k−1.

In the following, we present four different unfolding techniques for BISTA. We present a tied (weights are shared between different layers) and untied (individual weights per layer) case, which refers to different training approaches

Tied LBISTA: The idea of LBISTA is now to fix the matrices S and B for all layers but also include them in our set of trainable variables:

For LBISTA (Equation 16), we get trainable variables

where we initialize S = I − B^TD and B = γD. Algorithm (16) is also referred to as vanilla LISTA in the sparse case. Inspired by the LISTA-CP model, i.e., LISTA with coupled parameters, proposed in Liu and Chen [10], we will also consider LBISTA-CP

For LBISTA-CP (Equation 17), we get

where we initialize B = D.

Untied LBISTA: The idea of untied LBISTA is then to use in each layer different matrices S and B to train, i.e.,

For LBISTA (untied) Equation (18), we get trainable variables

where we initialize S^(k) = I − B^TD and B^(k) = γD for every k = 0, …, K − 1. Inspired by Algorithm (17), we will also consider Algorithm (19), which will be referred to as LBISTA-CP (untied), [9, 10]:

For LBISTA-CP (untied) (Equation 19), we get

where we initialize B^(k) = D for every k = 0, …, K − 1. Hence, compared to O(nynx+K) parameters in the tied case, now more training data and longer training time is required to train now O(Knynx+K) parameters.

We initialize the trainable variables with values from original BISTA. In Chen et al. [9], it has been shown that convergence of LISTA-CP (untied) can be guaranteed if the matrices B^(k) belong to a certain set and their proof can be extended to block-sparsity. The steps are very similar to the convergence proof for learned block analytical ISTA given in the next section.

This part of the study will focus on combining and extending several theoretical statements from Chen et al. [9] and Liu and Chen [10] and applying these for the block-sparse case. With the following two theorems, we are then going to present analytical LBISTA, by showing that this is as good as LBISTA-CP (untied) with a pre-computed B.

Analogously to Liu and Chen [10] and following the previous two theorems, decompose LBISTA-CP (untied), Algorithm 17, into two steps:

where, in the first step, B~ is pre-computed, such that

In the second step, the parameters θ=((α(k))k=0K-1,(γ(k))k=0K-1) are trained layer wise, as discussed in the previous section. This results in a comparable method, with only O(K) trainable parameters, instead of O(nynx+K) for LBISTA-CP (Equation 17) or even O(Knynx+K) for LBISTA (untied) Equation 18.

Since the objective in Equation (21) is not differentiable, one solves the following upper bound problem

This is derived from the following inequality

In Liu and Chen [10], this is solved by a projected gradient method, but since this is a constrained linear least squares problem in the vector space of matrices with Frobenius inner product the following Theorem states a closed form of the solution of Equation (31).

Theorem 3. The minimizer B∈ℝny×nx of Equation (31) is given as the concatenation

where the n blocks are given as

with

for i = 1, …, n.

The proof can be found in Appendix C.

Let d = 1 and the singular value decomposition of D given as D = VΣU^T and assume B=VΣ~UT. Then, the solution of Equation (31) is given in an even simpler form since

Choosing B=D+,Tdiag(d~)-1, where d~=diag(D+D), i.e., a normalized pseudo-inverse, yields also a solution for Equation (31). Here, diag follows the matlab/python notation, where diag of a matrix gives the vector of the main diagonal and gives a diagonal matrix with a given vector on its main diagonal. For d ≥ 2, the orthonormal block constraints in Equation (31) would not be met since with this construction we can only guarantee that the diagonal elements of B^TD are equal to one but cannot control what happens on the off-diagonal, thus not yielding a feasible solution.

In practice, often large data sets are obtained, i.e., by a large amount of measurements or measurements y^l, x^l representing pictures. For instance, 1, 000 pixels and 200 measurements lead to a matrix D with (1, 000·200)² elements. Applying the theory for analytical LBISTA could thus be difficult in practice. Although D is sparse, it can take a long time to calculate the analytical weight matrix. The following Theorem states the connection between the MMV setting and the block sparse setting for ALBISTA, showing that it is sufficient to minimize the generalized mutual coherence Equation (32), see Liu and Chen [10], for K instead of minimizing the generalized mutual block-coherence (Equation 21) for D = K ⊗ I_d. Moreover, in many applications the measurement model involves convolutions, see for example Ahmadi et al. [14], so we conclude this section by considering the cases where K is a circular or Toeplitz matrix.

Theorem 4. Let B~ attain the minimum in

where B~:,i refers to the ith column. Then the minimum of

is given as B=B~⊗Id, if D = K ⊗ I_d.

The proof can be found in Appendix D. Moreover, the following relation holds, if D = K ⊗ I_d,

thus, the block-coherence of D can be enhanced by increasing the number of measurements d. It is feasible to solve Equation (31) by the pseudo inverse in the MMV setting, since this is solved for K, i.e., d = 1

Consider now the following setting, where the measurements y^l ∈ ℝⁿ are obtained by a circular convolution of x^l with vector k, i.e.,

where K is a circular matrix generated by a vector k ∈ ℝⁿ, i.e., K = circ(k) ∈ ℝ^n×n. Applying Theorem 3 to the circular case yields the following lemma.

Lemma 1. Let k ∈ ℝⁿ and let B = circ(b) ∈ ℝ^n×n where b ∈ ℝⁿ is given by

λ ∈ ℝ. Then, B attains the minimum in Equation (31).

The latter statement implies a simpler way to compute b ∈ ℝⁿ by using singular value decomposition

where σ(B) ∈ ℝⁿ is the vector of singular values of B, respectively K, filled with zeros and U an unitary matrix. With U=1/nF, where F is the Discrete Fourier Transform (DFT) matrix, this leads to the conclusion σ(B)=Fb=b^, and thus

where k^=Fk. The expression 1/k^ should be interpreted point wise and to be zero if k^i=0. This also concludes that the computation of B tends to be difficult in practice if K has not full rank, since this means that k^ has at least one zero entry. In this case, b has to be scaled with nrank(K) since

Considering the more general convolutional setting,

where k∈ℝm~ and x ∈ ℝⁿ with m~<n. This results in a Toeplitz matrix K

The reasoning of the previous section cannot be applied to show that the solution of Equation (31) must be a Toeplitz matrix. But the following can be observed: Let b be constructed as discussed for K~=circ(k), where k is the concatenation of k and a zero vector of suitable dimensions, i.e., the first column of K. The analytical weight matrix B w.r.t. K can be constructed as

where T is the m × m cyclic shift matrix, i.e., only the m × n submatrix of B~=circ(b)∈ℝm×m is used. The columns of K can also be expressed through the cyclic shift of k. Hence,

i.e., B is a feasible solution of Equation (31) for K. On the other hand, the cross coherence is bounded since

where ‖D‖_max = max_i,j|d_ij| is the maximum norm and diag(σ(B)⊙σ(K)) = I_m derives from Equation (35). Note: To have this upper bound K~ needs to have full rank, which is the case if the full time continuous FT of k has no zero points. Or one has to adjust the discrete grid. Thus, constructing B by extending K to a circular matrix is a feasible approach.

It is known that using the Fast Fourier Transform (FFT) in CNNs can decrease the computation time if the convolutional filter is big [31, 32]. In Pratt et al. [31] showed that training weights in the Fourier domain can reduce training time while maintaining efficiency. By using this, the costs of O(n2) operations could be reduced to O(nlog(n)) operations. To connect the theory of FFT-CNNs and of unrolling ISTA in the context of deconvolution (Equations 33 or 36), the gradient step can be viewed as follows:

This can be interpreted as a convolutional layer with kernel f(γ) = (e − γb * k) and bias b~=b*y, where e = [1, 0, …, 0]. This means that ALISTA, with a Toeplitz matrix, can be interpreted as a CNN only two trainable parameters per layer, γ, λ. In the setting of FFT-CNN, the update rule can be formulated as

On the other hand, using FFT CNNs shows only a speed up if we deal with large data sets, i.e. by evaluating high-resolution images, or large filters, i.e., if m~ is greater than log(n).

The training data are sampled from an unknown distribution X and generated as follows. The signals x are generated for a given number of blocks n, given block length d and a possibility if a block x[i] is active or not, i.e., if ‖x[i]‖₂ ≠ 0 or ‖x[i]‖₂ = 0, called pnz (probability of non-zeros). If a block is active, the elements of this blocks are given by a normal Gaussian distribution with variance σ² = 1. The measurements y are obtained by Equation (3), where the elements of ϵ are given from a normal Gaussian distribution with variance σ2=pnz·nx/ny·10-SNRdB/10. SNR_dB is the signal-to-noise ratio given in decibel. We consider the following cases. In each case, we generate x with d = 15, n = 128, m = 32, and pnz = 10%.

Gaussian measurement matrix: In the Gaussian setting, we sample a m × n matrix K iid from a Gaussian distribution with variance σ² = 1. We normalize the columns, s.t. D has orthonormal blocks, as assumed in the beginning.

Circular convolution matrix: We construct the circular matrix as follows. At first, we generate a random iid sampled vector ã but set a certain amount of elements to zero. We define k = ℜ(F⁻¹ã), and thus we can generate a rank deficient Matrix D = K ⊗ I_d, where K = circ(k). Thus, y^l is obtained through a circular convolution with symmetric k^=Fk. It is important for this section to generate a rank deficient matrix to have compressive observations. Otherwise, we would get a trivial problem if K, respectively D, has full rank. Computing D = K ⊗ I_d yields the desired matrix. The properties of these two matrices can be found in Table 1.

The results of the proposed methods can be found in Figure 1A for the Gaussian Problem case in Figure 1B for the circular case. We also consider the performance of a version of AMP [33]. AMP can be viewed as an Bayesian extension of ISTA with an additional Onsager correction term, before applying the thresholding operator, i.e.,

where b^(k) = 𝔼[η′(x^(k))]. Here, we will train only γ, α, with the same procedure as already discussed and choose D^T = B^T. By using B^T instead of D^T we are resembling Orthogonal AMP (OAMP) [34], which follows a similar idea as ALISTA. In Ma and Ping [34], B is chosen to be de-correlated with respect to K, i.e.,

The analytical weight matrix B satisfies also this condition. Moreover, in Ma and Ping [34], the choice of different matrices is also discussed. Thus, untrained ALISTA with correction term can be viewed as a special case of OAMP. Note also, that the structure of Trainable ISTA, proposed in Ito et al. [35], is also based on OAMP and thus there are interrelations between AMP, unfolding ISTA, and analytical ISTA. We refer to learned AMP with analytical matrix as ALAMP. Different from Ma and Ping [34], we use the ℓ_2,1-regularizer, instead of only ℓ₁-regularization, since we consider the MMV setting. This is also discussed in [36, 37]. Every proposed algorithm performs better, in terms of NMSE, as their untrained original. Interestingly learned AMP, with analytical weight matrix B, has a performance almost as good as LBISTA (untied) with only a fraction of trainable parameters. This may come from the fact, that block-soft thresholding is not the correct MMSE estimator for the generated signals x and thus the correction of b^(k) yields the better estimation for x. As expected, we get almost the same performance of ALBISTA and LBISTA CP (untied). In Figure 2A, we show similar plots for the justification of Theorem 1, as also seen in Liu and Chen [10]. In particular, Figure 2A shows that α(k)γ(k) is proportional to the maximal ℓ_2,1-error over all training signals. An interesting behavior, which carries over from the sparse case, is that the learned ℓ_2,1-regularization parameters approach zero, as k increases. If α is close to zero, we approach a least-squares problem. This means that after LBISTA found the support of the unknown signal x^* the algorithm consist only of the least squares fitting. Figure 2B shows that the trained γ^(k) are bound in an interval. Note that, in contrast to Liu and Chen [10], Theorem 1 is based on a more general assumption, onto the thresholding parameters. One can take a suitable κ and obtains the upper bound for the ℓ₂-error and thus have convergence on the training set if the sparsity assumptions are met. Figure 3 shows the training loss over the training iterations for the results presented in Figure 1. One can see that ALBISTA needs less training iterations as LBISTA CP (untied) or LBISTA (untied) and thus less training data. The observed jumps occur when moving from one layer to the next due to the layer-wise training, Algorithm 1. A layer is defined to be optimized if NMSE converged within 1e − 5.