Trainable ISTA for Sparse Signal Recovery
Daisuke Ito∗ , Satoshi Takabe∗† , and Tadashi Wadayama∗
∗ Nagoya

arXiv:1801.01978v1 [] 6 Jan 2018

Institute of Technology, Gokiso, Nagoya, Aichi, 466-8555, Japan,
d.ito.480@stn.nitech.co.jp, {s_takabe, wadayama}@nitech.ac.jp
† RIKEN Center for Advanced Intelligence Project, Nihonbashi, Chuo-ku, Tokyo, 103-0027, Japan,

Abstract—In this paper, we propose a novel sparse signal
recovery algorithm called Trainable ISTA (TISTA). The proposed
algorithm consists of two estimation units such as a linear
estimation unit and a minimum mean squared error (MMSE)
estimator-based shrinkage unit. The estimated error variance
required in the MMSE shrinkage unit is precisely estimated
from a tentative estimate of the original signal. The remarkable
feature of the proposed scheme is that TISTA includes adjustable
variables controlling a step size and the error variance for the
MMSE shrinkage. The variables are adjusted by standard deep
learning techniques. The number of trainable variables of TISTA
is equal to the number of iteration rounds and it is much
smaller than that of known learnable sparse signal recovery
algorithms. This feature leads to highly stable and fast training
processes of TISTA. Computer experiments show that TISTA is
applicable to various classes of sensing matrices such as Gaussian
matrices, binary matrices and matrices with large condition
numbers. Numerical results also demonstrate that TISTA shows
significantly faster convergence than those of AMP and LISTA
in many cases.

I. I NTRODUCTION
The basic problem setup for compressed sensing [1], [2] is
as follows. A real vector x ∈ RN represents the source sparse
signal. It is assumed that we cannot directly observe x but we
observe y = Ax+w where A ∈ RM ×N (N ≥ M ) is a sensing
matrix and w ∈ RM is a Gaussian noise vector. Our goal is
to estimate x from y as correct as possible.
For a number of sparse reconstruction algorithms [3], Lasso
formulation [4] is fairly common for solving the sparse signal
recovery problems. In Lasso formulation, the original problem
is recast as a a convex optimization problem for minimizing
1
2
2 ||y − Ax||2 + λ||x||1 . The regularization term λ||x||1 promotes sparseness of a reconstruction vector where λ is the
regularization constant. In order to solve the Lasso problem
efficiently, a number of algorithms have been developed [5].
Iterative Shrinkage Thresholding Algorithm (ISTA) [6], [7]
is one of the most known algorithms for solving the Lasso
problem. ISTA is an iterative algorithm comprising of two
processes, i.e., a linear estimation process and a shrinkage
process based on a soft thresholding function. ISTA can be
seen as a proximal gradient descent algorithm [8] and it can
be directly derived from Lasso formulation.
Approximate Message Passing (AMP) [9], [10] which is a
variant of approximate belief propagation, shows much faster
convergence than that of ISTA in general. The remarkable
feature of AMP is that its asymptotic behavior is completely
described by the state evolution equations [11]. AMP is
derived based on the assumption that the sensing matrices

consist of i.i.d. Gaussian distributed components. Recently,
Ma and Ping proposed Orthogonal AMP (OAMP) [13] that can
handle various classes of sensing matrices including unitary invariant matrices. Rangan et al. proposed VAMP [14] for rightrotationally invariant matrices and provided a theoretical justification for its state evolution. Independently, Takeuchi [15]
also gave a rigorous analysis for a sparse recovery algorithm
for unitary invariant measurements based on the expectation
propagation framework.
The advent of recent powerful neural networks (NN) triggered explosive spread of research activities and applications
on deep neural networks (DNN) [16]. The DNN have found
a number of practical applications such as image recognition [17], [18], speech recognition [19] and robotics because
of their outstanding performance compared with the traditional methods. The advancement of DNN are also giving
impact on design of algorithms for communications and signal
processing [20], [21]. By unfolding an iterative process of a
sparse signal recovery algorithm, we can obtain a signal-flow
graph. The signal-flow graph includes trainable variables that
can be tuned with a supervised leaning method, i.e., standard
deep learning techniques such as stochastic gradient descent
algorithms based on the back propagation and mini-batchs can
be used to adjust the trainable parameters. Gregor and LeCun
presented Learned ISTA (LISTA) [22] that employs learnable
threshold variables for a shrinkage function. It is shown that
LISTA yields recovery performance superior to that of the
original ISTA. Borgerding et al. also presented variants of
AMP and VAMP with learnable capability [23] [24].
The goal of this work is to propose a simple sparse recovery
algorithm based on deep learning techniques. The proposed
algorithm, called Trainable ISTA (TISTA), borrows the basic
structure from ISTA, and an estimator of the squared error
between true signals and tentative estimations, i.e., the error
variance estimator, from OAMP. Thus, TISTA comprises
of the three parts: a linear estimator, an minimum mean
squared error (MMSE) estimator-based shrinkage function,
and the error variance estimator. The linear estimator of TISTA
includes trainable variables that can be adjusted via deep
learning techniques.
II. B RIEF REVIEW OF KNOWN RECOVERY ALGORITHMS
As preparation for describing the details of the proposed algorithm, several known sparse recovery algorithms are briefly
reviewed in this section. In the following, the observation vector is assumed to be y = Ax+w where A ∈ RM ×N (N ≥ M )

and x ∈ RN . Each entry of additive noise vector w ∈ RM
obeys zero mean Gaussian distribution with variance σ 2 .
A. ISTA
ISTA is a well-known sparse recovery algorithm [6] defined
by the following simple recursion: rt = st + βAT (y −
Ast ), st+1 = η(rt ; τ ), where β represents a step size and
η is the soft thresholding function defined by η(r; τ ) =
sign(r) max{|r| − τ, 0}. The parameter τ (> 0) indicates the
threshold value. After T -iterations, the estimate x̂ = sT of
the original sparse signal x is obtained. The initial value is
assumed to be s0 = 0. In order to have convergence, the
step parameter should be carefully determined [6]. Several
accelerated methods for ISTA using a momentum term have
been proposed [25], [26]. Since the proximal operator of L1 regularization term ||x||1 is the soft thresholding function,
ISTA can be seen as a proximal gradient descent algorithm [3].
B. AMP
AMP [10] is defined by the following recursion:
rt

=

y − Ast + bt rt−1 ,

(1)

=

C. OAMP
OAMP [13] is defined by the following recursive formula:
rt

=

st+1

=

vt2
τt2
Bt

st + Wt (y − Ast ),

ηdf (rt ; τt ),


||y − Ast ||22 − M σ 2
= max
,

,
trace(AT A)
1
1
=
trace(Bt BtT )vt2 + trace(Wt WtT )σ 2 ,
N
N
= I − Wt A,

III. D ETAILS OF TISTA
This section describes the details of TISTA and its training
process.
A. MMSE estimators for additive Gaussian noise channel

η(st + AT rt ; τt ),
(2)
θ
1
||st ||0 , τt = √ ||rt ||2
(3)
bt =
M
M
and it provides the final estimate x̂ = sT . It is assumed that
each entry of the sensing matrix A is generated according to
the Gaussian distribution N (0, 1/M ), i.e., Gaussian distribution with mean zero and variance 1/M . The recursive formula
of AMP seem similar to these of ISTA at a glance but there
are several critical differences. Due to the Onsager correction
term bt rt−1 in (1), the output of the linear estimator becomes
statistically decoupled and an error between each output signal
from the linear estimator and the true signal behaves as a white
Gaussian random variable in large system limit. This enable us
to use a scalar recursion called the state evolution to track the
evolution of the error variances. Another difference between
ISTA and AMP is the estimator of τt in (3), which is used as
the threshold value for the shrinkage function (2). In [10], it is
reported that AMP shows much faster convergence than that
of ISTA if the sensing matrix satisfies the above condition.
However, it is known that AMP cannot provide excellent
recovery performance for sensing matrices violating the above
condition such as non-Gaussian sensing matrices, Gaussian
matrices with large variance, Gaussian matrices with non-zero
mean, and matrices with large condition numbers [12].
st+1

for t = 0, 1, 2, . . . , T −1. To be precise, the estimator equations
on vt2 (6) and τt2 (7) (presented also in [27]) are not part of
OAMP (for example, we can use the state evolution to provide
vt2 and τt2 ) but these estimators were adopted numerical
evaluation in [13]. The matrix Wt in linear estimator (4) can be
chosen from the transpose of A, the pseudo inverse of A and
the LMMSE matrix. The nonlinear estimation unit (5) consists
of a divergence free function ηdf that replaces the Onsager
correction term. It is proved in [13] that the estimation errors
at the linear estimator (4) and non-linear estimator (5) are
statistically orthogonal if a sensing matrix is i.i.d. Gaussian
or unitary invariant. This fact provides a justification for the
state evolution of OAMP.

(4)
(5)
(6)
(7)
(8)

Let X be a real-valued random variable with probability
density function (PDF) PX (·). We assume an additive Gaussian noise channel defined by Y = X +N, where Y represents
a real-valued random variable as well. The random variable N
is a Gaussian random variable with mean 0 and variance σ 2 .
Consider the situation where a receiver can observe Y and
wish to estimate the value of X.
The MMSE estimator ηM M SE (y) is defined ηM M SE (y) =
E[X|y] where E[X|y] is the conditional expectation given by
Z ∞
E[X|y] =
xP (x|y)dx.
(9)
−∞

The posterior PDF P (x|y) is given by Bayes Theorem:
PX|Y (x|y) =

PX (x)PY |X (y|x)
,
PY (y)

where the conditional PDF is Gaussian:


1
−(y − x)2
PY |X (y|x) = √
exp
.
2σ 2
2πσ 2

(10)

(11)

In the case of the Bernoulli-Gaussian prior, PX (x) is given
by


x2
exp − 2 ,
(12)
2α
2πα2
where p represents the probability such that a non-zero element
occurs. The function δ(·) is the Dirac’s delta function. In this
case, a non-zero element follows the Gaussian PDF with mean
zero and variance α2 . The MMSE estimator for the BernoulliGaussian prior can be easily derived [29] by using Stein’s
formula:
d
ln PY (y)
(13)
ηM M SE (y; σ 2 ) = y + σ 2
dy
PX (x) = (1 − p)δ(x) + √

p

and we have
2

ηM M SE (y; σ ) =



yα2
ξ



pF (y; ξ)
,
(1 − p)F (y; σ 2 ) + pF (y; ξ)
(14)

where ξ = α2 + σ 2 and
F (z; v) = √

1
exp
2πv




2

−z
2v

.

(15)

These MMSE estimators are going to be used as a building
block of TISTA to be presented in the next subsection.
B. Recursive formula for TISTA
We assume that the sensing matrix A ∈ RM ×N is a full
rank matrix. The recursive formula of TISTA are summarized
as follows:
=

st + γt W (y − Ast ),

(16)

st+1

=

ηM M SE (rt ; τt2 ),

(17)

vt2

=

max

rt

τt2



=
+


||y − Ast ||22 − M σ 2
,

,
trace(AT A)

where the matrix W = AT (AAT )−1 is the pseudo inverse
matrix 1 of the sensing matrix A, and Z = W A. The initial
condition is s0 = 0 and the final estimate is given by x̂ = sT .
The scalar variables γt ∈ R(t = 0, 1, . . . , T − 1) are learnable
variables that are tuned in a training process. The number of
learnable variables is thus T that is much smaller than that of
LISTA [22] and LAMP [23].
An appropriate MMSE shrinkage (17) is chosen according
to the prior distribution of the original signal x. Note that the
MMSE shrinkage is also employed in [23]. The real constant
 is a sufficiently small value, e.g.,  = 10−9 . The max
operator in (18) is employed to prevent the estimate of the
variance being non-positive. The learnable variables γt in (16)
provide appropriate step sizes and control for the variance of
the MMSE shrinkage.
The true error variances τ̄t2 and v̄t2 are defined by
(20)

These error variances should be estimated as correct as possible in a sparse recovery process because the MMSE shrinkage
(17) requires to know τ̄t2 . As in the case of OAMP [13], we
make the following assumptions on the residual errors in order
to derive an error variance estimator.
The first assumption is that rt − x consists of i.i.d. zero
mean Gaussian entries. Due to this assumption, each entry of
the output from the linear estimator (16) can be seen as an
observation obtained from a virtual additive Gaussian noise
channel with the noise variance τ̄ 2 . This justifies the use of the
shrinkage function based on the MMSE estimator (17) with τ̄ 2 .
Another assumption is that st − x consists of zero mean i.i.d.
entries and satisfies E[(st −x)T AT w] = E[(st −x)T W w] = 0
for any t.
1 If

N < M , we let W = (AT A)−1 AT .

v̄t2 =

E[||y − Ast ||22 ] − M σ 2
trace(AT A)

(21)

holds.
The justification of the error variance estimator (19) for τ̄t2
is also provided by the following proposition.
Proposition 2: If each entry of st − x is i.i.d. with mean 0
and E[(st − x)T W w] = 0 is satisfied, then
τ̄t2

v̄t2
(N − 2γt trace(Z) + γt2 trace(ZZ T ))
N
γt2 σ 2
trace(W W T )
(22)
N

=

(18)

vt2
(N − 2γt trace(Z) + γt2 trace(ZZ T ))
N
γt2 σ 2
trace(W W T ),
(19)
N

τ̄t2 = E[||rt − x||22 ]/N, v̄t2 = E[||st − x||22 ]/N.

The error variance estimator for v̄t2 (18) is the same as that
of OAMP [13] and its justification comes from the following
proposition.
Proposition 1: If each entry of st − x is i.i.d. with mean 0
and E[(st − x)T AT w] = 0 is satisfied, then

+

holds.
(Proof) The residual error rt − x can be rewritten as
rt − x

= st + γt W (y − Ast ) − x
= st + γt W (Ax + w) − γt W Ast − x
(I − γt Z)(st − x) + γt W w.

=

From the definition τ̄t2 , we have
1
E[||(I − γt Z)(st − x) + γt W w||22 ]
τ̄t2 =
N
1
E[(st − x)T (I − γt Z)(I − γt Z)T (st − x)]
=
N
γt2
2γt
+
E[wT W T W w] +
E[(st − x)T (I − γt Z)T W w]
N
N
1
=
trace((I − γt Z)(I − γt Z)T )v̄t2
N
2(γt − γt2 )
γt2
trace(W W T )σ 2 +
E[(st − x)T W w].
+
N
N
The last term vanishes due to the assumption E[(st −
x)T W w] = 0 and the first term can be rewritten as
trace((I − γt W A)(I − γt W A)T )
X
X
=
(γt Zi,j )2 +
(1 − γt Zi,i )2
i

i,j:i6=j

=

γt2

X
i,j:i6=j

2
Zi,j

+

X

2
(1 − 2γt Zi,i + γt2 Zi,i
)

i

= N − 2γt trace(Z) + γt2 trace(ZZ T ).

(23)

The claim of the proposition is now obtained.
(QED)
These error variance estimators (18) and (19) play a crucial
role to give appropriate variance estimates required for the
MMSE shrinkage. Since the validity of these assumptions on
the residual errors cannot be proved, it will be experimentally
confirmed in the next section. It should be also remarked that
the TISTA recursive formula does not include neither an Onsager correction term nor a divergence free function. Thus, we
cannot expect stochastic orthogonality guaranteed in OAMP in
a process of TISTA. This means that the state evolution cannot
be applied to analyze the asymptotic performance of TISTA.

10-1

Fig. 1. The t-th iteration of the TISTA with learnable variable γt

C. Time complexity of TISTA
For treating a large scale problem, a sparse recovery algorithm should require low computational complexity for each
iteration. The time complexity required for evaluating the
recursive formula of TISTA per iteration is O(N 2 ), which
is the same time complexity as those of ISTA and AMP,
which means that TISTA has sufficient scalability for large
problems. The evaluation of the matrix-vector products, Ast
and W (y − Ast ) need O(N 2 )-time that are dominant in an
iteration. Although the evaluations of the scalar constants
trace(AT A), trace(W W T ), and trace(ZZ T ) and computation
of the pseudo inverse of A require O(N 3 )-time, they can be
pre-computed only once in advance.

error variance

10-2

10-3

10-4

10-5

0

2

4

6
8
iteration

10

12

14

Fig. 2. Comparison between the estimate τ̄ 2 and the true error variance τ 2 ;
A ∼ N (0, 1/M ), SNR = 40 dB. The optimized γt in this case are given by
γ0 = 1.67, γ1 = 4.42, γ2 = 1.41, γ3 = 1.35, γ4 = 5.66, γ5 = 1.51, γ6 =
2.83, γ7 = 1.38, γ8 = 5.84, γ9 = 0.92, γ10 = 1.13, γ11 = 1.49, γ12 =
1.70, γ13 = 1.91, γ14 = 2.18, γ15 = 2.21.

IV. P ERFORMANCE EVALUATION

D. Incremental training for TISTA

A. Details of experiments

In order to achieve reasonable recovery performance, the
trainable variables γt (t = 0, 1, . . . , T − 1) should be appropriately adjusted. By unfolding the recursive formula of
TISTA, we immediately have a signal-flow graph which is
similar to a multi-layer feedforward neural network. Figure 1
depicts a unit of the signal-flow graph corresponding to t-th
iteration of TISTA and we can stack the units to compose a
whole signal-flow graph. We here follow a standard recipe of
deep learning techniques, namely we apply mini-batch training
with a stochastic gradient descent algorithm to the signal-flow
graph of TISTA. From several experiments, we found that
the following incremental training is considerably effective
to learn appropriate values providing superior performance.
The training data consists of a number of randomly generated pairs (x, y) where y = Ax+w. The sample x follows the
prior distribution PX (x) and w is an i.i.d. Gaussian random
vector. The whole training data is divided into mini-batchs
to be used in a stochastic gradient descent algorithm such as
SGD, RMSprop or Adam.
In the t-th round of the incremental training (we call it a
generation), an optimizer tries to minimize E[||st − x||22 ] by
tuning γ0 , γ1 , . . . , γt−1 . The number of mini-batches used in
t-th generation is denoted by D (epochs). After processing D
epochs, the target of the optimizer is changed to E[||st+1 −
x||22 ]. Namely, after training the first to t-th layers, a new t + 1
layer is appended to the network and the whole network is
trained again for D epochs. Although the objective function
is changed, the values of the variables γ0 , . . . , γt−1 of the
previous generation is taken over as the initial values in the
optimization process for the new generation. In summary, the
incremental training updates the variables γt in a sequential
manner from the first layer to the last layer.

The basic conditions for computer experiments are summarized as follows. Each component of the sparse signal x is
assumed to be a realization of i.i.d. random variable following
the Bernoulli-Gaussian PDF (12) with p = 0.1, α2 = 1.
We thus use the MMSE estimator (17) for the BernoulliGaussian prior. Each component of the noise vector w follows
the zero mean Gaussian PDF with variance σ 2 . The signal
to noise ratio (SNR) of the system is defined as SN R =
E[||Ax||22 ]/E[||w||22 ]. The dimension of the sensing matrices
are set to be N = 500, M = 250. The size of the minibatch is set to 1000, and D = 200 epochs are allocated for
each generation. We used Adam optimizer with the initial
value 4.0 × 10−2 in the training phase. The experiment system
was implemented based on TensorFlow [28]. For comparison
purpose, we will include the NMSE performances of AMP
and LISTA in the following subsections. The hyper parameter
θ used in AMP is set to θ = 1.14. We used an implementation
of LISTA by Borgerding [30].
B. IID Gaussian matrix
This subsection describes the case where A ∼ N (0, 1/M ),
i.e., each component of the sensing matrix A obeys zero mean
Gaussian distribution with variance 1/M . Note that AMP is
designed for this matrix ensemble.
Figure 2 shows a comparison between the estimate τ 2 by
(19) and the empirically estimated values of the true error
variance τ̄ 2 . We can observe that the estimator τ 2 provides
accurate estimations and it justifies the use of (18) and (19)
and our assumptions on the residual errors. The caption of Fig.
2 includes a set of the optimized γt .
Figure 3 presents the average normalized MSE (NMSE)
of TISTA, LISTA and AMP as functions of iteration when

0

-5

TISTA
LISTA
AMP

-5
-10

-15
NMSE [dB]

NMSE [dB]

-15
-20
-25
-30

-20
-25
-30
-35

-35

-40

-40
2

4

6

8
10
iteration

12

14

-45

16

Fig. 3. NMSE of TISTA, LISTA and AMP; Ai,j ∼ N (0, 1/M ), SNR
= 40dB. The condition Ai,j ∼ N (0, 1/M ) is required for AMP to converge.

SNR = 40dB. The NMSE is defined by 10 log10 {||st+1 −
x||22 /||x||22 }(dB). From Fig. 3, we can observe that TISTA
provides the steepest NMSE curve than those of AMP and
LISTA at the first 16 rounds. For example, AMP and LISTA
require 16 and 10 rounds to achieve NMSE = −30dB, respectively, but TISTA needs only 6 rounds. It can be seen that
the NMSE curve of TISTA saturates around −43dB at which
TISTA converges. This means that TISTA shows significantly
faster convergence than that of AMP and LISTA in this
setting. Compared with the experimental results under the
same condition shown in [24], the NMSE values of TISTA is
almost comparable to those of LAMP [23]. CPU time required
for training a 7-layers TISTA graph was around 6 minutes by
using a PC with GPU NVIDIA GeForce GTX 1080. A large
Gaussian matrix with the size N = 5000, M = 2500 required
around 7 hours for training. It was observed that the NMSE
curve of N = 5000 is quite similar to that of N = 500.
In the next experiment, we made change on the variance of
sensing matrices to a larger number, i.e., each element in A
follows N (0, 1). Figure 4 shows the NMSE curves of TISTA
and LISTA. It should be noted that, under this condition,
AMP does not perform well, i.e, it cannot converge at all,
because the setting does not fit the required condition (Ai,j ∼
N (0, 1/M )) for achieving the guaranteed performance and the
convergence of AMP. As we can see from Fig. 4, TISTA
behaves soundly and shows faster convergence than that of
LISTA. This result suggests that TISTA is appreciably robust
against the change of the variance.
C. Binary matrix
In this subsection, we will discuss the case where sensing
matrices are binary, i.e., A ∈ {±1}M ×N . Each entry of
A is selected uniformly at random on {±1}. This situation
is closely related to multiuser detection in Coded Division
Multiple Access (CDMA) [9]. Figure 5 shows the NMSE
curves of TISTA and LISTA as a function of iteration. It can

2

4

6

8
10
iteration

12

14

16

Fig. 4. NMSE of TISTA and LISTA; Ai,j ∼ N (0, 1), SNR = 40dB. In this
case, AMP cannot converge because the variance of the matrix components
are too large.

-5

TISTA
LISTA

-10
-15
NMSE [dB]

-45

TISTA
LISTA

-10

-20
-25
-30
-35
-40
-45

2

4

6

8
10
iteration

12

14

16

Fig. 5. NMSE of TISTA and LISTA; Ai,j takes a value in {±1} uniformly
at random. SNR = 40dB. AMP is not applicable to this case.

be seen that the NMSE curves of TISTA almost coincide with
those of the Gaussian sensing matrices. This result can be
regarded as an evidence for the robustness of TISTA for nonGaussian sensing matrices.
D. Sensing matrices with large condition number
The condition number κ of a matrix is defined as the ratio of
the largest and the smallest singular values, i.e., κ = s1 /sM
where s1 ≥ s2 ≥ · · · ≥ sM are the singular values of the
matrix. It is known that regression problems regarding a matrix
with a large condition number are difficult to solve in an
accurate manner. In this subsection, we access the performance
of TISTA for sensing matrices with a large condition number.
Figure 6 presents NMSE of TISTA and AMP when SNR
= 60dB. AMP can converge up to κ = 4 but it diverges
when κ > 4. TISTA shows similar performance of AMP with
κ = 4 even when κ = 1000. This result is a strong evidence

0
-10

NMSE [dB]

-20
-30
-40
-50
-60
-70

κ=1
κ = 15
κ = 100
κ = 1000
AMP:κ = 1
AMP:κ = 4
2

4

6

8
10
iteration

12

14

16

Fig. 6. NMSE of TISTA and AMP; κ represents the condition number. SNR
= 60dB. The singular values are selected so that si /si−1 becomes a constant.
The sensing matrix A is normalized as ||A||2F = N .

that TISTA has a potential to achieve excellent NMSE even
for a sensing matrix with a large condition number.
V. C ONCLUDING S UMMARY
The crucial feature of TISTA is that it includes adjustable
variables which can be tuned by standard deep learning
techniques. The number of trainable variables of TISTA is
equal to the number of iterative rounds and it is much smaller
than those of the known learnable sparse signal recovery algorithms [22]–[24]. This feature leads to highly stable and fast
training processes of TISTA. Computer experiments indicate
that TISTA is applicable to various classes of sensing matrices
such as Gaussian matrices, binary matrices and matrices
with large condition numbers. Furthermore, numerical results
demonstrate that TISTA shows significantly faster convergence
than those of AMP and LISTA in many cases. By replacing the
MMSE shrinkage, we can expect that TISTA is also applicable
to non-sparse signal recovery problems such as detection of
BPSK signals in overloaded MIMO systems.
ACKNOWLEDGEMENT
This work is supported by JSPS Grant-in-Aid for Scientific
Research (B) Grant Number 16H02878 (TW) and Grant-inAid for Young Scientists (Start-up) Grant Number 17H06758
(ST). The last author is grateful to Keigo Takeuchi for an
inspiring seminar talk.
R EFERENCES
[1] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52,
no. 4, pp. 1289-1306, Apr. 2006.
[2] E. J. Candes and T. Tao, “Near-optimal signal recovery from random
projections: Universal encoding strategies?” IEEE Trans. Inf. Theory,
vol. 52, no. 12, pp. 5406-5425, Dec. 2006.
[3] Z. Zhang, Y. Xu, J. Yang, X. Li, and D. Zhang, “A Survey of Sparse
Representation: Algorithms and Applications,” IEEE Access, vol. 3, pp.
490-530, 2015.
[4] R. Tibshirani, “Regression shrinkage and selection via the lasso,” J.
Royal Stat. Society, Series B, vol. 58, pp. 267–288, 1996.

[5] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least Angle
Regression,” Ann. Stat., vol. 32, no. 2, pp. 407-499, Apr. 2004.
[6] A. Chambolle, R. A. DeVore, N. Lee, and B. J. Lucier, “Nonlinear
wavelet image processing: variational problems, compression, and noise
removal through wavelet shrinkage,” IEEE Trans. Image Process., vol.
7, no. 3, pp. 319–335, Mar, 1998.
[7] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding
algorithm for linear inverse problems with a sparsity constraint,” Comm.
Pure and Appl. Math., col. 57, no. 11, pp. 1413-1457, Nov. 2004.
[8] N. Parikh and S. Boyd, “Proximal Algorithms,” Foundations and Trends
in Optimization, vol. 1, no. 3, pp. 127-239, 2014.
[9] Y. Kabashima, “A CDMA multiuser detection algorithm on the basis of
belief propagation,” J. Phys. A: Math. Gen., vol. 36 11111–11121, Oct.
2003.
[10] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proceedings of the National Academy
of Sciences, vol. 106, no. 45, pp. 18914–18919, Nov. 2009.
[11] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. Motivation and construction,” IEEE
Information Theory Workshop 2010, ITW 2010, pp. 1-5, 2010.
[12] F. Caltagirone, L. Zdeborova, and F. Krzakala, “On convergence of
approximate message passing,” 2014 IEEE Int. Symp. Inf. Theory, Jun.
2014, pp. 1812-1816.
[13] J. Ma and L. Ping, “Orthogonal AMP,” IEEE Access, vol. 5, pp. 2020–
2033, 2017.
[14] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message
passing,” 2017 IEEE Int. Symp. Inf. Theory, Jun. 2017, vol. 65, no. 17,
pp. 1588-1592.
[15] K. Takeuchi, “Rigorous dynamics of expectation-propagation-based signal recovery from unitarily invariant measurements,” 2017 IEEE Int.
Symp. Inf. Theory, Jun. 2017, pp. 501-505.
[16] K. Fukushima, “Neocognitron: A self-organizing neural network model
for a mechanism of pattern recognition unaffected by shift in position,”
Bio. Cybern., vol. 36, no. 4, pp. 193-202, 1980.
[17] G. E. Hinton, R. R. Salakhutdinov, “Reducing the Dimensionality of
Data with Neural Networks,” Science, vol. 313, no. 5786, pp. 504-507,
Jun. 2006:
[18] A. Krizhevsky, I. Sutskever, G. E. Hinton, “Imagenet classification with
deep convolutional neural networks.” Advances in Neural Inf. Proc. Sys.
2012, pp. 1097-1105, 2012.
[19] G. Hinton et al., “Deep Neural Networks for Acoustic Modeling in
Speech Recognition: The Shared Views of Four Research Groups,” IEEE
Signal Processing Magazine, vol. 29, no. 6, pp. 82-97, Nov. 2012.
[20] B. Aazhang, B. P. Paris and G. C. Orsak, “Neural networks for
multiuser detection in code-division multiple-access communications,”
IEEE Trans. Comm., vol. 40, no. 7, pp. 1212-1222, jul. 1992.
[21] E. Nachmani, Y. Beéry and D. Burshtein, “Learning to decode linear
codes using deep learning,” 2016 54th Annual Allerton Conf. Comm.,
Control, and Computing, 2016, pp. 341-346.
[22] K. Gregor, and Y. LeCun, “Learning fast approximations of sparse
coding,” Proc. 27th Int. Conf. Machine Learning, pp. 399–406, 2010.
[23] M. Borgerding and P. Schniter, “Onsager-corrected deep learning for
sparse linear inverse problems,” 2016 IEEE Global Conf. Signal and
Inf. Proc. (GlobalSIP), Washington, DC, Dec. 2016, pp. 227-231.
[24] M. Borgerding, P. Schniter, and S. Rangan, “AMP-Inspired Deep
Networks for Sparse Linear Inverse Problems, ” arXiv:1612.01183v2
(2017).
[25] A. Beck, and M. Teboulle, “A fast iterative shrinkage-thresholding
algorithm for linear inverse problems,” SIAM J Imaging Sciences, vol.
2, no. 1, pp. 183-202, 2009.
[26] J. M. Bioucas-Dias and M. A. T. Figueiredo, “A New TwIST: Two-Step
Iterative Shrinkage/Thresholding Algorithms for Image Restoration,”
IEEE Trans. Image Proc., vol. 16, no. 12, pp. 2992-3004, Dec. 2007.
[27] J. Vila and P. Schniter, “Expectation-maximization gaussian-mixture
approximate message passing,” IEEE Trans. Signal Proc., vo. 61, no.
19, pp. 4658-4672, 2013.
[28] “TensorFlow: Large-scale machine learning on heterogeneous systems,”
http://tensorflow.org/ 2015. Software available from tensorflow.org.
[29] Rémi Gribonval “Should penalized least squares regression be interpreted as Maximum A Posteriori estimation?, ” IEEE Transactions on
Signal Processing, vol.59, no.5, pp. 2405-2410, May. 2011.
[30] https://github.com/mborgerding/onsager_deep_learning/blob/master/
README.md

