# Approximately Normalized Iterative Hard Thresholding for Nonlinear Compressive Sensing

Authors: Xunzhi Zhu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The nonlinear compressive sensing (NCS) is an extension of classical compressive sensing (CS) and the iterative hard thresholding (IHT) algorithm is a popular greedy-type method for solving CS. The normalized iterative hard thresholding (NIHT) is a modification of IHT and is more effective than IHT. In this paper, we propose an approximately normalized iterative hard thresholding (ANIHT) algorithm for NCS by using the approximate optimal stepsize combining with Armijo stepsize rule preiteration. Under the condition similar to restricted isometry property (RIP), we analyze the condition that can identify the iterative support sets in a finite number of iterations. Numerical experiments show the good performance of the new algorithm for the NCS.


---

## Body

## 1. Introduction

Compressed sensing (CS) [1, 2] deals with the problem of recovering sparse signals from underdetermined linear measurements. In recent years, it has attracted considerable attention in areas of signal processing, electrical engineering, computer science, and applied mathematics; see [3, 4]. However, many real-life applications in physics and biomedical sciences carry some strongly nonlinear structure, so that the linear model is not suited anymore, even as an approximation. It is of utmost interest to investigate compressive sensing with nonlinear measurements, which is called nonlinear compressive sensing (NCS) [5, 6]. So, it is necessary to consider the following NCS: to find a vector x
∈
R
n from the observations b
∈
R
m given by(1)
Φ
x
+
η
=
b
,
s
.
t
.
x
0
≤
s
,where Φ
(
x
)
=
(
ϕ
1
(
x
)
,
…
,
ϕ
m
(
x
)
)
⊤
:
R
n
→
R
m is nonlinear, η is some noise term and s
<
m
<
n, and x
0 is the l
0-norm of x, which refers to the number of nonzero elements in the vector x. The optimization problem associated with (1) is(2)
m
i
n
f
x
≜
1
2
Φ
x
-
b
2
,
s
.
t
.
x
0
≤
s
,where · is l
2-norm and r
(
x
)
≜
Φ
(
x
)
-
b is the residual function. Let Φ
(
x
) be continuously differentiable and let J
(
x
)
∈
R
m
×
n denote its Jacobian matrix. Then ∇
f
(
x
)
=
J
(
x
)
⊤
r
(
x
). Clearly, if Φ
(
x
) is linear function, model (1) and the optimization problem (2) reduce to classical CS,(3)
A
x
+
e
=
b
,
s
.
t
.
x
0
≤
s
,and problem(4)
m
i
n
1
2
A
x
-
b
2
,
s
.
t
.
x
0
≤
s
,respectively, where A
∈
R
m
×
n.Greedy methods have already proven useful and efficient to tackle (4) [7]. A variety of greedy methods have been proposed to solve (4), such as matching pursuit (MP) [8], orthogonal MP (OMP) [9], compressive sampling matching pursuit (CoSaMP) [10], subspace pursuit (SP) [11], hard thresholding pursuit (HTP) [12], and conjugate gradient iterative hard thresholding (CGIHT) [13].There is another greedy method, iterative hard thresholding (IHT) algorithm for problem (4), which was proposed by Blumensath and Davies in [14, 15]. When matrix A is row full-rank and the spectral norm A
2
<
1, IHT converges to a local minimum of (4) [14]. In [16], the authors showed that the numerical studies of IHT are not very promising and the algorithm often fails to converge when the conditions fail. Then they gave normalized IHT (NIHT) with an adaptive stepsize and line search and proved that it converges to a local minimum if A is row full-rank and s-regular, where s-regular means that any s columns of A are linear independent [17]. Cartis and Thompson [18] showed NIHT converges to local minimum if matrix A is 2
s-regular. Blumensath [5] showed that IHT can recover signals from NCS under conditions similar to those required in CS.Inspired by these works, we propose an approximately NIHT (ANIHT) algorithm to solve the NCS problem (2). Since problem (2) is in general a nonconvex programming, we can only expect to find out the stationary point rather than the local minimizer, which is different from NIHT for CS. Then we show that the accumulation point of the ANIHT algorithm is the stationary point. With the nondegeneracy and strict complementarity of stationary point, the support sets of the sequence are identified in a finite number of iterations. At last, we simulate several experiments to demonstrate the effectiveness of the algorithm.This paper is organized as follows. Section2 gives the ANIHT algorithm for (2) and proves its convergence properties. Numerical results are given in Section 3. The last section makes some concluding remarks.


## 2. Algorithm

In this section, we will present the approximately normalized iterative hard thresholding (ANIHT) algorithm for (2) and then analyze its convergence properties. Denote (5)
S
≜
x
∈
R
n
∣
x
0
≤
s
.For x
∈
R
n and set A
⊆
R
n, projector onto A is P
A
(
x
)
≜
a
r
g
m
i
n
y
∈
A
x
-
y. Note that the projection onto sparse set S, written as P
S
(
·
), sets all but s largest absolute value components of x to zero. The definition of L-stationarity was proposed in [17] based on the notion of fixed-point equation.Definition 1.
The pointx
∗
∈
S is called the stationary point of problem (2), if it holds that(6)
x
∗
∈
P
S
x
∗
-
α
∇
f
x
∗
.Note thatx
∗
∈
S is the stationary point of problem (2) if and only if [17] (7)
∇
i
f
x
∗
=
0
,
i
∈
Γ
∗
≤
1
α
M
s
x
∗
,
i
∉
Γ
∗
, where M
i
(
x
) denotes the ith largest element in absolute value of x
∈
S and Γ
∗
=
s
u
p
p
(
x
∗
)
=
{
i
∈
{
1
,
…
,
n
}
∣
x
i
∗
≠
0
}.In NIHT for CS [16], to guarantee the objective function a sufficient decrease per iteration, the authors added a stepsize strategy based on restricted isometry property (RIP) [1]. In ANIHT for NCS (2), we use the approximately optimal stepsize to accelerate the convergence and use Armijo-type stepsize and make a sufficient decrease of the objective function directly without RIP. Here (
·
)
Γ is the subvector (submatrix) obtained by discarding all but the elements (columns) in Γ. The framework of ANIHT is described as follows.Step 1.
Initializex
0, α
0
>
0, σ
>
0, 0
<
β
<
1.Step 2.
LetΓ
k
=
supp
⁡
(
x
k
) and compute (8)
x
~
k
+
1
=
P
S
x
k
-
α
k
0
J
x
k
⊤
r
x
k
,where (9)
α
k
0
=
min
⁡
α
0
,
J
x
k
⊤
r
x
k
Γ
k
2
J
x
k
Γ
k
J
x
k
⊤
r
x
k
Γ
k
2
.Step 3.
Ifsupp
⁡
(
x
~
k
+
1
)
=
Γ
k, then x
k
+
1
=
x
~
k
+
1, Γ
k
+
1
=
s
u
p
p
(
x
k
+
1
), and α
k
=
α
k
0; else compute (10)
x
k
+
1
=
P
S
x
k
-
α
k
J
x
k
⊤
r
x
kand Γ
k
+
1
=
s
u
p
p
(
x
k
+
1
), where α
k
=
α
k
0
β
m
k and m
k is the smallest positive integer m such that (11)
r
x
k
α
k
0
β
m
2
≤
r
x
k
2
-
σ
x
k
α
k
0
β
m
-
x
k
2
, and x
k
(
α
)
=
P
S
(
x
k
-
α
J
(
x
k
)
⊤
r
(
x
k
)
).Step 4.
If the stopping criterion is met, then stop.Remark 2.
We now briefly illustrate the algorithm:(i)
In Step2, a proposal point x
~
k
+
1 is calculated, whether accepting it depends on the relationship between its support set and the previous point’s; see Step 3. In addition, to accelerate the rate of convergence, we utilize the approximately optimal initial stepsize to compute x
~
k
+
1. The linear approximate of Φ
(
x
k
-
α
J
(
x
k
)
⊤
r
(
x
k
)
) at x
k is (12)
Φ
x
k
-
α
J
x
k
⊤
r
x
k
≈
Φ
x
k
-
α
J
x
k
J
x
k
⊤
r
x
k
.
Γ
k is the support of the best s term approximation to b at the current iteration. So we obtain that(13)
arg
⁡
min
⁡
Φ
x
k
-
α
J
x
k
Γ
k
J
x
k
⊤
r
x
k
Γ
k
-
b
2
,
α
>
0
=
J
x
k
⊤
r
x
k
Γ
k
2
J
x
k
Γ
k
J
x
k
⊤
r
x
k
Γ
k
2
.
This stepsize is in accordance with the optimal stepsize in NIHT for CS in [16]. Furthermore, by Assumption 4, when α is relatively small, the error introduced in the linear approximation is small; then the objective function decreases if the support set is not changed.
(ii)
Armijo-type stepsize rule in Step3 makes the choice of stepsize and support set adaptively and a sufficient decrease of the objective function meanwhile per iteration. It is well defined by Lemma 6.The following assumptions are chosen to ensure the descent property (14) of the objective function 1
/
2
r
(
x
)
2.Assumption 3.
There exists a constantU
2
s
>
0 such that J
(
x
)
(
x
-
y
)
≤
U
2
s
x
-
y for |
supp
⁡
x
∪
supp
⁡
(
y
)
|
≤
2
s.We also need the assumption that the JacobianJ
(
·
) of residual r
(
·
) is restricted Lipschitz continuous on R
n.Assumption 4.
There exists a constantJ
2
s
>
0 such that J
(
x
)
-
J
(
y
)
≤
J
2
s
x
-
y for supp
⁡
x
∪
supp
⁡
y
≤
2
s.Lemma 5.
Suppose that Assumptions3 and 4 hold. If there exists ρ
>
0 for x and the iterative point x
k of ANIHT satisfying supp
⁡
x
∪
supp
⁡
x
k
≤
2
s and x
-
x
k
≤
2
ρ, then(14)
r
x
2
≤
r
x
k
2
+
2
J
x
k
⊤
r
x
k
,
x
-
x
k
+
L
~
2
s
x
-
x
k
2
,where L
~
2
s
≜
U
2
s
2
+
J
2
s
r
x
0
+
ρ
U
2
s
J
2
s and x
0 is the initial point of ANIHT.Proof.
We first show(15)
J
x
⊤
r
x
-
J
x
k
⊤
r
x
k
,
x
-
x
k
≤
L
~
2
s
x
-
x
k
2
.Since ANIHT algorithm generates monotonically decreasing function values, then r
(
x
k
)
≤
r
(
x
0
) for all k. Direct calculation yields that (16)
J
x
⊤
r
x
-
J
x
k
⊤
r
x
k
,
x
-
x
k
=
J
x
⊤
r
x
-
J
x
⊤
r
x
k
+
J
x
⊤
r
x
k
-
J
x
k
⊤
r
x
k
,
x
-
x
k
≤
J
x
⊤
r
x
-
r
x
k
,
x
-
x
k
+
r
x
k
,
J
x
-
J
x
k
x
-
x
k
=
J
x
⊤
∫
0
1
J
x
k
+
t
x
-
x
k
d
t
x
-
x
k
,
x
-
x
k
+
r
x
k
,
J
x
-
J
x
k
x
-
x
k
≤
J
x
⊤
J
x
x
k
-
x
,
x
-
x
k
+
J
x
⊤
∫
0
1
J
x
k
+
t
x
-
x
k
d
t
-
J
x
x
-
x
k
,
x
-
x
k
+
r
x
0
J
2
s
x
-
x
k
2
=
J
x
x
-
x
k
,
J
x
x
-
x
k
+
∫
0
1
J
x
k
+
t
x
-
x
k
-
J
x
d
t
x
-
x
k
,
J
x
x
-
x
k
+
J
2
s
r
x
0
x
-
x
k
2
≤
U
2
s
2
x
-
x
k
2
+
∫
0
1
t
J
2
s
x
-
x
k
2
d
t
·
U
2
s
x
-
x
k
+
J
2
s
r
x
0
x
-
x
k
2
≤
U
2
s
2
+
J
2
s
r
x
0
+
ρ
U
2
s
J
2
s
x
-
x
k
2
=
L
~
2
s
x
-
x
k
2
. It follows from (15) that (17)
r
x
2
-
r
x
k
2
-
2
J
x
k
⊤
r
x
k
,
x
-
x
k
=
2
∫
0
1
J
x
k
+
t
x
-
x
k
⊤
r
x
k
+
t
x
-
x
k
-
J
x
k
⊤
r
x
k
,
x
-
x
k
d
t
≤
2
∫
0
1
t
L
~
2
s
x
-
x
k
2
d
t
=
L
~
2
s
x
-
x
k
2
, which completes the proof.Lemma 6.
Suppose that Assumptions3 and 4 hold and x
k is the iterative point of ANIHT with 0
<
α
0
<
1
/
L
~
2
s. Then(18)
r
x
k
α
2
≤
r
x
k
2
-
1
α
0
-
L
~
2
s
x
k
α
-
x
k
2
,
i
f
0
<
α
≤
α
0
,
r
x
k
2
-
σ
x
k
α
-
x
k
2
,
i
f
0
<
α
≤
1
L
~
2
s
+
σ
.Therefore, α
k is well defined.Proof.
According to the computation in Step2, we have (19)
x
k
α
∈
arg
⁡
min
⁡
x
-
x
k
+
α
J
x
k
⊤
r
x
k
2
,
x
0
≤
s
,which implies that (20)
x
k
α
-
x
k
+
α
J
x
k
⊤
r
x
k
2
≤
α
J
x
k
⊤
r
x
k
2
,that is,(21)
x
k
α
-
x
k
2
≤
-
2
α
J
x
k
⊤
r
x
k
,
x
k
α
-
x
k
.If supp
⁡
x
~
k
+
1
=
supp
⁡
(
x
k
), then the above inequality and the monotonicity of {
r
(
x
k
)
} yield that (22)
x
k
α
-
x
k
2
≤
-
2
α
r
x
k
,
J
x
k
x
k
α
-
x
k
≤
2
α
r
x
k
·
J
x
k
x
k
α
-
x
k
≤
2
α
0
r
x
0
U
2
s
x
k
α
-
x
k
. Otherwise, by the Armijo-type stepsize rule and the monotonicity of {
r
(
x
k
)
}, we have (23)
x
k
α
-
x
k
2
≤
1
σ
r
x
k
2
-
r
x
k
α
2
≤
1
σ
r
x
0
2
.Then x
k
(
α
)
-
x
k can be smaller than max
⁡
{
2
α
0
U
2
s
,
1
/
σ
}
·
r
(
x
0
).
By Lemma5 and (21), we get that(24)
r
x
k
α
2
≤
r
x
k
2
+
2
J
x
k
⊤
r
x
k
,
x
k
α
-
x
k
+
L
~
2
s
x
k
α
-
x
k
2
≤
r
x
k
2
-
1
α
x
k
α
-
x
k
2
+
L
~
2
s
x
k
α
-
x
k
2
=
r
x
k
2
-
1
α
-
L
~
2
s
x
k
α
-
x
k
2
. If supp
⁡
(
x
~
k
+
1
)
=
Γ
k, α
=
α
0
<
1
/
L
~
2
s, then 1
/
α
-
L
~
2
s
>
0. Otherwise, by letting 1
/
α
-
L
~
2
s
≥
σ, we can obtain the desired result by the definition of α.


### 2.1. Convergence

Combining Assumptions3 and 4 and Lemma 6, the convergence of ANIHT can be established in this subsection.Theorem 7.
Let Assumptions3 and 4 hold and let {
x
k
} be generated by ANIHT with 0
<
α
0
<
1
/
L
~
2
s. Then (i)
l
i
m
k
→
∞
x
k
+
1
-
x
k
/
α
k
=
0;
(ii)
any accumulation point of{
x
k
} is the stationary point of (2).Proof.
(i) It follows from (18) that r
(
x
k
)
2
-
r
x
k
+
1
2
≥
c
x
k
+
1
-
x
k
2, where c
=
min
⁡
{
1
/
α
0
-
L
~
2
s
,
σ
}
>
0. Then (25)
∑
k
=
0
∞
x
k
+
1
-
x
k
2
≤
1
c
∑
k
=
0
∞
r
x
k
2
-
r
x
k
+
1
2
<
1
c
r
x
0
2
<
+
∞
,which signifies l
i
m
k
→
∞
x
k
+
1
-
x
k
=
0.
It follows from Assumption3 that (26)
J
x
k
⊤
r
x
k
Γ
k
2
J
x
k
Γ
k
J
x
k
⊤
r
x
k
Γ
k
2
≥
1
U
2
s
2
,and then α
k
0
≥
min
⁡
{
α
0
,
1
/
U
2
s
2
}. By Lemma 6 and the definition of α
k in the algorithm, we have (27)
α
k
≥
min
⁡
α
k
0
,
β
L
~
2
s
+
σ
.Therefore, α
k is bounded from below by a positive constant min
⁡
α
0
,
β
/
L
~
2
s
+
σ (β
/
L
~
2
s
+
σ
<
1
/
U
2
s
2). We can conclude that (28)
l
i
m
k
→
∞
x
k
+
1
-
x
k
α
k
=
0
.
(ii) Suppose thatx
∗ is an accumulation point of the sequence {
x
k
}; then there exists a subsequence {
x
k
j
} converging to x
∗ and l
i
m
j
→
∞
x
k
j
+
1
=
x
∗ by (i). For the simplification of notation, denote α
k as α by its boundedness and Γ
∗
=
supp
⁡
(
x
∗
). Based on (29)
x
k
j
+
1
=
P
S
x
k
j
-
α
J
x
k
j
⊤
r
x
k
jin Step 2, we consider two cases.
Case  1 (i
∈
Γ
∗). The convergence of {
x
k
j
} and {
x
k
j
+
1
} guarantees that, for some n
1
>
0, x
i
k
j
>
0, x
i
k
j
+
1
>
0 for all j
>
n
1. The definition of the projection on S shows that (30)
x
i
k
j
+
1
=
x
i
k
j
-
α
J
x
k
j
⊤
r
x
k
j
i
. Taking j
→
∞, we have (
J
x
∗
⊤
r
x
∗
)
i
=
0.
Case  2 (i
∉
Γ
∗). If there exists an n
2
>
0 such that, for all j
>
n
2, x
i
k
j
+
1
=
0, the projection implies that (31)
x
i
k
j
-
α
J
x
k
j
⊤
r
x
k
j
i
≤
M
s
x
k
j
+
1
.Letting j
→
∞ and exploiting the continuity of the function M
s, we obtain that (32)
α
J
x
∗
⊤
r
x
∗
i
≤
M
s
x
∗
. On the other hand, if there exists an infinite number of indices of k
j for x
i
k
j
+
1
>
0, as the same proof in Case  1, it follows that (
J
(
x
∗
)
⊤
r
(
x
∗
)
)
i
=
0. Since α
k is bounded from below by a positive constant, we have(33)
J
x
∗
⊤
r
x
∗
i
=
0
,
i
f
i
∈
Γ
∗
,
≤
1
α
M
s
x
∗
,
i
f
i
∉
Γ
∗
, which means x
∗ is a stationary point of (2).We are now ready to show that under suitable conditions the support set of a point is identified in a finite number of iterations. We can easily verify that ifx
∗
0
=
s, then the support set of x in a sufficiently small neighborhood of x
∗ is identified. For x
∗
0
<
s, we introduce the concept of strict complementarity to identify the support set.Definition 8.
The pointx
∗
∈
S is called nondegenerate if x
∗
0
=
s. The condition(34)
J
x
∗
⊤
r
x
∗
Γ
∗
=
0
,
i
∈
Γ
∗
≠
0
,
i
∉
Γ
∗is called strict complementarity condition of (2), where Γ
∗
=
supp
⁡
(
x
∗
).Theorem 9.
For any sequence{
z
k
} converging to z
∗, we have the following:(i)
ifz
∗ is nondegenerate point, then(35)
supp
⁡
z
k
=
supp
⁡
z
∗for all k sufficiently large and(36)
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
k
=
0
,where Γ
z
k
=
supp
⁡
(
z
k
) and Γ
z
∗
=
supp
⁡
(
z
∗
).(ii)
ifz
∗ satisfies strict complementarity, then (35) holds if and only if (36) holds.Proof.
(i) Suppose thatz
∗ is a nondegenerate point. We have (35) when z
k is sufficiently close to z
∗. By (35) and the continuity of J
(
·
)
⊤
r
(
·
), we can easily get that (37)
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
k
=
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
∗
=
J
z
∗
⊤
r
z
∗
Γ
z
∗
=
0
.
(ii) Suppose thatz
∗ satisfies strict complementarity condition. The “only if” part can be obtained by the similar proof as in (i). Let (36) hold. Since z
k
→
z
∗, Γ
z
∗
⊆
Γ
z
k. Assume that there is an infinite subsequence K and an index j such that j
∉
Γ
z
∗ and j
∈
Γ
z
k for all k
∈
K. By (36), (38)
J
z
∗
⊤
r
z
∗
j
=
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
j
=
0
. While the strict complementarity condition (34) implies that [
J
(
z
∗
)
⊤
r
(
z
∗
)
]
j
≠
0. This contradiction proves that Γ
z
∗
=
Γ
z
k.


## 2.1. Convergence

Combining Assumptions3 and 4 and Lemma 6, the convergence of ANIHT can be established in this subsection.Theorem 7.
Let Assumptions3 and 4 hold and let {
x
k
} be generated by ANIHT with 0
<
α
0
<
1
/
L
~
2
s. Then (i)
l
i
m
k
→
∞
x
k
+
1
-
x
k
/
α
k
=
0;
(ii)
any accumulation point of{
x
k
} is the stationary point of (2).Proof.
(i) It follows from (18) that r
(
x
k
)
2
-
r
x
k
+
1
2
≥
c
x
k
+
1
-
x
k
2, where c
=
min
⁡
{
1
/
α
0
-
L
~
2
s
,
σ
}
>
0. Then (25)
∑
k
=
0
∞
x
k
+
1
-
x
k
2
≤
1
c
∑
k
=
0
∞
r
x
k
2
-
r
x
k
+
1
2
<
1
c
r
x
0
2
<
+
∞
,which signifies l
i
m
k
→
∞
x
k
+
1
-
x
k
=
0.
It follows from Assumption3 that (26)
J
x
k
⊤
r
x
k
Γ
k
2
J
x
k
Γ
k
J
x
k
⊤
r
x
k
Γ
k
2
≥
1
U
2
s
2
,and then α
k
0
≥
min
⁡
{
α
0
,
1
/
U
2
s
2
}. By Lemma 6 and the definition of α
k in the algorithm, we have (27)
α
k
≥
min
⁡
α
k
0
,
β
L
~
2
s
+
σ
.Therefore, α
k is bounded from below by a positive constant min
⁡
α
0
,
β
/
L
~
2
s
+
σ (β
/
L
~
2
s
+
σ
<
1
/
U
2
s
2). We can conclude that (28)
l
i
m
k
→
∞
x
k
+
1
-
x
k
α
k
=
0
.
(ii) Suppose thatx
∗ is an accumulation point of the sequence {
x
k
}; then there exists a subsequence {
x
k
j
} converging to x
∗ and l
i
m
j
→
∞
x
k
j
+
1
=
x
∗ by (i). For the simplification of notation, denote α
k as α by its boundedness and Γ
∗
=
supp
⁡
(
x
∗
). Based on (29)
x
k
j
+
1
=
P
S
x
k
j
-
α
J
x
k
j
⊤
r
x
k
jin Step 2, we consider two cases.
Case  1 (i
∈
Γ
∗). The convergence of {
x
k
j
} and {
x
k
j
+
1
} guarantees that, for some n
1
>
0, x
i
k
j
>
0, x
i
k
j
+
1
>
0 for all j
>
n
1. The definition of the projection on S shows that (30)
x
i
k
j
+
1
=
x
i
k
j
-
α
J
x
k
j
⊤
r
x
k
j
i
. Taking j
→
∞, we have (
J
x
∗
⊤
r
x
∗
)
i
=
0.
Case  2 (i
∉
Γ
∗). If there exists an n
2
>
0 such that, for all j
>
n
2, x
i
k
j
+
1
=
0, the projection implies that (31)
x
i
k
j
-
α
J
x
k
j
⊤
r
x
k
j
i
≤
M
s
x
k
j
+
1
.Letting j
→
∞ and exploiting the continuity of the function M
s, we obtain that (32)
α
J
x
∗
⊤
r
x
∗
i
≤
M
s
x
∗
. On the other hand, if there exists an infinite number of indices of k
j for x
i
k
j
+
1
>
0, as the same proof in Case  1, it follows that (
J
(
x
∗
)
⊤
r
(
x
∗
)
)
i
=
0. Since α
k is bounded from below by a positive constant, we have(33)
J
x
∗
⊤
r
x
∗
i
=
0
,
i
f
i
∈
Γ
∗
,
≤
1
α
M
s
x
∗
,
i
f
i
∉
Γ
∗
, which means x
∗ is a stationary point of (2).We are now ready to show that under suitable conditions the support set of a point is identified in a finite number of iterations. We can easily verify that ifx
∗
0
=
s, then the support set of x in a sufficiently small neighborhood of x
∗ is identified. For x
∗
0
<
s, we introduce the concept of strict complementarity to identify the support set.Definition 8.
The pointx
∗
∈
S is called nondegenerate if x
∗
0
=
s. The condition(34)
J
x
∗
⊤
r
x
∗
Γ
∗
=
0
,
i
∈
Γ
∗
≠
0
,
i
∉
Γ
∗is called strict complementarity condition of (2), where Γ
∗
=
supp
⁡
(
x
∗
).Theorem 9.
For any sequence{
z
k
} converging to z
∗, we have the following:(i)
ifz
∗ is nondegenerate point, then(35)
supp
⁡
z
k
=
supp
⁡
z
∗for all k sufficiently large and(36)
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
k
=
0
,where Γ
z
k
=
supp
⁡
(
z
k
) and Γ
z
∗
=
supp
⁡
(
z
∗
).(ii)
ifz
∗ satisfies strict complementarity, then (35) holds if and only if (36) holds.Proof.
(i) Suppose thatz
∗ is a nondegenerate point. We have (35) when z
k is sufficiently close to z
∗. By (35) and the continuity of J
(
·
)
⊤
r
(
·
), we can easily get that (37)
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
k
=
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
Γ
z
∗
=
J
z
∗
⊤
r
z
∗
Γ
z
∗
=
0
.
(ii) Suppose thatz
∗ satisfies strict complementarity condition. The “only if” part can be obtained by the similar proof as in (i). Let (36) hold. Since z
k
→
z
∗, Γ
z
∗
⊆
Γ
z
k. Assume that there is an infinite subsequence K and an index j such that j
∉
Γ
z
∗ and j
∈
Γ
z
k for all k
∈
K. By (36), (38)
J
z
∗
⊤
r
z
∗
j
=
l
i
m
k
→
∞
J
z
k
⊤
r
z
k
j
=
0
. While the strict complementarity condition (34) implies that [
J
(
z
∗
)
⊤
r
(
z
∗
)
]
j
≠
0. This contradiction proves that Γ
z
∗
=
Γ
z
k.


## 3. Numerical Experiments

In this part, sensor localization problem and phase retrieval problem will be stimulated. In both examples, the stop criteria will be set as(
J
x
k
⊤
r
(
x
k
)
)
Γ
k
≤
ϵ
2, where ϵ
2 is pretty small in different cases or the maximum iterative times being equal to 5000.Sensor localization problem can be described as follows: given M known anchors p
1
,
p
2
,
…
,
p
M
∈
R
n, the purpose is to find a sensor x
∈
R
n satisfying (39)
x
-
p
i
2
+
η
i
=
b
i
,
i
=
1
,
…
,
M
,
x
0
≤
s
, where η
i, i
=
1
,
…
,
M, is the noise (which obeys the normal distribution with zero expectation and σ
0
2 variance here). The problem of finding an x
∈
R
n satisfying above equalities is the same as finding an optimal solution to the optimization problem (4) with f
(
x
)
=
∑
i
=
1
M
(
x
-
p
i
2
-
b
i
)
2. We first compare the A
N
I
H
T on an example with M
=
80, N
=
120, and s
=
1,2
,
…
,
10. Each component of the M vectors p
1
,
…
,
p
M was randomly and independently generated from a standard normal distribution. Then the true vector x
o
r
i
g and b are designed as following MATLAB codes:(40)
x
o
r
i
g
=
z
e
r
o
s
n
,
1
;
T
=
r
a
n
d
p
e
r
m
n
;
x
o
r
i
g
T
1
:
s
=
10
∗
r
a
n
d
s
,
1
;
b
i
=
x
o
r
i
g
-
p
i
2
+
σ
0
*
r
a
n
d
n
1,1
.For each value of s
(
s
=
1,2
,
…
,
10
), we ran A
N
I
H
T algorithm from 100 different and randomly generated initial data sets. The numbers of runs of 100 in which the methods found the “correct” solution are given in Table 1. Here, the “correct” solution x
(
x
0
≤
s
) means that f
(
x
)
≤
f
(
x
o
r
i
g
), or x and x
o
r
i
g are very close, say (41)
x
-
x
o
r
i
g
x
<
10
-
3
. As can be clearly seen by the results in the table, the A
N
I
H
T performs well in terms of the success probability. For more details, when the “true” solutions are quite sparse compared to the dimension n, A
N
I
H
T can almost recover all the “ture” solutions, while the performance is becoming worse as s rises.Table 1
The success numbers over 100 runs under different noise.
s
1
2
3
4
5
6
7
8
9
10
σ
0
=
0.00
100
99
100
99
94
90
86
82
82
75
σ
0
=
0.0
1
100
100
99
96
99
93
91
91
82
75
σ
0
=
0.1
0
100
98
94
83
81
68
66
56
54
46Then we runA
N
I
H
T algorithm in a higher dimensional data set, where N
=
2
M, s
=
⌈
0.01
N
⌉, and M
=
100,200
,
…
,
2000. For each data set, we run 40 times and record the average results (in which the unsuccessful recoveries are expelled). Figure 1 shows the performance of A
N
I
H
T when addressing this problem.Figure 1
Average results over40 simulations with different noise.Phase retrieval is to recover a signal from the magnitude of its Fourier transform, or of any other linear transform. Due to the loss of Fourier phase information, the problem is generally ill-posed. The phase retrieval problem can be described as follows: given M known measurement vectors p
1
,
p
2
,
…
,
p
M
∈
R
n, the purpose is to reconstruct a signal x
∈
R
n satisfying (42)
x
,
p
i
2
+
η
i
=
b
i
,
i
=
1
,
…
,
M
,
x
0
≤
s
, where p
i is the ith column of the general matrix or the discrete Fourier transform (DFT) and η
i, i
=
1
,
…
,
M, is the noise (which obeys the normal distribution with zero expectation and σ
0
2 variance here). Also this problem is equivalent to recover an optimal solution to the optimization problem (4) with f
(
x
)
=
∑
i
=
1
M
(
x
⊤
A
i
x
-
b
i
)
2, where A
i
=
p
i
p
i
⊤, i
=
1
,
…
,
M.There are some other methods for sparse phase retrieval and the codes are available. So we can compare ourA
N
I
H
T algorithm with them. We first compare the A
N
I
H
T algorithm with the partial sparse-simplex method (P
S
S) and greedy sparse-simplex (G
S
S) method in [17] with M
=
80, N
=
120, and s
=
2,3
,
…
,
9 which is identical to those in [17]. The true vector x
o
r
i
g and the measurement vectors p
1
,
p
2
,
…
,
p
M are generated as that produced in sensor localization problem; b is designed as following MATLAB codes:(43)
b
i
=
x
o
r
i
g
,
p
i
2
+
σ
0
∗
r
a
n
d
n
1,1
.For each value of s
(
s
=
2,3
,
…
,
9
), we ran A
N
I
H
T algorithm from 100 different and randomly generated initial data sets. The numbers of runs of 100 in which the methods found the “correct” solution are given in Table 2. As can be clearly seen by the results in the table, the A
N
I
H
T outperforms P
S
S and G
S
S in terms of the success probability. What is more, the data in the row with σ
0
=
0 are higher than those with σ
0
=
0.01 and σ
0
=
0.1.Table 2
The success numbers of four methods over 100 runs under different noise.
s
2
3
4
5
6
7
8
9
A
NIHT
σ
0
=
0.00
100
98
91
85
76
70
52
41
σ
0
=
0.01
92
98
93
88
68
60
47
34
σ
0
=
0.10
78
85
82
85
68
54
43
31
PSS
σ
0
=
0.00
35
28
22
9
5
9
4
3
σ
0
=
0.01
33
25
20
7
8
7
5
3
σ
0
=
0.10
30
24
16
4
5
6
3
2
GSS
σ
0
=
0.00
80
75
69
20
17
13
7
6
σ
0
=
0.01
57
52
35
30
27
15
6
5
σ
0
=
0.10
65
54
27
60
25
12
5
3We also compare ourA
N
I
H
T algorithm with G
E
S
P
A
R in [19] to recover a signal from the magnitude of its Fourier transform. Namely, it is to find a real-valued discrete time signal x
∈
R
N from its magnitude-squared value of an N point discrete Fourier transform (DFT): (44)
b
j
=
∑
k
=
1
n
x
k
e
-
2
π
i
j
-
1
k
-
1
/
N
2
,
j
=
1
,
…
,
N
. We denote by F the DFT matrix; then elements F
j
k
=
e
-
2
π
i
(
j
-
1
)
(
k
-
1
)
/
N and b
=
F
x
2, where ·
2 denotes the element-wise absolute-squared value. We get b by the pseudo MATLAB codes: (45)
b
=
a
b
s
f
f
t
x
o
r
i
g
.
∧
2
+
σ
0
∗
r
a
n
d
n
N
,
1
.To see the accuracy of the solutions and the speed of these two methods, we run the two methods forn increasing from 512 to 3072 and keeping N
=
2
n, s
=
1
%
n. We also test them under noiseless and two noise levels, σ
0
=
0.01 and σ
0
=
0.1. From Table 3, we can see that A
N
I
H
T outperforms G
E
S
P
A
R in terms of both average CPU time and average relative error for large n
(
n
≥
2048
).Table 3
Average results withN
=
2
n
,
s
=
1
%
n.
n
σ
0
=
0
Time
Relative error
A
NIHT
GESPAR
ANIHT
GESPAR
256
0.1869
0.09
3.013
e
-
09
4.97
e
-
06
512
0.541
0.18
2.816
e
-
11
1.936
e
-
06
1024
1.021
1.03
9.234
e
-
11
2.235
e
-
06
2048
2.732
21.83
1.587
e
-
06
4.683
e
-
06
3072
5.443
491.22
1.067
e
-
06
1.300
e
-
06
n
σ
0
=
0.01
Time
Relative error
ANIHT
GESPAR
ANIHT
GESPAR
256
0.183
0.0960
1.756
e
-
06
1.6470
e
-
05
512
0.456
0.2022
5.404
e
-
06
1.0901
e
-
06
1024
0.976
1.0755
1.663
e
-
06
1.0517
e
-
06
2048
2.701
22.64
1.184
e
-
06
1.0579
e
-
06
3072
5.251
472.48
3.780
e
-
07
1.3010
e
-
06
n
σ
0
=
0.1
Time
Relative error
ANIHT
GESPAR
ANIHT
GESPAR
256
0.179
0.0914
5.912
e
-
06
4.0600
e
-
05
512
0.167
0.2056
7.974
e
-
05
1.0901
e
-
06
1024
1.013
0.7589
3.178
e
-
05
1.0517
e
-
06
2048
2.750
37.36
9.063
e
-
06
9.9184
e
-
06
3072
5.417
299.62
4.403
e
-
06
1.3010
e
-
06


## 4. Conclusion

Nonlinear CS (NCS) not only is of academic interest but also might be important in many real-world applications when the measurements cannot be designed to be perfectly linear. In this paper, we have proposed an ANIHT algorithm for NCS and studied its convergence. We have showed that any accumulation point of the algorithm is the stationary point. The support set of the sequence can be identified with the assumption of nondegeneracy and strict complementarity of stationary point. The numerical experiments show that ANIHT algorithm is effective for NCS. In the future, we will further consider other methods for nonlinear least square problem to improve the rate of convergence, such as L-M method or cubic regularization methods [20].



---

*Source: 2594752-2016-08-15.xml*