# On a Gradient-Based Algorithm for Sparse Signal Reconstruction in the Signal/Measurements Domain

Authors: Ljubiša Stanković; Miloš Daković
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Sparse signals can be recovered from a reduced set of samples by using compressive sensing algorithms. In common compressive sensing methods the signal is recovered in the sparsity domain. A method for the reconstruction of sparse signals which reconstructs the missing/unavailable samples/measurements is recently proposed. This method can be efficiently used in signal processing applications where a complete set of signal samples exists. The missing samples are considered as the minimization variables, while the available samples are fixed. Reconstruction of the unavailable signal samples/measurements is preformed using a gradient-based algorithm in the time domain, with an adaptive step. Performance of this algorithm with respect to the step-size and convergence are analyzed and a criterion for the step-size adaptation is proposed in this paper. The step adaptation is based on the gradient direction angles. Illustrative examples and statistical study are presented. Computational efficiency of this algorithm is compared with other two commonly used gradient algorithms that reconstruct signal in the sparsity domain. Uniqueness of the recovered signal is checked using a recently introduced theorem. The algorithm application to the reconstruction of highly corrupted images is presented as well.


---

## Body

## 1. Introduction

A signal is sparse in a transformation domain if the number of nonzero coefficients is much lower than the number of signal samples. For linear signal transforms the signal samples can be considered as linear combinations (measurements) of the signal transform coefficients. Sparse signals can be fully recovered from a reduced set of samples/measurements. Compressive sensing theory is dealing with sparse signal reconstruction [1–22]. In common signal processing problems, for signals sparse in a transformation domain, signal samples can be considered as observations since they are linear combinations of the transformation domain coefficients. The goal of compressive sensing is to reduce the set of signal measurements while preserving complete information about the whole signal. In some applications unavailable signal samples result from the system physical constraints including unavailability to perform all signal measurements. Signals may also contain arbitrarily positioned samples that are heavily corrupted. In these cases it is better to remove such samples and consider them as unavailable in the signal recovery [6]. Compressive sensing theory based reconstruction algorithms may be also used when the unavailable signal samples are not a result of an intentional strategy to compress the data. Note that one of the initial compressive sensing theory successes in applications (computed tomography reconstruction) was not related to the intentional compressive sensing strategy but to the physical problem constraints, restricting the set and positions of available data.Various algorithms for reconstruction of sparse signals from a reduced set of observations are introduced [8–17]. Two of the most important approaches are based either on the gradient of the sparsity measure [8] or on the orthogonal matching pursuit methods [16].The topic of this paper is the reconstruction of signals with some arbitrary positioned unavailable samples. These samples may result from a compressive sensing strategy or from their unavailability due to various reasons. A method for the unavailable signal samples reconstruction, considering them as variables, has been proposed in [22]. In contrast to the other common methods that recover sparse signals in their sparsity domain this method reconstructs missing samples/measurements to make the set of samples/measurements complete. The available samples are fixed (invariant) during the reconstruction process. An analysis of the algorithm performance is presented in this paper. A relation between reconstructed signal bias and the algorithm step-size is derived and statistically confirmed. A new criterion for the algorithm step-size adaptation is proposed. The improved version of this algorithm is compared with two other common gradient-based algorithms proving its accuracy and computational efficiency. The discrete Fourier transform domain is used as a case study, although the algorithm application is not restricted to this transform [23]. The solution uniqueness is checked after signal reconstruction by using the recently proposed theorem [24].The paper is organized as follows. After the definitions in the next section, the adaptive gradient algorithm is presented. In the comments to the algorithm the step-size and bias are considered and illustrated. A new criterion for the algorithm parameter adaptation is proposed in Section3 as well. Reconstruction is illustrated on examples and uniqueness of the obtained solution is analyzed. Statistical performance are checked in Section 4. An illustration of the algorithm application to the image reconstruction is given as well.


## 2. Gradient-Based Reconstruction of Missing Samples

Discrete-time signalx
(
n
) with N samples will be considered (with vector notation x
=
x
0
x
1
⋯
x
N
-
1
T). Its transform coefficients will be denoted by X
(
k
)
=
T
[
x
(
n
)
]. If the number of nonzero coefficients X
(
k
)
≠
0 for k
∈
{
k
1
,
k
2
,
…
,
k
s
} is such that s
≪
N then the signal x
(
n
) is sparse in this transformation domain. Next assume that only M
<
N signal samples are available in the time domain at the discrete-time positions:(1)
n
i
∈
N
A
=
n
1
,
n
2
,
…
,
n
M
⊂
N
=
0,1
,
2
,
…
,
N
-
1
.Under some conditions, studied within the compressive sensing theory [1–3], full signal reconstruction is possible based on the reduced set of available samples if the signal is sparse in a transformation domain. In the signal theory this problem can be formulated as a full recovery of the unavailable samples as well. In general, a reconstruction algorithm is based on the sparsity measure minimization. Nonzero transform coefficients’ counting is the simplest sparsity measure. Counting can be achieved by the mathematical form X
0, referred to as the “l
0-norm” [1, 20]. The signal reconstruction problem statement is then(2)
m
i
n
X
0
s
u
b
j
e
c
t
t
o
y
=
A
X
.Vector y
=
x
n
1
x
n
2
⋯
x
(
n
M
)
T contains the available signal samples, while the transform vector X
=
X
1
X
2
⋯
X
N
T elements are the unknown transform coefficients. Matrix A is the measurement matrix. In common signal transforms it corresponds to the inverse transformation matrix where the rows corresponding to the unavailable samples are omitted. If x
=
Ψ
X then the matrix A is a submatrix of Ψ containing its rows corresponding to the discrete-time instants (1).Although there are some approaches to reconstruct the signal using thel
0-norm based formulation [25] in principle, this is an NP-hard combinatorial optimization problem. It is the reason why the closest convex l
1-norm of the signal transform is used instead of the “l
0-norm” in the minimization process:(3)
m
i
n
X
1
s
u
b
j
e
c
t
t
o
y
=
A
X
,where X
1
=
∑
k
=
0
N
-
1
‍
X
(
k
). Under the conditions defined using the restricted isometry property (RIP) [3, 4], minimization of the l
1-norm can produce the same result as (2) [1, 22, 26]. Common gradient-based algorithms reformulate the problem defined by (3) into a form suitable for application of standard minimization tools. Minimization is done in the sparsity domain. One such a method is briefly reviewed in Appendix.


### 2.1. Algorithm

A gradient-based algorithm that minimizes the sparsity measure by varying the missing sample values is presented next. In this gradient-based minimization approach the missing samples are considered as variables [22]. The missing sample values are varied in an iterative way to produce lower sparsity measure values and finally to approach the minimum of the convex l
1-norm based sparsity measure (3) with an acceptable accuracy. Note that if the reconstruction conditions are met for the l
1-norm [3] then the l
1-norm minimum position will correspond to the true values of the missing samples.The signal(4)
x
a
0
n
=
0
for  missing  samples
,
n
∈
N
Q
x
n
for  available  samples
,
n
∈
N
Awill be used as the initial estimate. Note that this signal would follow as a result of the l
2-norm based minimization of the signal transform. The set of missing sample positions is N
Q. Values of the available samples are fixed and will be considered as constants. The missing sample values are varied through iterations. If we denote by x
a
(
m
) the vector of signal values reconstructed after m iterations, then the goal is to achieve(5)
m
i
n
X
a
1
s
u
b
j
e
c
t
t
o
x
a
m
n
=
x
n
for
n
∈
N
A
, where the components of X
a are X
a
(
k
)
=
T
[
x
a
(
m
)
(
n
)
]. In order to find the position of function z
=
X
a
1 minimum, the relation for iterative missing samples (variables) calculation is defined by using the gradient-based descend on the sparsity measure X
a
1 as(6)
y
c
m
+
1
=
y
c
m
-
α
∂
X
a
1
N
∂
y
c
y
c
=
y
c
m
.Vector y
c
(
m
) contains the minimization variables (missing signal sample values) in the mth iteration. The gradient vector coordinates g
(
n
i
)
=
∂
X
a
1
/
∂
y
c
N can be estimated, in the mth iteration, by using the finite differences of the sparsity measure. They are calculated for each missing sample at instants n
i
∈
N
Q:(7)
g
n
i
=
X
a
+
1
-
X
a
-
1
2
Δ
N
, where (8)
X
a
+
k
=
T
x
a
+
n
X
a
-
k
=
T
x
a
-
n
,
(9)
x
a
+
n
=
x
a
m
n
+
Δ
δ
n
-
n
i
x
a
-
n
=
x
a
m
n
-
Δ
δ
n
-
n
i
.Again for n
i
∈
N
A the available signal values are not changed, g
(
n
i
)
=
0. All gradient vector g
(
m
) values, including the gradient zero and nonzero coordinates, in the mth iteration are denoted by g
(
m
)
(
n
).Through statistical study [22] it has been concluded that an appropriate step parameter value in (6) is related to the finite difference step as α
=
2
Δ.


## 2.1. Algorithm

A gradient-based algorithm that minimizes the sparsity measure by varying the missing sample values is presented next. In this gradient-based minimization approach the missing samples are considered as variables [22]. The missing sample values are varied in an iterative way to produce lower sparsity measure values and finally to approach the minimum of the convex l
1-norm based sparsity measure (3) with an acceptable accuracy. Note that if the reconstruction conditions are met for the l
1-norm [3] then the l
1-norm minimum position will correspond to the true values of the missing samples.The signal(4)
x
a
0
n
=
0
for  missing  samples
,
n
∈
N
Q
x
n
for  available  samples
,
n
∈
N
Awill be used as the initial estimate. Note that this signal would follow as a result of the l
2-norm based minimization of the signal transform. The set of missing sample positions is N
Q. Values of the available samples are fixed and will be considered as constants. The missing sample values are varied through iterations. If we denote by x
a
(
m
) the vector of signal values reconstructed after m iterations, then the goal is to achieve(5)
m
i
n
X
a
1
s
u
b
j
e
c
t
t
o
x
a
m
n
=
x
n
for
n
∈
N
A
, where the components of X
a are X
a
(
k
)
=
T
[
x
a
(
m
)
(
n
)
]. In order to find the position of function z
=
X
a
1 minimum, the relation for iterative missing samples (variables) calculation is defined by using the gradient-based descend on the sparsity measure X
a
1 as(6)
y
c
m
+
1
=
y
c
m
-
α
∂
X
a
1
N
∂
y
c
y
c
=
y
c
m
.Vector y
c
(
m
) contains the minimization variables (missing signal sample values) in the mth iteration. The gradient vector coordinates g
(
n
i
)
=
∂
X
a
1
/
∂
y
c
N can be estimated, in the mth iteration, by using the finite differences of the sparsity measure. They are calculated for each missing sample at instants n
i
∈
N
Q:(7)
g
n
i
=
X
a
+
1
-
X
a
-
1
2
Δ
N
, where (8)
X
a
+
k
=
T
x
a
+
n
X
a
-
k
=
T
x
a
-
n
,
(9)
x
a
+
n
=
x
a
m
n
+
Δ
δ
n
-
n
i
x
a
-
n
=
x
a
m
n
-
Δ
δ
n
-
n
i
.Again for n
i
∈
N
A the available signal values are not changed, g
(
n
i
)
=
0. All gradient vector g
(
m
) values, including the gradient zero and nonzero coordinates, in the mth iteration are denoted by g
(
m
)
(
n
).Through statistical study [22] it has been concluded that an appropriate step parameter value in (6) is related to the finite difference step as α
=
2
Δ.


## 3. Comments on the Algorithm

(i)
The algorithm inputs are the signal lengthN, the set of available samples N
A, and the available signal values x
(
n
i
), n
i
∈
N
A.
(ii)
For the DFT instead of using signals defined by (9) and their transforms for each n
i
∈
N
Q it is possible to calculate (10)
X
a
+
k
=
X
a
m
k
+
Δ
D
n
i
k
X
a
-
k
=
X
a
m
k
-
Δ
D
n
i
kwith X
a
(
m
)
(
k
)
=
D
F
T
[
x
a
(
m
)
(
n
)
] and D
n
i
(
k
)
=
D
F
T
[
δ
(
n
-
n
i
)
]
=
e
x
p
(
-
j
2
π
n
i
k
/
N
). Note that D
n
i
(
k
) values are not dependent on the signal and the iteration number m. They can be calculated only once. Similar simplification can be done for any linear signal transform T
[
·
].


### 3.1. Adaptive Step

The influence of step-sizeΔ to the solution (minimum position) precision is analyzed next. Assume that we have obtained the exact solution for the missing samples. The change of sparsity measure is tested by the change of signal sample x
(
n
0
) for ±
Δ. For a signal of sparsity s in the DFT domain, whose form is x
(
n
)
=
∑
i
=
1
s
‍
A
i
e
j
2
π
n
0
k
i
/
N, the signals with the steps used for the gradient estimate calculation are x
a
+
(
n
)
=
x
(
n
)
+
Δ
δ
(
n
-
n
0
) and x
a
-
(
n
)
=
x
(
n
)
-
Δ
δ
(
n
-
n
0
). Their sparsity measures are (11)
X
a
+
1
=
∑
i
=
1
s
A
i
+
Δ
e
-
j
2
π
n
0
k
i
/
N
+
N
-
s
Δ
X
a
-
1
=
∑
i
=
1
s
A
i
-
Δ
e
-
j
2
π
n
0
k
i
/
N
+
N
-
s
Δ
. In the worst case, amplitudes A
i are in phase with e
-
j
2
π
n
0
k
i
/
N and Δ
≤
A
i when (12)
X
a
+
1
=
∑
i
=
1
s
A
i
+
s
Δ
+
N
-
s
Δ
=
X
1
+
N
Δ
X
a
-
1
=
∑
i
=
1
s
A
i
-
s
Δ
+
N
-
s
Δ
=
X
1
+
N
-
2
s
Δ
.It means that g
(
n
0
)
=
X
a
+
1
-
X
a
-
1
/
(
2
N
Δ
)
=
s
/
N
≠
0. This is an important fact leading to the conclusion that the stationary state of the iterative algorithm is biased. The algorithm moves the solution x
(
n
0
) to a biased value x
(
n
0
)
-
b in order to produce g
(
n
0
)
=
0 in the algorithm stationary point. For the zero value of the gradient we have x
a
+
(
n
0
)
=
x
(
n
0
)
-
b
+
Δ and x
a
-
(
n
0
)
=
x
(
n
)
-
b
-
Δ. By replacing Δ with Δ
-
b in X
a
+
1 and Δ with Δ
+
b in X
a
-
1 the worst case sparsity measures are X
a
+
1
=
X
1
+
N
(
Δ
-
b
) and X
a
-
1
=
X
1
+
(
N
-
2
s
)
(
Δ
+
b
). In this case the bias b value can be obtained from X
a
+
1
=
X
a
-
1 as(13)
N
Δ
-
b
=
N
-
2
s
Δ
+
b
b
=
s
N
-
s
Δ
≅
s
N
Δ
for
s
≪
N
.The bias limit is proportional to the step-size Δ. Obviously, the bias limit can be reduced by using small values of the step Δ. Calculation with a very small step Δ would be inefficient and time consuming with an extremely large number of iterations. An efficient way to use the gradient-based algorithm is in adapting its step-size. In the initial iteration the step-size of an order of signal amplitude is used:(14)
Δ
=
m
a
x
n
∈
N
A
x
n
.When the algorithm reaches a stationary point with this step Δ then the bias will dominate and the mean squared error will be almost constant. The algorithm can improve the solution precision by reducing the step-size.This simple relation between the bias and step-size has been confirmed through statistical analysis on more complex cases, with a large number of missing samples.


### 3.2. Adaptation Criterion

The gradient-based algorithm convergence improvement by changing steps is analyzed in detail in [22]. A criterion that can efficiently detect the event that the algorithm has reached the stationary point with regard to the mean square error (the vicinity of the sparsity measure minimum defined by the bias) described in the previous subsection is proposed next. The presented criterion is based on the direction change of the gradient vector. When the vicinity of the minimum sparsity measure point is reached within the region defined by the bias value, then the gradient estimate oscillates and changes its direction for almost 180 degrees. The angle between two successive gradient vector estimations g
(
m
-
1
) and g
(
m
) denoted as β
m can be calculated as(15)
β
m
=
arccos
⁡
∑
n
=
0
N
-
1
g
m
-
1
n
g
m
n
∑
n
=
0
N
-
1
g
m
-
1
n
2
∑
n
=
0
N
-
1
g
m
n
2
. If this angle β
m is, for example, above 17
0
∘ then we can conclude that the signal missing samples (variables) reached oscillatory nature around the minimal sparsity measure position. When these values of the angles are detected the algorithm step Δ is reduced, for example, for Δ
/
10
→
Δ or Δ
/
10
→
Δ. The iterations are continued starting from the reconstructed signal values reached in the previous step.When the minimum of the sparsity measure is reached with a sufficiently smallΔ, the value of Δ can also be used as an indicator of the solution precision. For example, if the signal amplitudes are of order 1 and s
/
N
=
0.1 taking Δ
=
1 in the first iteration will produce the solution with a precision better than 20 [dB]. Then if the step Δ is reduced to Δ
=
0.1 a precision better than 40 [dB] will be obtained and so on. Value of Δ can be used as the algorithm stopping criterion.


### 3.3. Stopping Criterion

The precision of the result in iterative algorithms is commonly estimated based on the change of the result values in the last iterations. Therefore an average of changes in a large number of variables is a good estimate of the achieved precision. The value of(16)
T
r
=
10
log
10
⁡
∑
n
∈
N
Q
x
a
p
n
-
x
a
m
n
2
∑
n
∈
N
Q
x
a
m
n
2can be used as a rough estimate of the reconstruction error to signal ratio. In this relation x
a
p
(
n
) is the signal reconstructed before Δ reduction and x
a
(
m
)
(
n
) is the reconstructed signal after iterations done with the reduced step. Value of T
r can be used as the algorithm stopping criterion. If the value of T
r is above the required precision T
m
a
x (e.g., if T
r
>
-
100 dB), the calculation procedure should be continued with reduced values of step Δ.


### 3.4. Large Step Influence

A possible divergence of a gradient-based algorithm is related to large stepsα
=
2
Δ. We can write (17)
X
a
+
k
-
X
a
-
k
=
X
a
m
k
+
Δ
D
n
i
k
-
X
a
m
k
-
Δ
D
n
i
k
=
Δ
D
n
i
k
1
+
X
a
m
k
Δ
D
n
i
k
-
1
-
X
a
m
k
Δ
D
n
i
k
. In general for a complex number a
=
X
a
(
m
)
(
k
)
/
(
Δ
D
n
i
(
k
)
), with |
a
|
≪
1 and a large Δ, from the problem geometry we can show that the following bounds hold 0
≤
1
+
a
-
1
-
a
≤
2
a. Therefore,(18)
0
≤
X
a
+
k
-
X
a
-
k
≤
2
X
a
m
k
. Lower limit 0 is obtained if a is imaginary-valued, while the upper limit 2
X
a
(
m
)
(
k
) follows if a is real-valued. If Δ is large the missing signal values will be adapted for a value independent on Δ. The missing sample values will oscillate within the range of the original signal values until the step Δ is reduced in iterations. The variable values will be arbitrary changed within the signal amplitude order as far as Δ is too large. It will not influence further convergence of the algorithm, when the step Δ assumes appropriately small values.A pseudocode of this algorithm is presented in Algorithm1.Algorithm 1:Reconstruction.
Require:
(i) Set of missing/omitted sample positionsN
Q
(ii) Available samplesx
(
n
), n
∉
N
Q
(1)    Set x
a
(
0
)
(
n
)
←
x
(
n
)                   ⊳ for n
∉
N
Q
(2)    Set x
a
(
0
)
(
n
)
←
0                      ⊳ for n
∈
N
Q
(3)    Set m
←
0
(4)    Set Δ
←
max
⁡
x
a
(
0
)
(
n
)
(5)    repeat
(6)     Set x
a
p
(
n
)
=
x
a
(
m
)
(
n
)                     ⊳  for each n
(7)     repeat
(8)      m
←
m
+
1
(9)      for  n
i
←
0 to N
-
1  do
(10)    if  n
i
∈
N
Q  then
(11)     X
a
+
(
k
)
←
D
F
T
x
a
(
m
)
(
n
)
+
Δ
δ
(
n
-
n
i
)
(12)     X
a
-
(
k
)
←
D
F
T
x
a
(
m
)
(
n
)
-
Δ
δ
(
n
-
n
i
)
(13)     g
m
(
n
i
)
←
(
1
/
N
)
∑
k
=
0
N
-
1
X
a
+
(
k
)
-
X
a
-
(
k
)
(14)    else
(15)     g
m
(
n
i
)
←
0
(16)    end if
(17)    x
a
m
+
1
(
n
i
)
←
x
a
m
(
n
i
)
-
g
m
(
n
i
)
(18)   end for
(19)   β
m
=
arccos
⁡
(
∑
n
=
0
N
-
1
g
m
-
1
(
n
)
g
m
(
n
)
/
∑
n
=
0
N
-
1
(
g
(
m
-
1
)
(
n
)
)
2
∑
n
=
0
N
-
1
(
g
(
m
)
(
n
)
)
2
)
(20)  until  β
m
<
17
0
∘
(21)  Δ
←
Δ
/
10
(22)  T
r
=
10
log
10
⁡
(
∑
n
∈
N
Q
x
a
p
n
-
x
a
m
n
2
/
∑
n
∈
N
Q
x
a
m
n
2
)
(23) until  T
r
<
T
m
a
x
(24) return  x
a
(
m
)
(
n
)
Output:
(i) Reconstructed signalx
R
(
n
)
=
x
a
(
m
)
(
n
)Illustrative Example
1. The effects analyzed within the comments on the algorithm will be illustrated on a simple signal:(19)
x
n
=
3
cos
⁡
2
π
n
N
-
π
8
-
0.2
sin
⁡
4
π
n
N
+
π
2 with N
=
8 and sparsity s
=
4. Assume that the samples at positions n
∈
N
Q
=
{
1,6
} are missing. Reconstruction of the signal is done in 60 iterations. The initial value of the step Δ
=
1 and the initial values of missing samples x
(
1
)
=
0 and x
(
6
)
=
0 are used. The missing sample values in the first 20 iterations are given in Figure 1(a) by dots (connected by a line). After 6 iterations the algorithm does not significantly change the missing sample values with given step-size Δ
=
1 (lower subplot within the figure presents changes of variables in a highly zoomed scale). When the stationary point vicinity is reached with Δ
=
1 the gradient estimate is almost zero-valued. Its direction change for almost 180 degrees. The measure values are on the contour with approximately the same value (circles on contour plot). The algorithm resumes its fast approach toward the exact missing sample values after the algorithm step is reduced to Δ
=
0.1 in the 20
t
h iteration. When a new stationary state is reached change of Δ to Δ
=
0.01 is done and the approach to the correct missing sample values is continued.Figure 1
(a) Illustration of a signal reconstruction using variable step gradient algorithm. (b) Original signal. (c) Signal with two missing samples. (d) Reconstructed signal. (e) Reconstructed signal MSE as a function of the iteration number.
(a)
(b)
(c)
(d)
(e)The bias upper limit in the stationary state for stepΔ
=
1, according to (13), is s
/
N
-
s
Δ
=
1. After each reduction of the step-size from Δ to Δ
/
10 the bias caused MSE will be reduced for almost 20 [dB]. The result for the signal reconstruction and the obtained MSE for the calculated missing values x
(
1
) and x
(
6
) are presented in Figure 1.Illustrative Example
2. The reconstruction process is repeated using the signal(20)
x
n
=
3
sin
⁡
20
π
n
N
+
cos
⁡
60
π
n
N
+
0.7
sin
⁡
46
π
n
N and N
=
32. Its sparsity is s
=
6. Missing sample positions are at n
∈
N
Q
=
{
0,1
,
2,5
,
7,8
,
11,15,18,27,28,30
}. The reconstruction results for this signal are presented in Figure 2. After the signal is reconstructed the problem is to establish if this reconstruction is unique in the sense that there is no other signal (missing sample values) with the same or lower sparsity. A very simple check of the solution uniqueness, using the missing sample positions and the positions of the reconstructed signal nonzero components in the sparsity domain, is given in [24] by the following.Figure 2
Gradient-based reconstruction of a sparse signal with reconstructed MSE (a). Angle between successive gradient estimationsβ
m and the MSE as a function of the number of iterations with step changes based on the gradient angles (b).
(a)
(b)Theorem 1.
Consider a signalx
(
n
) that is sparse in the DFT domain with unknown sparsity. Assume that the signal length is N
=
2
r samples and that Q samples are missing at the instants q
m
∈
N
Q. Also assume that the reconstruction is performed and that the DFT of reconstructed signal is of sparsity s. Assume that the positions of the reconstructed nonzero values in the DFT are k
0
i
∈
K
s
=
{
k
01
,
k
02
,
…
,
k
0
s
}
. Reconstruction result is unique if the inequality(21)
s
<
N
-
m
a
x
h
=
0,1
,
…
,
r
-
1
2
h
Q
2
h
-
1
-
s
+
2
S
2
r
-
hholds. Integers Q
2
h and S
2
r
-
h are calculated as (22)
Q
2
h
=
m
a
x
b
=
0,1
,
…
,
2
h
-
1
card
⁡
q
:
q
∈
N
Q
,
mod
⁡
q
,
2
h
=
b
S
2
r
-
h
=
∑
l
=
1
Q
2
h
-
1
P
h
l
P
h
l
=
s
o
r
t
b
=
0,1
,
…
,
2
r
-
h
-
1
card
⁡
k
:
k
∈
K
s
,
mod
⁡
k
,
2
r
-
h
=
b
,where P
h
(
1
)
≤
P
h
(
2
)
≤
⋯
≤
P
h
(
2
r
-
h
).
Signal independent uniqueness corresponds to the worst case signal form, whenS
2
r
-
h
=
0.The answer is obtained almost immediately, since the computational complexity of the Theorem is of orderO
(
N
). The proof is given in [24].For the considered set of missing samplesn
∈
N
Q
=
{
0,1
,
2,5
,
7,8
,
11,15,18,27,28,30
}, from the previous example, we easily get (23)
Q
2
0
=
8
,
Q
2
1
=
m
a
x
6,6
=
6
,
Q
2
2
=
m
a
x
3,2
,
3,4
=
4
,
Q
2
3
=
m
a
x
2,1
,
2,2
,
1,1
,
1,2
=
4
,
Q
2
4
=
2
. For signal independent uniqueness, (24)
s
<
32
-
max
⁡
1
×
8
-
1
,
2
×
2
-
1
,
4
×
4
-
1
,
8
×
2
-
1
,
16
×
2
-
1
-
s
=
16
-
s
. It means that s
<
8, with considered N
Q, guarantees uniqueness of the reconstructed signal for any signal parameters. When a sparse signal is reconstructed, for some sets of K
s uniqueness may be guaranteed for higher values of sparsity s. Additional details and examples may be found in [24].


## 3.1. Adaptive Step

The influence of step-sizeΔ to the solution (minimum position) precision is analyzed next. Assume that we have obtained the exact solution for the missing samples. The change of sparsity measure is tested by the change of signal sample x
(
n
0
) for ±
Δ. For a signal of sparsity s in the DFT domain, whose form is x
(
n
)
=
∑
i
=
1
s
‍
A
i
e
j
2
π
n
0
k
i
/
N, the signals with the steps used for the gradient estimate calculation are x
a
+
(
n
)
=
x
(
n
)
+
Δ
δ
(
n
-
n
0
) and x
a
-
(
n
)
=
x
(
n
)
-
Δ
δ
(
n
-
n
0
). Their sparsity measures are (11)
X
a
+
1
=
∑
i
=
1
s
A
i
+
Δ
e
-
j
2
π
n
0
k
i
/
N
+
N
-
s
Δ
X
a
-
1
=
∑
i
=
1
s
A
i
-
Δ
e
-
j
2
π
n
0
k
i
/
N
+
N
-
s
Δ
. In the worst case, amplitudes A
i are in phase with e
-
j
2
π
n
0
k
i
/
N and Δ
≤
A
i when (12)
X
a
+
1
=
∑
i
=
1
s
A
i
+
s
Δ
+
N
-
s
Δ
=
X
1
+
N
Δ
X
a
-
1
=
∑
i
=
1
s
A
i
-
s
Δ
+
N
-
s
Δ
=
X
1
+
N
-
2
s
Δ
.It means that g
(
n
0
)
=
X
a
+
1
-
X
a
-
1
/
(
2
N
Δ
)
=
s
/
N
≠
0. This is an important fact leading to the conclusion that the stationary state of the iterative algorithm is biased. The algorithm moves the solution x
(
n
0
) to a biased value x
(
n
0
)
-
b in order to produce g
(
n
0
)
=
0 in the algorithm stationary point. For the zero value of the gradient we have x
a
+
(
n
0
)
=
x
(
n
0
)
-
b
+
Δ and x
a
-
(
n
0
)
=
x
(
n
)
-
b
-
Δ. By replacing Δ with Δ
-
b in X
a
+
1 and Δ with Δ
+
b in X
a
-
1 the worst case sparsity measures are X
a
+
1
=
X
1
+
N
(
Δ
-
b
) and X
a
-
1
=
X
1
+
(
N
-
2
s
)
(
Δ
+
b
). In this case the bias b value can be obtained from X
a
+
1
=
X
a
-
1 as(13)
N
Δ
-
b
=
N
-
2
s
Δ
+
b
b
=
s
N
-
s
Δ
≅
s
N
Δ
for
s
≪
N
.The bias limit is proportional to the step-size Δ. Obviously, the bias limit can be reduced by using small values of the step Δ. Calculation with a very small step Δ would be inefficient and time consuming with an extremely large number of iterations. An efficient way to use the gradient-based algorithm is in adapting its step-size. In the initial iteration the step-size of an order of signal amplitude is used:(14)
Δ
=
m
a
x
n
∈
N
A
x
n
.When the algorithm reaches a stationary point with this step Δ then the bias will dominate and the mean squared error will be almost constant. The algorithm can improve the solution precision by reducing the step-size.This simple relation between the bias and step-size has been confirmed through statistical analysis on more complex cases, with a large number of missing samples.


## 3.2. Adaptation Criterion

The gradient-based algorithm convergence improvement by changing steps is analyzed in detail in [22]. A criterion that can efficiently detect the event that the algorithm has reached the stationary point with regard to the mean square error (the vicinity of the sparsity measure minimum defined by the bias) described in the previous subsection is proposed next. The presented criterion is based on the direction change of the gradient vector. When the vicinity of the minimum sparsity measure point is reached within the region defined by the bias value, then the gradient estimate oscillates and changes its direction for almost 180 degrees. The angle between two successive gradient vector estimations g
(
m
-
1
) and g
(
m
) denoted as β
m can be calculated as(15)
β
m
=
arccos
⁡
∑
n
=
0
N
-
1
g
m
-
1
n
g
m
n
∑
n
=
0
N
-
1
g
m
-
1
n
2
∑
n
=
0
N
-
1
g
m
n
2
. If this angle β
m is, for example, above 17
0
∘ then we can conclude that the signal missing samples (variables) reached oscillatory nature around the minimal sparsity measure position. When these values of the angles are detected the algorithm step Δ is reduced, for example, for Δ
/
10
→
Δ or Δ
/
10
→
Δ. The iterations are continued starting from the reconstructed signal values reached in the previous step.When the minimum of the sparsity measure is reached with a sufficiently smallΔ, the value of Δ can also be used as an indicator of the solution precision. For example, if the signal amplitudes are of order 1 and s
/
N
=
0.1 taking Δ
=
1 in the first iteration will produce the solution with a precision better than 20 [dB]. Then if the step Δ is reduced to Δ
=
0.1 a precision better than 40 [dB] will be obtained and so on. Value of Δ can be used as the algorithm stopping criterion.


## 3.3. Stopping Criterion

The precision of the result in iterative algorithms is commonly estimated based on the change of the result values in the last iterations. Therefore an average of changes in a large number of variables is a good estimate of the achieved precision. The value of(16)
T
r
=
10
log
10
⁡
∑
n
∈
N
Q
x
a
p
n
-
x
a
m
n
2
∑
n
∈
N
Q
x
a
m
n
2can be used as a rough estimate of the reconstruction error to signal ratio. In this relation x
a
p
(
n
) is the signal reconstructed before Δ reduction and x
a
(
m
)
(
n
) is the reconstructed signal after iterations done with the reduced step. Value of T
r can be used as the algorithm stopping criterion. If the value of T
r is above the required precision T
m
a
x (e.g., if T
r
>
-
100 dB), the calculation procedure should be continued with reduced values of step Δ.


## 3.4. Large Step Influence

A possible divergence of a gradient-based algorithm is related to large stepsα
=
2
Δ. We can write (17)
X
a
+
k
-
X
a
-
k
=
X
a
m
k
+
Δ
D
n
i
k
-
X
a
m
k
-
Δ
D
n
i
k
=
Δ
D
n
i
k
1
+
X
a
m
k
Δ
D
n
i
k
-
1
-
X
a
m
k
Δ
D
n
i
k
. In general for a complex number a
=
X
a
(
m
)
(
k
)
/
(
Δ
D
n
i
(
k
)
), with |
a
|
≪
1 and a large Δ, from the problem geometry we can show that the following bounds hold 0
≤
1
+
a
-
1
-
a
≤
2
a. Therefore,(18)
0
≤
X
a
+
k
-
X
a
-
k
≤
2
X
a
m
k
. Lower limit 0 is obtained if a is imaginary-valued, while the upper limit 2
X
a
(
m
)
(
k
) follows if a is real-valued. If Δ is large the missing signal values will be adapted for a value independent on Δ. The missing sample values will oscillate within the range of the original signal values until the step Δ is reduced in iterations. The variable values will be arbitrary changed within the signal amplitude order as far as Δ is too large. It will not influence further convergence of the algorithm, when the step Δ assumes appropriately small values.A pseudocode of this algorithm is presented in Algorithm1.Algorithm 1:Reconstruction.
Require:
(i) Set of missing/omitted sample positionsN
Q
(ii) Available samplesx
(
n
), n
∉
N
Q
(1)    Set x
a
(
0
)
(
n
)
←
x
(
n
)                   ⊳ for n
∉
N
Q
(2)    Set x
a
(
0
)
(
n
)
←
0                      ⊳ for n
∈
N
Q
(3)    Set m
←
0
(4)    Set Δ
←
max
⁡
x
a
(
0
)
(
n
)
(5)    repeat
(6)     Set x
a
p
(
n
)
=
x
a
(
m
)
(
n
)                     ⊳  for each n
(7)     repeat
(8)      m
←
m
+
1
(9)      for  n
i
←
0 to N
-
1  do
(10)    if  n
i
∈
N
Q  then
(11)     X
a
+
(
k
)
←
D
F
T
x
a
(
m
)
(
n
)
+
Δ
δ
(
n
-
n
i
)
(12)     X
a
-
(
k
)
←
D
F
T
x
a
(
m
)
(
n
)
-
Δ
δ
(
n
-
n
i
)
(13)     g
m
(
n
i
)
←
(
1
/
N
)
∑
k
=
0
N
-
1
X
a
+
(
k
)
-
X
a
-
(
k
)
(14)    else
(15)     g
m
(
n
i
)
←
0
(16)    end if
(17)    x
a
m
+
1
(
n
i
)
←
x
a
m
(
n
i
)
-
g
m
(
n
i
)
(18)   end for
(19)   β
m
=
arccos
⁡
(
∑
n
=
0
N
-
1
g
m
-
1
(
n
)
g
m
(
n
)
/
∑
n
=
0
N
-
1
(
g
(
m
-
1
)
(
n
)
)
2
∑
n
=
0
N
-
1
(
g
(
m
)
(
n
)
)
2
)
(20)  until  β
m
<
17
0
∘
(21)  Δ
←
Δ
/
10
(22)  T
r
=
10
log
10
⁡
(
∑
n
∈
N
Q
x
a
p
n
-
x
a
m
n
2
/
∑
n
∈
N
Q
x
a
m
n
2
)
(23) until  T
r
<
T
m
a
x
(24) return  x
a
(
m
)
(
n
)
Output:
(i) Reconstructed signalx
R
(
n
)
=
x
a
(
m
)
(
n
)Illustrative Example
1. The effects analyzed within the comments on the algorithm will be illustrated on a simple signal:(19)
x
n
=
3
cos
⁡
2
π
n
N
-
π
8
-
0.2
sin
⁡
4
π
n
N
+
π
2 with N
=
8 and sparsity s
=
4. Assume that the samples at positions n
∈
N
Q
=
{
1,6
} are missing. Reconstruction of the signal is done in 60 iterations. The initial value of the step Δ
=
1 and the initial values of missing samples x
(
1
)
=
0 and x
(
6
)
=
0 are used. The missing sample values in the first 20 iterations are given in Figure 1(a) by dots (connected by a line). After 6 iterations the algorithm does not significantly change the missing sample values with given step-size Δ
=
1 (lower subplot within the figure presents changes of variables in a highly zoomed scale). When the stationary point vicinity is reached with Δ
=
1 the gradient estimate is almost zero-valued. Its direction change for almost 180 degrees. The measure values are on the contour with approximately the same value (circles on contour plot). The algorithm resumes its fast approach toward the exact missing sample values after the algorithm step is reduced to Δ
=
0.1 in the 20
t
h iteration. When a new stationary state is reached change of Δ to Δ
=
0.01 is done and the approach to the correct missing sample values is continued.Figure 1
(a) Illustration of a signal reconstruction using variable step gradient algorithm. (b) Original signal. (c) Signal with two missing samples. (d) Reconstructed signal. (e) Reconstructed signal MSE as a function of the iteration number.
(a)
(b)
(c)
(d)
(e)The bias upper limit in the stationary state for stepΔ
=
1, according to (13), is s
/
N
-
s
Δ
=
1. After each reduction of the step-size from Δ to Δ
/
10 the bias caused MSE will be reduced for almost 20 [dB]. The result for the signal reconstruction and the obtained MSE for the calculated missing values x
(
1
) and x
(
6
) are presented in Figure 1.Illustrative Example
2. The reconstruction process is repeated using the signal(20)
x
n
=
3
sin
⁡
20
π
n
N
+
cos
⁡
60
π
n
N
+
0.7
sin
⁡
46
π
n
N and N
=
32. Its sparsity is s
=
6. Missing sample positions are at n
∈
N
Q
=
{
0,1
,
2,5
,
7,8
,
11,15,18,27,28,30
}. The reconstruction results for this signal are presented in Figure 2. After the signal is reconstructed the problem is to establish if this reconstruction is unique in the sense that there is no other signal (missing sample values) with the same or lower sparsity. A very simple check of the solution uniqueness, using the missing sample positions and the positions of the reconstructed signal nonzero components in the sparsity domain, is given in [24] by the following.Figure 2
Gradient-based reconstruction of a sparse signal with reconstructed MSE (a). Angle between successive gradient estimationsβ
m and the MSE as a function of the number of iterations with step changes based on the gradient angles (b).
(a)
(b)Theorem 1.
Consider a signalx
(
n
) that is sparse in the DFT domain with unknown sparsity. Assume that the signal length is N
=
2
r samples and that Q samples are missing at the instants q
m
∈
N
Q. Also assume that the reconstruction is performed and that the DFT of reconstructed signal is of sparsity s. Assume that the positions of the reconstructed nonzero values in the DFT are k
0
i
∈
K
s
=
{
k
01
,
k
02
,
…
,
k
0
s
}
. Reconstruction result is unique if the inequality(21)
s
<
N
-
m
a
x
h
=
0,1
,
…
,
r
-
1
2
h
Q
2
h
-
1
-
s
+
2
S
2
r
-
hholds. Integers Q
2
h and S
2
r
-
h are calculated as (22)
Q
2
h
=
m
a
x
b
=
0,1
,
…
,
2
h
-
1
card
⁡
q
:
q
∈
N
Q
,
mod
⁡
q
,
2
h
=
b
S
2
r
-
h
=
∑
l
=
1
Q
2
h
-
1
P
h
l
P
h
l
=
s
o
r
t
b
=
0,1
,
…
,
2
r
-
h
-
1
card
⁡
k
:
k
∈
K
s
,
mod
⁡
k
,
2
r
-
h
=
b
,where P
h
(
1
)
≤
P
h
(
2
)
≤
⋯
≤
P
h
(
2
r
-
h
).
Signal independent uniqueness corresponds to the worst case signal form, whenS
2
r
-
h
=
0.The answer is obtained almost immediately, since the computational complexity of the Theorem is of orderO
(
N
). The proof is given in [24].For the considered set of missing samplesn
∈
N
Q
=
{
0,1
,
2,5
,
7,8
,
11,15,18,27,28,30
}, from the previous example, we easily get (23)
Q
2
0
=
8
,
Q
2
1
=
m
a
x
6,6
=
6
,
Q
2
2
=
m
a
x
3,2
,
3,4
=
4
,
Q
2
3
=
m
a
x
2,1
,
2,2
,
1,1
,
1,2
=
4
,
Q
2
4
=
2
. For signal independent uniqueness, (24)
s
<
32
-
max
⁡
1
×
8
-
1
,
2
×
2
-
1
,
4
×
4
-
1
,
8
×
2
-
1
,
16
×
2
-
1
-
s
=
16
-
s
. It means that s
<
8, with considered N
Q, guarantees uniqueness of the reconstructed signal for any signal parameters. When a sparse signal is reconstructed, for some sets of K
s uniqueness may be guaranteed for higher values of sparsity s. Additional details and examples may be found in [24].


## 4. Statistical Analysis and Comparison

Consider a general real-valued form of a signal sparse in the DFT domain(25)
x
n
=
∑
i
=
1
K
A
i
cos
⁡
2
π
n
k
i
N
+
φ
i
.The total number of signal samples is N. The signal sparsity s
=
2
K is varied from s
=
2 to s
=
N
/
2. The signal amplitudes A
i, frequencies k
i, and phases φ
i are random variables. Signal amplitudes A
i are Gaussian random variables with unity variance; the frequency indices k
i are random numbers within 1
≤
k
i
≤
N
-
1. The phases are uniform random variables 0
≤
φ
i
≤
2
π, in each signal realization. Signal is reconstructed in 100 different realizations for each sparsity s and random positions of missing samples. The statistical calculations are performed for N
=
128. The reconstructed signals are denoted by x
R
(
n
). The reconstruction results are presented in Figure 3(a). The signal-to-reconstruction-error ratio (SRR) is given in [dB](26)
S
R
R
=
10
log
10
⁡
∑
n
=
0
N
-
1
x
n
2
∑
n
=
0
N
-
1
x
n
-
x
R
n
2
.Figure 3
Signal-to-reconstruction-error (SRR) obtained by using Algorithm1, averaged over 100 realizations for various values of sparsity s and number of available samples M. (a) The total number of samples is N
=
128. (b) With a Gaussian noise in the input signal, S
N
R
=
20 [dB] and N
=
64
.
(a)
(b)Bright colors indicate the region where the algorithm had fully recovered missing samples (compared to the original samples) in all realizations, while dark colors indicate the region where the algorithm could not recover missing samples in any realization. In the transition region forM slightly greater than 2
s we have both the cases when the signal recovery is not achieved and the cases of full signal recovery. A stopping criterion for the accuracy of 120 [dB] is used. It corresponds to a precision in the recovered signal of the same order as in input samples, if they are acquired by a 20-bit A/D converter. The case with N
=
64 is repeated with an additive input Gaussian noise such that the input signal-to-noise ratio is 20 [dB] in each realization (Figure 3(b)). The reconstruction error in this case is limited by the input signal-to-noise ratio.The average reconstruction error in the noise-free cases is related to the number of the full recovery events. ForN
=
64 the number of the full recovery events is checked and presented in Figures 4(a) and 4(b). The average number of the algorithm iterations to produce the required precision, as a function of the number of missing samples and signal sparsity s, is presented as well (Figure 4(c)) along with the corresponding average computation time (in seconds) for the Windows PC with Intel Dual Core processor (Figure 4(d)). The average computation time is proportional to the average number of iterations multiplied by the number of missing samples (variables) Q
=
N
-
M.Figure 4
((a), (b)) The percentage of the full recovery events as a function of the number of available samplesM and the sparsity s in the case of N
=
64. (c) The average number of iterations as a function of the number of missing samples and sparsity. (d) The average computation time.
(a)
(b)
(c)
(d)The efficiency of the presented algorithm is compared with the standard routines used for thel
1-norm based minimization problems. The performance of the proposed algorithm are compared with the algorithm that recasts the recovery problem (3) into regression framework (LASSO-ISTA) and a linear programming framework with the primal-dual interior point method (L1-magic code in MATLAB). A simple review of the LASSO-ISTA algorithm is presented in Appendix. All there algorithms are run using 100 sparse signals with random parameters. The algorithm parameters are set in such a way that the average computational time for all cases is about 25ms. The results are presented in Table 1. Columns’ notation in the table is s for sparsity and Q
=
N
-
M for the number of missing samples and MAE stands for the mean absolute error. Calculation time using MATLAB is presented in all cases.Table 1
MAE and computational time for the L1-magic (LP-DP), LASSO-ISTA, and the proposed algorithm.
s
N
-
M
L1-magic
LASSO
Proposed
MAE
Time
MAE
Time
MAE
Time
8
8
7.94
×
1
0
-
11
29.7
1.24
×
1
0
-
4
25.9
6.05
×
1
0
-
12
7.6
16
8
6.70
×
1
0
-
10
28.4
1.72
×
1
0
-
4
25.9
7.56
×
1
0
-
12
7.9
32
8
1.28
×
1
0
-
9
29.8
2.25
×
1
0
-
4
31.5
4.96
×
1
0
-
12
13.1
8
16
1.64
×
1
0
-
10
31.2
1.30
×
1
0
-
4
30.1
6.92
×
1
0
-
12
15.2
16
16
1.22
×
1
0
-
5
25.4
1.91
×
1
0
-
4
30.6
6.11
×
1
0
-
12
18.9
32
16
5.44
×
1
0
-
5
22.7
7.16
×
1
0
-
4
26.4
1.07
×
1
0
-
11
21.5
8
24
2.83
×
1
0
-
10
23.2
1.37
×
1
0
-
4
25.7
5.25
×
1
0
-
12
20.6
16
24
4.62
×
1
0
-
6
21.6
2.51
×
1
0
-
4
25.4
1.01
×
1
0
-
11
22.8
32
24
2.58
×
1
0
-
4
24.2
1.17
×
1
0
-
3
25.1
1.88
×
1
0
-
11
35.8
8
32
6.66
×
1
0
-
10
20.4
1.51
×
1
0
-
4
25.4
7.83
×
1
0
-
12
28.1
16
32
1.94
×
1
0
-
4
22.4
4.09
×
1
0
-
4
25.0
1.89
×
1
0
-
11
32.8
32
32
1.08
×
1
0
-
4
25.0
8.08
×
1
0
-
4
25.1
1.95
×
1
0
-
11
74.1We can conclude that in most of the considered cases the proposed algorithm outperforms other two algorithms in both the calculation time and accuracy.Example
3. The algorithm can be applied to other signal transforms as well [23]. An example with missing (salt-and-pepper noise) pixels in an image with the 2D-DCT as its sparsity domain will be presented. Consider a standard MATLAB test image canoe.tif (Figure 5) and assume that 50% of the samples are corrupted (noisy image). These pixels with extreme values are removed and considered as unavailable. Since the coefficients in the 2D-DCT converge fast we will assume that the sparsity condition is satisfied (without imposing it in an explicit way by making some 2D-DCT coefficients equal zero). The image is reconstructed by minimizing the 2D-DCT sparsity measure and varying missing pixel values according to the presented algorithm. The original, corrupted, and reconstructed images after 4th and 50th iteration are shown in Figure 5.Figure 5
Original MATLAB image canoe.tif (a). Corrupted image with salt-and-pepper noise in 50% of the pixels (b). Reconstructed images after 4 and 50 iterations (c).
(a)
(b)
(c)


## 5. Conclusion

Performance of the recently proposed gradient-based signal reconstruction algorithm in the domain of the signal samples/measurements is studied. Application of this algorithm is very efficient in signal processing problems with reduced set of data, since it is easy to define a set of missing signal samples/measurements. The values of signal samples are considered as variables, while the available samples are fixed and not changed in the algorithm. An efficient criterion, based on the gradient directions, is proposed for the algorithm step-size adaptation. Statistical analysis confirms calculation efficiency. The algorithm can easily be adapted to the cases on nonuniform signal samples, when only a reduced set of linear combinations of the signal samples are available.



---

*Source: 6212674-2016-06-19.xml*