# A General Solution to Least Squares Problems with Box Constraints and Its Applications

Authors: Yueyang Teng; Shouliang Qi; Dayu Xiao; Lisheng Xu; Jianhua Li; Yan Kang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The main contribution of this paper is presenting a flexible solution to the box-constrained least squares problems. This solution is applicable to many existing problems, such as nonnegative matrix factorization, support vector machine, signal deconvolution, and computed tomography reconstruction. The key concept of the proposed algorithm is to replace the minimization of the cost function at each iteration by the minimization of a surrogate, leading to a guaranteed decrease in the cost function. In addition to the monotonicity, the proposed algorithm also owns a few good features including the self-constraint in the feasible region and the absence of a predetermined step size. This paper theoretically proves the global convergence for a special case of below-bounded constraints. Using the proposed mechanism, some valuable algorithms can be derived. The simulation results demonstrate that the proposed algorithm provides performance that is comparable to that of other commonly used methods in numerical experiment and computed tomography reconstruction.


---

## Body

## 1. Introduction

Solving the linear systemY
=
A
X is a classic inverse problem, where X
∈
R
N and Y
∈
R
M are vectors and A
∈
R
M
×
N is a matrix. This problem is applicable in many fields, including nonnegative matrix factorization (NMF), support vector machine (SVM), signal deconvolution, and medical image reconstruction (e.g., computed tomography (CT)) [1–5]. In many cases, we need to impose a box constraint on the problem. A classical approach to choosing X is to minimize the least squares (LS) error between Y and A
X: (1)
min
C
1
≤
X
≤
C
2
⁡
F
X
=
A
X
-
Y
2
,where C
1 and C
2 are given constant vectors and · denotes the Euclidean distance.In fact, (1) can be easily converted into an equivalent form by a linear transformation ξ
=
X
-
C
1, which is convenient to study. In the following, we only focus on such LS problems and still denote variables with X for simple notation; subsequently, the readers may generalize all the results of this paper to the generic form by themselves:(2)
min
0
≤
X
≤
C
⁡
F
X
=
A
X
-
Y
2
,where C is a constant vector.If we denoteC
=
+
∞, then it is the widely used case, that is, the nonnegativity constrained LS problem: (3)
min
X
≥
0
⁡
F
X
=
A
X
-
Y
2
.There are many methods for solving nonnegativity constrained problems. A simple method for solving such problems is by using(4)
X
=
A
T
A
-
1
A
T
Y
,
X
=
max
⁡
X
,
0
.Of course, we must truncate the negative pixel values to satisfy the nonnegativity constraint. However, we will encounter two difficulties: perhaps A
T
A is not invertible or A
T
A is too large. A viable example for the minimization of (3) is the gradient-based method. The steepest descent method is perhaps the simplest technique to implement, which takes the negative gradient as the descent direction: (5)
X
t
+
1
=
X
t
-
α
t
∇
F
X
t
,
X
t
+
1
=
max
⁡
X
t
+
1
,
0
,where the superscript (t) denotes the tth iteration and α
(
t
) is the step size.Letα
(
t
)
=
α be a constant, where 0
<
α
<
2
/
∇
2
F
X
t
E, and ·
E denotes the maximum eigenvalue for a matrix; then, (5) becomes the projected Landweber (PL) method [6]: (6)
X
t
+
1
=
X
t
-
α
∇
F
X
t
,
X
t
+
1
=
max
⁡
X
t
+
1
,
0
.In addition, Lantéri et al. [7, 8] provided a general multiplicative algorithm, which was also a gradient-based method. Let a
j
(
X
t
) be a positive function with positive values for any X
j
t
>
0; then, -
diag
⁡
[
X
j
t
a
j
(
X
t
)
]
∇
F
(
X
t
) is also a descent direction. Therefore, the algorithm in a modified gradient form can be written as (7)
X
j
t
+
1
=
X
j
t
-
X
j
t
a
j
X
t
∇
j
F
X
t
.Again, let U
j
(
X
t
) and V
j
(
X
t
) be two positive functions for any X
j
t
>
0, which satisfy -
∇
j
F
(
X
t
)
=
U
j
(
X
t
)
-
V
j
(
X
t
). Taking a
j
(
X
t
)
=
1
/
V
j
(
X
t
), Lantéri et al. obtained a multiplicative update (8)
X
j
t
+
1
=
X
j
t
U
j
X
t
V
j
X
t
.The conjugate gradient (CG) method is also a popular method, and it is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation.The above approaches are not effective for processing the upper-bounded constraints, and, generally, a truncation is imposed to the update rule, which may lead to a divergent modification.In this paper, we consider a specific application of the surrogate-based methods [9] to a specific type of LS problem with box constraints. We derive a multiplicative update rule to iteratively and monotonically (in the sense of decreasing the cost function) solve the problem, similar to the EM algorithm [10]. The nonnegativity constraints will be satisfied automatically, and the upper-bounded constraints are performed by an upper truncation. Meanwhile, the algorithm still has the monotonicity and does not require an adjustable step size. We provide a rigorous global convergence proof for the case of only nonnegativity constraints, which can be easily generalized to the below-constrained case. We demonstrate the power of this mechanism by deriving many existing and new algorithms. We also present some computer simulation results to show the desirable behavior of the proposed algorithm with respect to the convergence rate and the stability compared with existing methods.


## 2. Methodology

A note about our notation: all vectors will be column vectors unless transposed to a row vector by a prime superscriptT. For a matrix or vector X, X
≥
0 means that any component of X is equal to or greater than 0. For a matrix A, A
i
∗ and A
∗
i represent the ith row and ith column of A, respectively.As mentioned in many articles [1, 2, 9, 11], a surrogate as defined below is useful in algorithm derivations and convergence proofs.Definition 1 (surrogate).
Denoteψ
(
X
∣
X
t
) as a surrogate for Ψ
(
X
) at X
t (fixed) if ψ
(
X
t
∣
X
t
)
=
Ψ
(
X
t
) and ψ
(
X
∣
X
t
)
≥
Ψ
(
X
).Clearly,Ψ
(
X
) is decreasing under the update X
t
+
1
=
m
i
n
X
ψ
(
X
∣
X
t
) because (9)
Ψ
X
t
+
1
≤
ψ
X
t
+
1
∣
X
t
≤
ψ
X
t
∣
X
t
=
Ψ
X
t
.There are two obvious and important properties for a surrogate: additivity and transitivity. For the former, the sum of two surrogates is a surrogate of the sum of the two original functions. For the latter, the surrogate of a surrogate of a function is a surrogate of this function.In the following, we proceed step by step: firstly, assume thatA
≥
0 and C
=
+
∞; secondly, remove the restriction of A
≥
0; and thirdly, remove the restriction of C
=
+
∞.


### 2.1.A
≥
0 and C
=
+
∞

Let-
Y
=
Y
+
-
Y
-, where Y
+ and Y
- are nonnegative vectors; then, F
(
X
)
=
A
X
+
Y
+
-
Y
-
2. Note that we decompose -
Y instead of Y for convenient derivation, which can be seen below. We construct a surrogate f
(
X
∣
X
t
) by the convexity of F
(
X
). Denote (10)
λ
i
∗
=
Y
i
+
A
X
t
+
Y
+
i
,
λ
i
j
=
A
i
j
X
j
t
A
X
t
+
Y
+
ithat satisfy λ
i
∗, λ
i
j
≥
0 and λ
i
∗
+
∑
j
=
1
N
‍
λ
i
j
=
1. They can be the convex combination coefficients such that (11)
f
X
∣
X
t
=
∑
i
=
1
M
λ
i
∗
Y
i
+
λ
i
∗
-
Y
i
-
2
+
∑
j
=
1
N
λ
i
j
A
i
j
X
j
λ
i
j
-
Y
i
-
2
.It is easy to verify thatf
(
X
t
∣
X
t
)
=
F
(
X
t
). If considering Jensen’s inequality and the convex combination coefficients λ
i
j, then f
(
X
∣
X
t
)
≥
F
(
X
) is proven by the following inequality: (12)
λ
i
∗
Y
i
+
λ
i
∗
-
Y
i
-
2
+
∑
j
=
1
N
λ
i
j
A
i
j
X
j
λ
i
j
-
Y
i
-
2
≥
A
X
i
+
Y
i
+
-
Y
i
-
2
.Take the partial derivatives off
(
X
∣
X
t
); then, we can solve the one-dimensional equations ∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 to obtain a multiplicative update rule: (13)
X
j
t
+
1
=
X
j
t
A
T
Y
-
j
A
T
A
X
t
+
Y
+
j
.


### 2.2. AnyA but C
=
+
∞

LetA
=
A
+
-
A
-, w
1
+
w
2
=
1, and w
1
,
w
2
>
0, where A
+ and A
- are two nonnegative matrices; then, we can construct the surrogate f
m
i
d
(
X
∣
X
t
) for F
(
X
) as follows: (14)
f
m
i
d
X
∣
X
t
=
f
-
X
∣
X
t
+
f
^
X
∣
X
t
,where(15)
f
-
X
∣
X
t
=
w
1
1
w
1
A
+
X
+
Y
+
-
1
w
1
-
1
A
+
X
t
+
Y
+
-
A
-
X
t
+
Y
-
2
,
f
^
X
∣
X
t
=
w
2
1
w
2
A
-
X
+
Y
-
-
A
+
X
t
+
Y
+
-
1
w
2
-
1
A
-
X
t
+
Y
-
2
.It is easy to verify thatf
m
i
d
(
X
t
∣
X
t
)
=
F
(
X
t
). By the convexity of F
(
X
), it is clear that(16)
f
m
i
d
X
∣
X
t
=
w
1
1
w
1
A
+
X
+
Y
+
-
1
w
1
-
1
A
+
X
t
+
Y
+
-
A
-
X
t
+
Y
-
2
+
w
2
-
1
w
2
A
-
X
+
Y
-
+
A
+
X
t
+
Y
+
+
1
w
2
-
1
A
-
X
t
+
Y
-
2
≥
F
X
n
o
t
e
t
h
a
t
w
1
+
w
2
=
1
.Following the same process as in Section2.1, we can construct surrogates for f
-
(
X
∣
X
t
) and f
^
(
X
∣
X
t
). Let (17)
λ
i
j
+
=
A
i
j
+
X
j
t
A
+
X
t
+
Y
+
i
,
λ
i
∗
+
=
Y
i
+
A
+
X
t
+
Y
+
i
,
λ
i
∗
-
=
Y
i
-
A
-
X
t
+
Y
-
i
,
λ
i
j
-
=
A
i
j
-
X
j
t
A
-
X
t
+
Y
-
i
,
B
-
=
1
w
1
-
1
A
+
X
t
+
Y
+
+
A
-
X
t
+
Y
-
,
B
^
=
A
+
X
t
+
Y
+
+
1
w
2
-
1
A
-
X
t
+
Y
-
;then, (18)
f
̿
X
∣
X
t
=
w
1
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
+
1
w
1
A
i
j
+
X
j
λ
i
j
+
-
B
-
i
2
+
c
o
n
s
t
,
f
^
^
X
∣
X
t
=
w
2
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
-
1
w
2
A
i
j
-
X
j
λ
i
j
-
-
B
^
i
2
+
c
o
n
s
t
.Note that λ
i
∗
+ and λ
i
∗
- are relative to the constant terms, which do not work for the optimization problems; thus, we ignore them. We can now obtain a surrogate for F
(
X
) at X
t: (19)
f
X
∣
X
t
=
f
̿
X
∣
X
t
+
f
^
^
X
∣
X
t
.Take the partial derivatives off
̿
(
X
∣
X
t
) and f
^
^
(
X
∣
X
t
): (20)
∂
f
̿
X
∣
X
t
∂
X
j
=
2
X
j
w
1
X
j
t
A
+
T
A
+
X
t
+
Y
+
j
-
2
A
+
T
B
-
j
,
∂
f
^
^
X
∣
X
t
∂
X
j
=
2
X
j
w
2
X
j
t
A
-
T
A
-
X
t
+
Y
-
j
-
2
A
-
T
B
^
j
.Then, we solve the one-dimensional equations∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 to obtain a multiplicative update rule: (21)
X
j
t
+
1
=
X
j
t
w
1
w
2
A
+
T
B
-
+
A
-
T
B
^
j
f
f
w
2
A
+
T
A
+
X
t
+
Y
+
j
+
w
1
A
-
T
A
-
X
t
+
Y
-
j
=
X
j
t
w
2
A
+
+
w
1
A
-
T
w
2
A
+
X
t
+
Y
+
+
w
1
A
-
X
t
+
Y
-
j
f
f
w
2
A
+
T
A
+
X
t
+
Y
+
j
+
w
1
A
-
T
A
-
X
t
+
Y
-
j
.


### 2.3. AnyA and C

As shown, (21) will ensure the monotonic decrease of the cost function and satisfy the nonnegativity constraints. In the following, it will be modified by a truncation at the end of each iteration to satisfy the box constraints: (22)
X
=
min
⁡
X
,
C
.We will show that the truncation still ensures the monotonic decrease of the cost function.It is easy to see the separability of the variables of (19); that is, f
(
X
,
X
t
)
=
∑
j
=
1
N
‍
f
j
(
X
j
,
X
t
), and every separated function f
j
(
X
j
,
X
t
) has a quadratic form: (23)
f
j
X
j
,
X
t
=
w
1
∑
i
=
1
M
λ
i
j
+
1
w
1
A
i
j
+
X
j
λ
i
j
+
-
B
-
i
2
+
w
2
∑
i
=
1
M
λ
i
j
-
1
w
2
A
i
j
-
X
j
λ
i
j
-
-
B
^
i
2
+
c
o
n
s
t
.We further simplify it into a general form:(24)
g
X
j
=
a
X
j
-
b
2
+
c
o
n
s
t
,where a
>
0, b
>
0, and X
j
t
∈
0
C. It is clear that X
j
=
b minimizes the quadratic function. If b
≤
C, we change nothing. If b
>
C, then it is easy to obtain that g
(
X
j
t
)
≥
g
(
C
) because g
′
(
X
j
)
=
2
a
(
X
j
-
b
)
<
0 on the open interval (X
j
t
b
) such that g
(
X
j
) strictly monotonously decreases on that open interval. Thus, the proof is established.In fact, we know thatX
=
m
i
n
{
X
,
+
∞
}; therefore, we may provide a generic solution to (2), which is the main result of this paper.Algorithm 2.
Start from an initial pointX
0
>
0; then, (1)
updateX by (21);
(2)
truncateX by (22).It is easy to generalize the constraints from0
≤
X
≤
C to C
1
≤
X
≤
C
2 by the transformation ξ
=
X
-
C
1, which is left for readers to complete by themselves.


## 2.1.A
≥
0 and C
=
+
∞

Let-
Y
=
Y
+
-
Y
-, where Y
+ and Y
- are nonnegative vectors; then, F
(
X
)
=
A
X
+
Y
+
-
Y
-
2. Note that we decompose -
Y instead of Y for convenient derivation, which can be seen below. We construct a surrogate f
(
X
∣
X
t
) by the convexity of F
(
X
). Denote (10)
λ
i
∗
=
Y
i
+
A
X
t
+
Y
+
i
,
λ
i
j
=
A
i
j
X
j
t
A
X
t
+
Y
+
ithat satisfy λ
i
∗, λ
i
j
≥
0 and λ
i
∗
+
∑
j
=
1
N
‍
λ
i
j
=
1. They can be the convex combination coefficients such that (11)
f
X
∣
X
t
=
∑
i
=
1
M
λ
i
∗
Y
i
+
λ
i
∗
-
Y
i
-
2
+
∑
j
=
1
N
λ
i
j
A
i
j
X
j
λ
i
j
-
Y
i
-
2
.It is easy to verify thatf
(
X
t
∣
X
t
)
=
F
(
X
t
). If considering Jensen’s inequality and the convex combination coefficients λ
i
j, then f
(
X
∣
X
t
)
≥
F
(
X
) is proven by the following inequality: (12)
λ
i
∗
Y
i
+
λ
i
∗
-
Y
i
-
2
+
∑
j
=
1
N
λ
i
j
A
i
j
X
j
λ
i
j
-
Y
i
-
2
≥
A
X
i
+
Y
i
+
-
Y
i
-
2
.Take the partial derivatives off
(
X
∣
X
t
); then, we can solve the one-dimensional equations ∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 to obtain a multiplicative update rule: (13)
X
j
t
+
1
=
X
j
t
A
T
Y
-
j
A
T
A
X
t
+
Y
+
j
.


## 2.2. AnyA but C
=
+
∞

LetA
=
A
+
-
A
-, w
1
+
w
2
=
1, and w
1
,
w
2
>
0, where A
+ and A
- are two nonnegative matrices; then, we can construct the surrogate f
m
i
d
(
X
∣
X
t
) for F
(
X
) as follows: (14)
f
m
i
d
X
∣
X
t
=
f
-
X
∣
X
t
+
f
^
X
∣
X
t
,where(15)
f
-
X
∣
X
t
=
w
1
1
w
1
A
+
X
+
Y
+
-
1
w
1
-
1
A
+
X
t
+
Y
+
-
A
-
X
t
+
Y
-
2
,
f
^
X
∣
X
t
=
w
2
1
w
2
A
-
X
+
Y
-
-
A
+
X
t
+
Y
+
-
1
w
2
-
1
A
-
X
t
+
Y
-
2
.It is easy to verify thatf
m
i
d
(
X
t
∣
X
t
)
=
F
(
X
t
). By the convexity of F
(
X
), it is clear that(16)
f
m
i
d
X
∣
X
t
=
w
1
1
w
1
A
+
X
+
Y
+
-
1
w
1
-
1
A
+
X
t
+
Y
+
-
A
-
X
t
+
Y
-
2
+
w
2
-
1
w
2
A
-
X
+
Y
-
+
A
+
X
t
+
Y
+
+
1
w
2
-
1
A
-
X
t
+
Y
-
2
≥
F
X
n
o
t
e
t
h
a
t
w
1
+
w
2
=
1
.Following the same process as in Section2.1, we can construct surrogates for f
-
(
X
∣
X
t
) and f
^
(
X
∣
X
t
). Let (17)
λ
i
j
+
=
A
i
j
+
X
j
t
A
+
X
t
+
Y
+
i
,
λ
i
∗
+
=
Y
i
+
A
+
X
t
+
Y
+
i
,
λ
i
∗
-
=
Y
i
-
A
-
X
t
+
Y
-
i
,
λ
i
j
-
=
A
i
j
-
X
j
t
A
-
X
t
+
Y
-
i
,
B
-
=
1
w
1
-
1
A
+
X
t
+
Y
+
+
A
-
X
t
+
Y
-
,
B
^
=
A
+
X
t
+
Y
+
+
1
w
2
-
1
A
-
X
t
+
Y
-
;then, (18)
f
̿
X
∣
X
t
=
w
1
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
+
1
w
1
A
i
j
+
X
j
λ
i
j
+
-
B
-
i
2
+
c
o
n
s
t
,
f
^
^
X
∣
X
t
=
w
2
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
-
1
w
2
A
i
j
-
X
j
λ
i
j
-
-
B
^
i
2
+
c
o
n
s
t
.Note that λ
i
∗
+ and λ
i
∗
- are relative to the constant terms, which do not work for the optimization problems; thus, we ignore them. We can now obtain a surrogate for F
(
X
) at X
t: (19)
f
X
∣
X
t
=
f
̿
X
∣
X
t
+
f
^
^
X
∣
X
t
.Take the partial derivatives off
̿
(
X
∣
X
t
) and f
^
^
(
X
∣
X
t
): (20)
∂
f
̿
X
∣
X
t
∂
X
j
=
2
X
j
w
1
X
j
t
A
+
T
A
+
X
t
+
Y
+
j
-
2
A
+
T
B
-
j
,
∂
f
^
^
X
∣
X
t
∂
X
j
=
2
X
j
w
2
X
j
t
A
-
T
A
-
X
t
+
Y
-
j
-
2
A
-
T
B
^
j
.Then, we solve the one-dimensional equations∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 to obtain a multiplicative update rule: (21)
X
j
t
+
1
=
X
j
t
w
1
w
2
A
+
T
B
-
+
A
-
T
B
^
j
f
f
w
2
A
+
T
A
+
X
t
+
Y
+
j
+
w
1
A
-
T
A
-
X
t
+
Y
-
j
=
X
j
t
w
2
A
+
+
w
1
A
-
T
w
2
A
+
X
t
+
Y
+
+
w
1
A
-
X
t
+
Y
-
j
f
f
w
2
A
+
T
A
+
X
t
+
Y
+
j
+
w
1
A
-
T
A
-
X
t
+
Y
-
j
.


## 2.3. AnyA and C

As shown, (21) will ensure the monotonic decrease of the cost function and satisfy the nonnegativity constraints. In the following, it will be modified by a truncation at the end of each iteration to satisfy the box constraints: (22)
X
=
min
⁡
X
,
C
.We will show that the truncation still ensures the monotonic decrease of the cost function.It is easy to see the separability of the variables of (19); that is, f
(
X
,
X
t
)
=
∑
j
=
1
N
‍
f
j
(
X
j
,
X
t
), and every separated function f
j
(
X
j
,
X
t
) has a quadratic form: (23)
f
j
X
j
,
X
t
=
w
1
∑
i
=
1
M
λ
i
j
+
1
w
1
A
i
j
+
X
j
λ
i
j
+
-
B
-
i
2
+
w
2
∑
i
=
1
M
λ
i
j
-
1
w
2
A
i
j
-
X
j
λ
i
j
-
-
B
^
i
2
+
c
o
n
s
t
.We further simplify it into a general form:(24)
g
X
j
=
a
X
j
-
b
2
+
c
o
n
s
t
,where a
>
0, b
>
0, and X
j
t
∈
0
C. It is clear that X
j
=
b minimizes the quadratic function. If b
≤
C, we change nothing. If b
>
C, then it is easy to obtain that g
(
X
j
t
)
≥
g
(
C
) because g
′
(
X
j
)
=
2
a
(
X
j
-
b
)
<
0 on the open interval (X
j
t
b
) such that g
(
X
j
) strictly monotonously decreases on that open interval. Thus, the proof is established.In fact, we know thatX
=
m
i
n
{
X
,
+
∞
}; therefore, we may provide a generic solution to (2), which is the main result of this paper.Algorithm 2.
Start from an initial pointX
0
>
0; then, (1)
updateX by (21);
(2)
truncateX by (22).It is easy to generalize the constraints from0
≤
X
≤
C to C
1
≤
X
≤
C
2 by the transformation ξ
=
X
-
C
1, which is left for readers to complete by themselves.


## 3. Convergence

The convergence proof with box constraints is very difficult; thus, we only focus on the below-bounded constraints, which can be viewed as a special case of box constraints withC
=
+
∞. Additionally, we consider the equivalence between the general below-constrained form and the nonnegative constrained one by a linear transformation ξ
=
X
-
C
1; then, we will theoretically prove the global convergence of the latter case for simplification, which is easy to extend to that of the former. In theory, the Kuhn-Tucker (KT) point will be a global solution if the cost function is convex. It is easy to see the convexity of F
(
X
). By Theorem 2.19 in [12], the KT conditions of (3) (X
≥
0) are as follows: (25)
∂
F
X
∂
X
j
=
0
i
f
X
j
>
0
,
(26)
∂
F
X
∂
X
j
≥
0
i
f
X
j
=
0
.We entertain several important and reasonable assumptions.Assumption 3.
For the iteration sequence{
X
t
}, we assume that (1)
the algorithm starts from a positive image;
(2)
(
A
T
A
)
j
j
>
0 for all j;
(3)
F is a strictly convex function.The first assumption forces the iterations to be positive, but the limit may be zero. The second condition is reasonable because(
A
T
A
)
j
j
=
0 suggests A
i
j
=
0 for any i. Thus, the equation Y
=
A
X is irrelevant to X
j, and, then, X
j in X is removable. The third one is indeed restrictive; however, it is important for the following derivation.We will prove the global convergence along the lines of [13–17]. First, we provide several useful lemmas.Lemma 4.
The set of accumulation points of a bounded sequence{
Z
t
} with {
Z
t
+
1
-
Z
t
}
→
0 is connected and compact.Proof.
This is Theorem 28.1 from Ostrowski [18]. The reader is kindly referred to this paper for the proof.Lemma 5.
The iteration sequence{
X
t
} is bounded.Proof.
BecauseF is a strictly convex function, A
T
A is a positive definite symmetric matrix, and, thus, it can be factorized into A
T
A
=
U
T
U by Cholesky’s method, where U is an invertible matrix. We assume that X
∗ is the sole minimum, and, then, we expand F on X
∗ using a Taylor series: (27)
F
X
-
F
X
∗
=
X
-
X
∗
T
A
T
A
X
-
X
∗
=
X
-
X
∗
T
U
T
U
X
-
X
∗
≥
0
.Let ξ
=
U
(
X
-
X
∗
); then, for any constant L, it is easy to know the boundedness of S
=
{
ξ
:
ξ
T
ξ
<
L
} such that {
X
:
X
=
U
-
1
ξ
+
X
∗
,
ξ
∈
S
} is bounded, which is equivalent to {
X
:
F
(
X
)
-
F
(
X
∗
)
<
L
} being bounded. Let L
0
=
F
(
X
0
)
-
F
(
X
∗
); then, {
X
t
}
⊂
{
X
:
F
(
X
)
-
F
(
X
∗
)
<
L
0
} because {
F
(
X
t
)
} monotonically decreases. Therefore, {
X
t
} is bounded.Lemma 6.
The sequence{
X
t
-
X
t
+
1
}
→
0.Proof.
Take into account that(28)
∇
j
j
2
f
X
∣
X
t
=
2
w
1
X
j
t
A
+
T
A
+
X
t
+
Y
+
j
+
2
w
2
X
j
t
A
-
T
A
-
X
t
+
Y
-
j
≥
2
w
1
X
j
t
A
+
T
A
+
X
t
j
+
2
w
2
X
j
t
A
-
T
A
-
X
t
j
≥
2
w
1
∑
i
=
1
M
A
i
j
+
2
+
2
w
2
∑
i
=
1
M
A
i
j
-
2
=
γ
j
>
0
s
e
e
t
h
e
s
e
c
o
n
d
o
n
e
o
f
A
s
s
u
m
p
t
i
o
n
3
.Let γ
=
m
i
n
j
{
γ
j
}. Because ∇
f
(
X
t
+
1
∣
X
t
)
=
0, (29)
F
X
t
-
F
X
t
+
1
≥
f
X
t
∣
X
t
-
f
X
t
+
1
∣
X
t
=
1
2
X
t
-
X
t
+
1
T
∇
2
f
X
t
+
1
∣
X
t
X
t
-
X
t
+
1
≥
γ
2
X
t
-
X
t
+
1
2
.Because {
F
(
X
t
)
} monotonically decreases and is bounded from below, {
F
(
X
t
)
-
F
(
X
t
+
1
)
}
→
0; therefore, {
X
t
-
X
t
+
1
}
→
0.Lemma 7.
If a subsequence{
X
t
s
}
→
X
∗, then {
X
t
s
+
1
}
→
X
∗ as well.Proof.
By contradiction, if{
X
t
s
+
1
} diverges, then it must have a convergent subsequence {
X
t
s
s
+
1
}
→
X
∗
∗
≠
X
∗ because of the boundness by Lemma 5. Let ϵ
0
=
X
∗
-
X
∗
∗
>
0, and if we consider the two convergent subsequences {
X
t
s
s
} and {
X
t
s
s
+
1
}, then there must be a positive integer T to make X
t
s
s
-
X
∗
<
ϵ
0
/
4 and X
t
s
s
+
1
-
X
∗
∗
<
ϵ
0
/
4 when t
s
s
>
T. By the triangle inequality, we can obtain the contradictive result to {
X
t
-
X
t
+
1
}
→
0 as (30)
X
t
s
s
-
X
t
s
s
+
1
+
X
t
s
s
-
X
∗
+
X
t
s
s
+
1
-
X
∗
∗
≥
X
∗
-
X
∗
∗
⟹
X
t
s
s
-
X
t
s
s
+
1
>
ϵ
0
2
.Lemma 8.
At each iteration, one knows that(31)
X
j
t
+
1
=
X
j
t
-
α
j
∇
F
X
t
j
,
α
j
=
w
1
w
2
X
j
t
2
f
f
w
2
A
+
T
A
+
X
t
+
Y
+
j
+
w
1
A
-
T
A
-
X
t
+
Y
-
j
.Proof.
We consider that(32)
∇
F
X
t
=
2
A
+
-
A
-
T
A
+
X
t
+
Y
+
-
A
-
X
t
+
Y
-
=
2
A
+
T
A
+
X
t
+
Y
+
+
2
A
-
T
A
-
X
t
+
Y
-
-
2
A
+
T
A
-
X
t
+
Y
-
-
2
A
-
T
A
+
X
t
+
Y
+
=
2
w
1
A
+
T
A
+
X
t
+
Y
+
+
2
w
2
A
-
T
A
-
X
t
+
Y
-
-
2
A
+
T
1
w
1
-
1
A
+
X
t
+
Y
+
+
A
-
X
t
+
Y
-
-
2
A
-
T
A
+
X
t
+
Y
+
+
1
w
2
-
1
A
-
X
t
+
Y
-
=
2
w
1
w
2
w
2
A
+
T
A
+
X
t
+
Y
+
+
w
1
A
-
T
A
-
X
t
+
Y
-
-
2
w
1
w
2
w
2
A
+
+
w
1
A
-
T
w
2
A
+
X
t
+
Y
+
+
w
1
A
-
X
t
+
Y
-
(33)
=
∇
f
X
t
∣
X
t
c
o
m
p
a
r
e
32
w
i
t
h
20
.It is easy to verify the correctness of (31) if replacing ∇
F
(
X
t
) by (32).The global convergence will be proven from the three theorems below.Theorem 9.
Let{
X
t
s
}
→
X
∗ be any convergent subsequence; then, X
∗ meets the first KT condition (25).Proof.
WhenX
j
∗
>
0, it is easy to obtain that [
∇
F
(
X
∗
)
]
j
=
[
∇
f
(
X
∗
∣
X
∗
)
]
j by (33). In addition, ∂
F
(
X
∣
X
t
s
)
/
∂
X
j
|
X
t
s
+
1
=
0 and {
X
t
s
}
→
X
∗ such that {
X
t
s
+
1
}
→
X
∗ (Lemma 7); thus, (34)
∂
F
X
∣
X
∗
∂
X
j
X
∗
=
lim
t
s
→
+
∞
⁡
∂
F
X
∣
X
s
t
∂
X
j
X
t
s
+
1
=
0
.Theorem 10.
The entire sequence{
X
t
} converges.Proof.
According to Lemmas4, 5, and 6, the set of accumulation points of {
X
t
} is connected and compact. If we can prove that the number of accumulation points is finite, then the desired result follows because a finite set can be connected only if it consists of a single point [16].
To prove the existence of a finite number of accumulation points, we consider any accumulation pointX
∗. Given an integer set Ω
=
{
1,2
,
…
,
T
}, where T is the total number of components of X, Ω
∗
=
{
j
:
X
j
∗
=
0
} is a subset of Ω. Let F
Ω
∗ be the restrictions of F to the set {
X
:
X
j
=
0
,
j
∈
Ω
∗
}, which is a strictly convex function of the reduced variables. It follows that F
Ω
∗ has a unique minimum (Theorem 9: ∂
F
(
X
∗
)
/
∂
X
j
=
0
i
f
X
j
∗
>
0). This means that an accumulation point must correspond to a subset of Ω. The number of subsets of Ω is finite, and, thus, the number of accumulation points is also finite.In Theorem9, we prove that every accumulation point meets the first KT condition, by which the full sequence convergence is provided in Theorem 10. Naturally, the limit of {
X
t
} satisfies the first KT condition. In the following, we will show that the second KT condition is satisfied.Theorem 11.
The limitX
∗ of {
X
t
} satisfies the second KT condition (26).Proof.
WhenX
j
∗
=
0, by contradiction, we assume that there is an X
j
∗
=
0 that satisfies [
∇
F
(
X
∗
)
]
j
<
0. Because {
X
t
}
→
X
∗, there exists an ϵ
<
0 and a positive integer T such that [
∇
F
(
X
t
)
]
j
<
ϵ for t
>
T; then, (35)
X
j
t
+
1
-
X
j
t
=
-
α
j
∇
F
X
t
j
>
-
α
j
ϵ
>
0
b
y
L
e
m
m
a
8
.Thus, we can obtain that X
j
t
+
1
>
X
j
t, which is a contradiction to {
X
j
t
}
→
0.


## 4. Solutions to Some Existing and New Examples

### 4.1. NMF

NMF has been widely used in pattern recognition, machine learning, and data mining. It decomposes the nonnegative matrixY by the product of two other nonnegative matrices A and X. Lee and Seung [1] have proposed using the square of the Euclidean distance to measure the similarity between Y and A
X: (36)
min
A
≥
0
,
X
≥
0
⁡
F
A
,
X
=
Y
-
A
X
2
.Minimizing them under the constraintsA
≥
0 and X
≥
0, Lee and Seung [1] alternated between solving an optimization problem in the variables A and then solving another one in the variables X. As can be seen, A and X have the same roles.To takeX, for instance, with X
t
>
0 and A
>
0 given, denote λ
i
j
k
=
A
i
k
t
X
k
j
t
/
(
A
X
t
)
i
j such that the surrogate function can be constructed: (37)
f
X
∣
X
t
=
∑
i
j
k
λ
i
j
k
A
i
k
X
k
j
λ
i
j
k
-
Y
i
j
2
.The update rule can be obtained by solving ∂
f
(
X
,
X
t
)
/
∂
X
k
j
=
0. Reversing the roles of A and X, we can similarly construct the surrogate for A and acquire the update rule:(38)
X
k
j
t
+
1
=
X
k
j
t
A
T
Y
k
j
A
T
A
X
t
k
j
,
A
i
k
t
+
1
=
A
i
k
t
Y
X
T
i
k
A
t
X
X
T
i
k
.


### 4.2. Variation of Linearized SVM

SVM attempts to separate points belonging to two given sets in real Euclidean spaceR
n by a surface. Several practical applications of SVMs use nonlinear kernels, such as the polynomial and radial basis function kernels. However, in applications such as text classification, linear SVMs are still used because it has been observed that many text classification problems are linearly separable [19]. Most literature on large-scale SVM training have targeted the linear SVM problem, citing this fact. A dual form of linear SVM is usually used because of the simple structure, which can be formulated as [20](39)
min
⁡
⁡
1
2
X
T
D
A
T
A
D
X
-
1
T
X
s
.
t
.
0
≤
X
≤
C
,
D
T
X
=
0
,where A is the input data, D is the class label, C is the penalty parameter, and X represents the Lagrangian multiplier that needs to be optimized. Hsieh et al. [21] handled the constraint D
T
X
=
0 by removing it, which corresponded to removing the “intercept” from the classifier. For a simplified expression, we denote B
=
A
D, and, then, we obtain the following simple formula:(40)
min
⁡
1
2
X
T
B
T
B
X
-
1
T
X
s.t.
0
≤
X
≤
C
.We can solve this optimization problem using the proposed method. Letw
1
=
w
2
=
1
/
2 and B
=
B
+
-
B
-, where B
+ and B
- are two nonnegative matrices; then, (41)
f
m
i
d
X
∣
X
t
=
f
-
X
∣
X
t
+
f
^
X
∣
X
t
-
1
T
X
,where(42)
f
-
X
∣
X
t
=
1
4
2
B
+
X
-
B
+
X
t
-
B
-
X
t
2
,
f
^
X
∣
X
t
=
1
4
2
B
-
X
-
B
+
X
t
-
B
-
X
t
2
.We, respectively, construct surrogates forf
-
(
X
∣
X
t
) and f
^
(
X
∣
X
t
). Let (43)
λ
i
j
+
=
B
i
j
+
X
j
t
B
+
X
t
i
,
λ
i
j
-
=
B
i
j
-
X
j
t
B
-
X
t
i
;then, (44)
f
̿
X
∣
X
t
=
1
4
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
+
2
B
i
j
+
X
j
λ
i
j
+
-
B
+
X
t
+
B
-
X
t
i
2
,
f
^
^
X
∣
X
t
=
1
4
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
-
2
B
i
j
-
X
j
λ
i
j
-
-
B
+
X
t
+
B
-
X
t
i
2
.We now obtain a surrogate forF
(
X
) at X
t: (45)
f
X
∣
X
t
=
f
̿
X
∣
X
t
+
f
^
^
X
∣
X
t
.We solve the one-dimensional equation∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 and truncate it to obtain an update: (46)
X
j
t
+
1
=
X
j
t
1
+
B
+
+
B
-
T
B
+
+
B
-
X
t
j
2
B
+
T
B
+
+
B
-
T
B
-
X
t
j
,
X
j
t
+
1
=
min
⁡
X
j
t
+
1
,
C
.


### 4.3. Nonnegative Image Deblurring

In many optical devices, the process of image blurring can be considered to be the result of convolution by a point spread function (PSF) [22, 23]. It is assumed that the degraded image Y is in the form Y
=
A
X, where X is the true image and A is the PSF matrix consisting of PSFs at every pixel. As noted in [24, 25], a simple way to approach the deconvolution problem is to find the least squares (LS) estimation between Y and A
X. It is well known that the problem of restoring the original image from the noisy and degraded version is an ill-posed inverse problem: small perturbations in the data may result in an unacceptable result [26]. The Tikhonov regularization method [27, 28] is a popular method that generally leads to a unique solution, which is formulated as (47)
min
X
≥
0
⁡
F
X
=
A
X
-
Y
2
+
β
R
X
2
.Note that we must approve negative values existing in the matrix R.In [17], we assume that A and Y are a nonnegative matrix and vector and that there is no restriction on R; then, we, respectively, construct surrogate functions for A
X
-
Y
2 and R
X
2 using the method proposed in this paper such that we obtain a multiplicative update rule as follows:(48)
X
j
t
+
1
=
X
j
t
A
T
Y
j
+
β
R
+
+
R
-
T
R
+
+
R
-
X
t
j
A
T
A
X
t
j
+
2
β
R
+
T
R
+
+
R
-
T
R
-
X
t
j
.If β
=
0, it is De Pierro’s ISRT (Image Space Reconstruction Technique) [4]: (49)
X
j
t
+
1
=
X
j
t
A
T
Y
j
A
T
A
X
t
j
.


### 4.4. CT Reconstruction withL
1 Regularization

In CT reconstruction, letA be the system probability matrix, Y projections, and R a predetermined matrix. Then, a classical approach is to choose an X (X-ray attenuation image) such that the LS error [5] between Y and A
X is minimized: (50)
min
X
≥
0
⁡
F
X
=
A
X
-
Y
2
+
β
R
X
1
,where β
≥
0 serves as a penalty parameter. Many R matrices have been used, such as the first- or second-order derivative matrix or wavelet basis matrix. Note that both A and Y are a nonnegative matrix and vector; however, we must accept negative values existing in the matrix R.The alternating direction method of multipliers (ADMM) [29, 30] has been developed for solving optimization problems. This method decomposes the original problem into three subproblems and then sequentially solves them at each iteration.Introducing the additional variableV
=
R
X, a fixed penalty parameter ρ, and the Lagrangian multiplier μ, a unified framework can then be given to solve the L
1-norm regularized LS reconstruction problem [29].Algorithm 12 (ADMM general framework).
(1) V
i
t
+
1
=
(
|
(
R
X
t
+
μ
t
)
i
|
-
β
/
ρ
)
+
sgn
⁡
[
(
R
X
t
+
μ
t
)
i
]
,
i
=
1
,
…
,
M.
(2) X
t
+
1
=
a
r
g
m
i
n
X
≥
0
L
(
X
,
V
t
+
1
,
μ
t
)
=
A
X
-
Y
2
+
(
ρ
/
2
)
R
X
-
V
t
+
1
+
μ
t
2
2.
(3) μ
t
+
1
=
μ
t
+
R
X
t
+
1
-
V
t
+
1.TheX-update is the most difficult problem that minimizes L
(
X
,
V
t
+
1
,
μ
t
) with nonnegativity constraints. We can, respectively, construct surrogate functions for A
X
-
Y
2 and R
X
-
V
t
+
1
+
μ
t
2
2; then minimize the sum of the surrogates to update X; however, there is a high overlap ratio with the above, so we leave them to the readers.


## 4.1. NMF

NMF has been widely used in pattern recognition, machine learning, and data mining. It decomposes the nonnegative matrixY by the product of two other nonnegative matrices A and X. Lee and Seung [1] have proposed using the square of the Euclidean distance to measure the similarity between Y and A
X: (36)
min
A
≥
0
,
X
≥
0
⁡
F
A
,
X
=
Y
-
A
X
2
.Minimizing them under the constraintsA
≥
0 and X
≥
0, Lee and Seung [1] alternated between solving an optimization problem in the variables A and then solving another one in the variables X. As can be seen, A and X have the same roles.To takeX, for instance, with X
t
>
0 and A
>
0 given, denote λ
i
j
k
=
A
i
k
t
X
k
j
t
/
(
A
X
t
)
i
j such that the surrogate function can be constructed: (37)
f
X
∣
X
t
=
∑
i
j
k
λ
i
j
k
A
i
k
X
k
j
λ
i
j
k
-
Y
i
j
2
.The update rule can be obtained by solving ∂
f
(
X
,
X
t
)
/
∂
X
k
j
=
0. Reversing the roles of A and X, we can similarly construct the surrogate for A and acquire the update rule:(38)
X
k
j
t
+
1
=
X
k
j
t
A
T
Y
k
j
A
T
A
X
t
k
j
,
A
i
k
t
+
1
=
A
i
k
t
Y
X
T
i
k
A
t
X
X
T
i
k
.


## 4.2. Variation of Linearized SVM

SVM attempts to separate points belonging to two given sets in real Euclidean spaceR
n by a surface. Several practical applications of SVMs use nonlinear kernels, such as the polynomial and radial basis function kernels. However, in applications such as text classification, linear SVMs are still used because it has been observed that many text classification problems are linearly separable [19]. Most literature on large-scale SVM training have targeted the linear SVM problem, citing this fact. A dual form of linear SVM is usually used because of the simple structure, which can be formulated as [20](39)
min
⁡
⁡
1
2
X
T
D
A
T
A
D
X
-
1
T
X
s
.
t
.
0
≤
X
≤
C
,
D
T
X
=
0
,where A is the input data, D is the class label, C is the penalty parameter, and X represents the Lagrangian multiplier that needs to be optimized. Hsieh et al. [21] handled the constraint D
T
X
=
0 by removing it, which corresponded to removing the “intercept” from the classifier. For a simplified expression, we denote B
=
A
D, and, then, we obtain the following simple formula:(40)
min
⁡
1
2
X
T
B
T
B
X
-
1
T
X
s.t.
0
≤
X
≤
C
.We can solve this optimization problem using the proposed method. Letw
1
=
w
2
=
1
/
2 and B
=
B
+
-
B
-, where B
+ and B
- are two nonnegative matrices; then, (41)
f
m
i
d
X
∣
X
t
=
f
-
X
∣
X
t
+
f
^
X
∣
X
t
-
1
T
X
,where(42)
f
-
X
∣
X
t
=
1
4
2
B
+
X
-
B
+
X
t
-
B
-
X
t
2
,
f
^
X
∣
X
t
=
1
4
2
B
-
X
-
B
+
X
t
-
B
-
X
t
2
.We, respectively, construct surrogates forf
-
(
X
∣
X
t
) and f
^
(
X
∣
X
t
). Let (43)
λ
i
j
+
=
B
i
j
+
X
j
t
B
+
X
t
i
,
λ
i
j
-
=
B
i
j
-
X
j
t
B
-
X
t
i
;then, (44)
f
̿
X
∣
X
t
=
1
4
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
+
2
B
i
j
+
X
j
λ
i
j
+
-
B
+
X
t
+
B
-
X
t
i
2
,
f
^
^
X
∣
X
t
=
1
4
∑
i
=
1
M
∑
j
=
1
N
λ
i
j
-
2
B
i
j
-
X
j
λ
i
j
-
-
B
+
X
t
+
B
-
X
t
i
2
.We now obtain a surrogate forF
(
X
) at X
t: (45)
f
X
∣
X
t
=
f
̿
X
∣
X
t
+
f
^
^
X
∣
X
t
.We solve the one-dimensional equation∂
f
(
X
∣
X
t
)
/
∂
X
j
=
0 and truncate it to obtain an update: (46)
X
j
t
+
1
=
X
j
t
1
+
B
+
+
B
-
T
B
+
+
B
-
X
t
j
2
B
+
T
B
+
+
B
-
T
B
-
X
t
j
,
X
j
t
+
1
=
min
⁡
X
j
t
+
1
,
C
.


## 4.3. Nonnegative Image Deblurring

In many optical devices, the process of image blurring can be considered to be the result of convolution by a point spread function (PSF) [22, 23]. It is assumed that the degraded image Y is in the form Y
=
A
X, where X is the true image and A is the PSF matrix consisting of PSFs at every pixel. As noted in [24, 25], a simple way to approach the deconvolution problem is to find the least squares (LS) estimation between Y and A
X. It is well known that the problem of restoring the original image from the noisy and degraded version is an ill-posed inverse problem: small perturbations in the data may result in an unacceptable result [26]. The Tikhonov regularization method [27, 28] is a popular method that generally leads to a unique solution, which is formulated as (47)
min
X
≥
0
⁡
F
X
=
A
X
-
Y
2
+
β
R
X
2
.Note that we must approve negative values existing in the matrix R.In [17], we assume that A and Y are a nonnegative matrix and vector and that there is no restriction on R; then, we, respectively, construct surrogate functions for A
X
-
Y
2 and R
X
2 using the method proposed in this paper such that we obtain a multiplicative update rule as follows:(48)
X
j
t
+
1
=
X
j
t
A
T
Y
j
+
β
R
+
+
R
-
T
R
+
+
R
-
X
t
j
A
T
A
X
t
j
+
2
β
R
+
T
R
+
+
R
-
T
R
-
X
t
j
.If β
=
0, it is De Pierro’s ISRT (Image Space Reconstruction Technique) [4]: (49)
X
j
t
+
1
=
X
j
t
A
T
Y
j
A
T
A
X
t
j
.


## 4.4. CT Reconstruction withL
1 Regularization

In CT reconstruction, letA be the system probability matrix, Y projections, and R a predetermined matrix. Then, a classical approach is to choose an X (X-ray attenuation image) such that the LS error [5] between Y and A
X is minimized: (50)
min
X
≥
0
⁡
F
X
=
A
X
-
Y
2
+
β
R
X
1
,where β
≥
0 serves as a penalty parameter. Many R matrices have been used, such as the first- or second-order derivative matrix or wavelet basis matrix. Note that both A and Y are a nonnegative matrix and vector; however, we must accept negative values existing in the matrix R.The alternating direction method of multipliers (ADMM) [29, 30] has been developed for solving optimization problems. This method decomposes the original problem into three subproblems and then sequentially solves them at each iteration.Introducing the additional variableV
=
R
X, a fixed penalty parameter ρ, and the Lagrangian multiplier μ, a unified framework can then be given to solve the L
1-norm regularized LS reconstruction problem [29].Algorithm 12 (ADMM general framework).
(1) V
i
t
+
1
=
(
|
(
R
X
t
+
μ
t
)
i
|
-
β
/
ρ
)
+
sgn
⁡
[
(
R
X
t
+
μ
t
)
i
]
,
i
=
1
,
…
,
M.
(2) X
t
+
1
=
a
r
g
m
i
n
X
≥
0
L
(
X
,
V
t
+
1
,
μ
t
)
=
A
X
-
Y
2
+
(
ρ
/
2
)
R
X
-
V
t
+
1
+
μ
t
2
2.
(3) μ
t
+
1
=
μ
t
+
R
X
t
+
1
-
V
t
+
1.TheX-update is the most difficult problem that minimizes L
(
X
,
V
t
+
1
,
μ
t
) with nonnegativity constraints. We can, respectively, construct surrogate functions for A
X
-
Y
2 and R
X
-
V
t
+
1
+
μ
t
2
2; then minimize the sum of the surrogates to update X; however, there is a high overlap ratio with the above, so we leave them to the readers.


## 5. Experimental Results

We compare the performance of the proposed surrogate method, with those of the PL method (6) and the CG method on simulated data and CT projection data. For the PL method, we use α
=
1
/
(
A
1
A
∞
) [31], where ·
1 and ·
∞ denote the 1- and ∞-norm of a matrix. The code for the CG method comes from [32], which is slightly modified to meet our criteria.The experiments are performed on a HP Compaq PC with a 3.00 GHz Core i5 CPU and 4 GB memory. The algorithms are implemented in MATLAB 7.0. All of the algorithms are initiated by the same uniform image for a fair comparison.The mean square error (MSE) is used to measure the similarity to the true solution, as given below:(51)
M
S
E
t
=
1
N
X
t
-
X
T
r
u
e
2
.The following criterion is applied to stop the iterative process:(52)
χ
t
=
X
t
+
1
-
X
t
X
t
<
ϵ
,where ϵ is a difference tolerance.


### 5.1. Numerical Simulation

Here, we randomly generateX
∈
R
500 and A
∈
R
100
×
500, where the former are uniformly distributed on 0
1 and the latter are on -
1
1. Then we obtain Y by Y
=
A
X. We will minimize F
(
X
)
=
Y
-
A
X
2 under the box constraints 0
≤
X
≤
1. In the experiment, we generate X
,
Y, and A once, but we randomly generate start points (uniformly distributed on 0
1) 50 times to obtain a reliable averaged result.Figure1 presents a comparison of MSE and χ versus iteration number. As shown, CG rapidly finishes iteration to a static solution, so that χ
(
t
)
=
0 after that, which can not be shown with a log-scaled y-axis. However, such a phenomenon does not mean a good result. Instead, it is the worst result among the three algorithms by observing the MSE curve. We can conclude that CG is not suitable to the box-constrained LS problem. PL and the proposed method always decrease the curves of MSE and χ. However, the proposed method shows a superior performance than the PL method.Figure 1
Simulated data:M
S
E and χ versus iteration number: (a) M
S
E and (b) χ.
(a)
(b)Table1 shows the quantitative comparison of the three algorithms; then we can draw the same conclusion as above. This table further proves that the proposed method is rapid and stable, and CG is not suitable to the box-constrained LS problem.Table 1
Performance comparison of the three algorithms.
Iteration
1
50
150
200
250
300
MSE
Proposed
0.2860
0.2226
0.1961
0.1822
0.1737
0.1682
CG
0.3441
0.2240
0.3183
0.3183
0.3183
0.3183
PL
0.3441
0.2398
0.2136
0.2033
0.1978
0.1943
χ
Proposed
0.0031
0.0016
0.0009
0.0006
0.0004
0.0003
CG
1.0000
0.0001
0
0
0
0
PL
1.0000
0.0117
0.0037
0.0018
0.0010
0.0007Table2 shows the computational time of the algorithms. As can be seen, they have similar computational complexity; however, the proposed algorithm requires only a little more time because of the nonnegative decomposition of Y and A. In fact, from the above theoretical analysis, we can observe that the proposed algorithm is also some gradient-based method; thus, it has a similar computational complexity with the CG and PL methods.Table 2
Total running time of the three algorithms with 300 iterations and 50 random start points.
Algorithm
Proposed
CG
PL
Time (second)
2.22
2.13
1.95


### 5.2. CT Reconstruction

CT reconstruction is a medical imaging technique for creating a meaningful diagnostic image. A computer can take the input from the CT machine, run it through the formula, and return a set of images for a physician to examine. Here, we only consider the simplest mathematical modelm
i
n
X
≥
0
Y
-
A
X
2 to pursue the structural image. A simulated thorax phantom with 128
×
128 grids and 0.5 mm pixels, as shown in Figure 2, is used in the following experiments. There are many advantages to using simulated phantoms, including prior knowledge of the pixel values and the ability to control noise. The total attenuation value is approximately 5
×
1
0
7. For this case, the proposed algorithm is called ISRT, which is widely used in the field of CT reconstruction [4].Figure 2
Shepp-Logan phantom with128
×
128 grids: (a) true phantom and (b) noisy projections.
(a)
(b)The system matrix is obtained using the “angle of view” method [10]. From the system matrix, we forward project the phantom on the sinogram with 128 radial bins (bin size of 0.5 mm) and 128 angular views evenly spaced over π. The noisy projections are obtained by the Poisson-random formula y
i
=
P
o
i
s
s
o
n
(
y
i
∗
), where y
i
∗ is the noise-free projection.Figure3 shows the reconstructions with 50 iterations for every algorithm. As shown, the CG’s image suffers from serious noise artifacts. The PL algorithm provides a reconstructed image that is slightly blurred compared with that of ISRT. The proposed ISRT provides the sharpest edge and clearest interior region.Figure 3
Reconstructed images with 50 iterations by (a) the proposed ISRT, (b) CG, and (c) PL.
(a)
(b)
(c)Figure4 quantitatively compares the histograms of the reconstructed images. The CG’s result still shows a clear turbulence around the true value. PL and ISRT provide two comparable curves, and neither is significantly better than the other.Figure 4
Intensity of the true and reconstructed images: (a) horizontal central line and (b) vertical central line.
(a)
(b)We also present a comparison of the MSE curves andχ with increasing iteration number in Figure 5. As shown, CG exhibits an unsteady reconstruction property in which more iterations may lead to a worse result. Both PL and ISRT overcome this shortcoming so that more iterations always provide considerably better results. In addition, ISRT will provide a more rapid convergence than PL.Figure 5
M
S
E and χ versus iteration number: (a) M
S
E and (b) χ.
(a)
(b)


## 5.1. Numerical Simulation

Here, we randomly generateX
∈
R
500 and A
∈
R
100
×
500, where the former are uniformly distributed on 0
1 and the latter are on -
1
1. Then we obtain Y by Y
=
A
X. We will minimize F
(
X
)
=
Y
-
A
X
2 under the box constraints 0
≤
X
≤
1. In the experiment, we generate X
,
Y, and A once, but we randomly generate start points (uniformly distributed on 0
1) 50 times to obtain a reliable averaged result.Figure1 presents a comparison of MSE and χ versus iteration number. As shown, CG rapidly finishes iteration to a static solution, so that χ
(
t
)
=
0 after that, which can not be shown with a log-scaled y-axis. However, such a phenomenon does not mean a good result. Instead, it is the worst result among the three algorithms by observing the MSE curve. We can conclude that CG is not suitable to the box-constrained LS problem. PL and the proposed method always decrease the curves of MSE and χ. However, the proposed method shows a superior performance than the PL method.Figure 1
Simulated data:M
S
E and χ versus iteration number: (a) M
S
E and (b) χ.
(a)
(b)Table1 shows the quantitative comparison of the three algorithms; then we can draw the same conclusion as above. This table further proves that the proposed method is rapid and stable, and CG is not suitable to the box-constrained LS problem.Table 1
Performance comparison of the three algorithms.
Iteration
1
50
150
200
250
300
MSE
Proposed
0.2860
0.2226
0.1961
0.1822
0.1737
0.1682
CG
0.3441
0.2240
0.3183
0.3183
0.3183
0.3183
PL
0.3441
0.2398
0.2136
0.2033
0.1978
0.1943
χ
Proposed
0.0031
0.0016
0.0009
0.0006
0.0004
0.0003
CG
1.0000
0.0001
0
0
0
0
PL
1.0000
0.0117
0.0037
0.0018
0.0010
0.0007Table2 shows the computational time of the algorithms. As can be seen, they have similar computational complexity; however, the proposed algorithm requires only a little more time because of the nonnegative decomposition of Y and A. In fact, from the above theoretical analysis, we can observe that the proposed algorithm is also some gradient-based method; thus, it has a similar computational complexity with the CG and PL methods.Table 2
Total running time of the three algorithms with 300 iterations and 50 random start points.
Algorithm
Proposed
CG
PL
Time (second)
2.22
2.13
1.95


## 5.2. CT Reconstruction

CT reconstruction is a medical imaging technique for creating a meaningful diagnostic image. A computer can take the input from the CT machine, run it through the formula, and return a set of images for a physician to examine. Here, we only consider the simplest mathematical modelm
i
n
X
≥
0
Y
-
A
X
2 to pursue the structural image. A simulated thorax phantom with 128
×
128 grids and 0.5 mm pixels, as shown in Figure 2, is used in the following experiments. There are many advantages to using simulated phantoms, including prior knowledge of the pixel values and the ability to control noise. The total attenuation value is approximately 5
×
1
0
7. For this case, the proposed algorithm is called ISRT, which is widely used in the field of CT reconstruction [4].Figure 2
Shepp-Logan phantom with128
×
128 grids: (a) true phantom and (b) noisy projections.
(a)
(b)The system matrix is obtained using the “angle of view” method [10]. From the system matrix, we forward project the phantom on the sinogram with 128 radial bins (bin size of 0.5 mm) and 128 angular views evenly spaced over π. The noisy projections are obtained by the Poisson-random formula y
i
=
P
o
i
s
s
o
n
(
y
i
∗
), where y
i
∗ is the noise-free projection.Figure3 shows the reconstructions with 50 iterations for every algorithm. As shown, the CG’s image suffers from serious noise artifacts. The PL algorithm provides a reconstructed image that is slightly blurred compared with that of ISRT. The proposed ISRT provides the sharpest edge and clearest interior region.Figure 3
Reconstructed images with 50 iterations by (a) the proposed ISRT, (b) CG, and (c) PL.
(a)
(b)
(c)Figure4 quantitatively compares the histograms of the reconstructed images. The CG’s result still shows a clear turbulence around the true value. PL and ISRT provide two comparable curves, and neither is significantly better than the other.Figure 4
Intensity of the true and reconstructed images: (a) horizontal central line and (b) vertical central line.
(a)
(b)We also present a comparison of the MSE curves andχ with increasing iteration number in Figure 5. As shown, CG exhibits an unsteady reconstruction property in which more iterations may lead to a worse result. Both PL and ISRT overcome this shortcoming so that more iterations always provide considerably better results. In addition, ISRT will provide a more rapid convergence than PL.Figure 5
M
S
E and χ versus iteration number: (a) M
S
E and (b) χ.
(a)
(b)


## 6. Conclusion

We present a special application of a surrogate-based method for solving box-constrained LS problems. A new update rule is developed to iteratively update variables, and this rule exhibits desirable properties, including monotonic decrease of the cost function, self-constraining in the feasible region, and no need to impose a step size. This algorithm covers many existing algorithms, as well as new ones, as the special examples. Targeting the special case of only below-bounded constraints, we provide a rigorous theoretical global convergence proof. We use the simulated data to evaluate the performance of the algorithm, demonstrating that the proposed algorithm provides a stabler and faster convergence than the PL and CG approaches.In this paper, we only provide a global convergence proof for the below-bounded case, not the box-constrained case; thus, proving the convergence for the latter case will be the focus of future work.



---

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