# An Efficient Method for Convex Constrained Rank Minimization Problems Based on DC Programming

Authors: Wanping Yang; Jinkai Zhao; Fengmin Xu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The constrained rank minimization problem has various applications in many fields including machine learning, control, and signal processing. In this paper, we consider the convex constrained rank minimization problem. By introducing a new variable and penalizing an equality constraint to objective function, we reformulate the convex objective function with a rank constraint as a difference of convex functions based on the closed-form solutions, which can be reformulated as DC programming. A stepwise linear approximative algorithm is provided for solving the reformulated model. The performance of our method is tested by applying it to affine rank minimization problems and max-cut problems. Numerical results demonstrate that the method is effective and of high recoverability and results on max-cut show that the method is feasible, which provides better lower bounds and lower rank solutions compared with improved approximation algorithm using semidefinite programming, and they are close to the results of the latest researches.


---

## Body

## 1. Introduction

In recent decades, with the increase of the system acquisition and data services, the explosion of data poses challenges to storage, transmission, and processing, as well as device design. The burgeoning theory of sparse recovery reveals the outstanding sensing ability on the large scales of high-dimensional data. Recently, a new theory called compressed sensing (or compressive sampling (CS)) has appeared. And this theory is praised as a big deal in signal processing area. This approach can acquire the signals while compressing data properly. Its sampling frequency is lower than that of Nyquist, which makes the collections of high resolution signal possible. One noticeable merit of this approach is its ability to combine traditional data collection and data compression into a union when facing sparse representation signals. This merit means that sampling frequency of signals, time and computational cost consuming on data processing, and expense for data storage and transmission are all greatly reduced. Thus this approach leads signal processing to a new age. However, same with CS, Low Rank Matrix theory can also solve signal reconstruction problems that are related to theory and practice. It is also carried out through exerting dimension reduction treatments on unknown multidimension signals validly. And then use sparse concepts in a more broad sense to the problems. That is to say, reconstruct the original high-dimensional signals array completely or approximately by a small number of values measured by sparse (or reducing dimensions) sampling. Considering that both of these two theories are closely related to intrinsic sparsity of signals, we generally view them as sparse recovery problems. The main goal of sparse recovery problems is using less linear measurement data to recovery high-dimensional sparse signal. Its core mainly includes three steps: signal sparse representation, linear dimension reduction measure, and the nonlinear reconstruction of signal. The three parts for the sparse recovery theory are closely related; each part will have a significant impact on the quality of signal reconstruction. However, the nonlinear reconstruction of signal can be seen as a special case of rank minimization problem to be solved.Constrained rank minimization problems have attracted a great deal of attention over recent years, which aim to find sparse solutions of a system. The problems are widely applicable in a variety of fields including machine learning [1, 2], control [3, 4], Euclidean embedding [5], image processing [6], and finance [7, 8], to name but a few. There are some typical examples of constrained rank minimization problems; for example, one is affine matrix rank minimization problems [5, 9, 10](1)
m
i
n
1
2
A
X
-
b
F
2
s
.
t
.
rank
⁡
X
≤
k
,where X
∈
R
m
×
n is the decision variable and the linear map A
:
R
m
×
n
→
R
p and vector b
∈
R
p are known. Equation (1) includes matrix completion, based on affine rank minimization problem and compressed sensing. Matrix completion is widely used in collaborative filtering [11, 12], system identification [13, 14], sensor network [15, 16], image processing [17, 18], sparse channel estimation [19, 20], spectrum sensing [21, 22], and multimedia coding [23, 24]. Another is max-cut problems [25] in combinatorial optimization(2)
m
i
n
Tr
⁡
L
X
s
.
t
.
diag
⁡
X
=
e
rank
⁡
X
=
1
X
⪰
0
,where X
=
x
x
T, x
=
(
x
1
,
…
,
x
n
)
T, L
=
(
1
/
4
)
(
diag
⁡
(
W
e
)
-
W
), e
=
(
1,1
,
…
,
1
)
T, W is the weight matrix, and so on. As a matter of fact, these problems can be written as a common convex constrained rank minimization model; that is, (3)
m
i
n
f
X
s
.
t
.
rank
⁡
X
≤
k
X
∈
Γ
∩
Ω
,where f
:
R
m
×
n
→
R is a continuously differentiable convex function, rank
⁡
(
X
) is the rank of a matrix X, k
≥
0 is a given integer, Γ is a closed convex set, and Ω is a closed unitarily invariant convex set. According to [26], a set C
⊆
R
n
×
n is a unitarily invariant set if(4)
U
X
V
:
U
∈
U
m
,
V
∈
U
n
,
X
∈
C
=
C
,where U
n denotes the set of all unitary matrices in R
n
×
n. Various types of methods have been proposed to solve (3) or its special cases with different types of Γ
∩
Ω.WhenΓ
∩
Ω is R
m
×
n, (3) is the unconstrained rank minimization problem. One approach is to solve this case by convex relaxation of rank
⁡
(
·
). Cai et al. [9] proposed singular value thresholding (SVT) algorithm; Ma et al. [10] proposed fixed point continuation (FPC) algorithm. Another approach is to solve rank constraint directly. Jain et al. [27] proposed simple and efficient algorithm Singular Value Projection (SVP) based on the projected gradient algorithm. Haldar and Hernando [28] proposed an alternating least squares approach. Keshavan et al. [29] proposed an algorithm based on optimization over a Grassmann manifold.WhenΓ
∩
Ω is a convex set, (3) is a convex constrained rank minimization problem. One approach is to solve it by replacing rank constraint with other constraints, such as trace constraint and norm constraint. Based on this, several heuristic methods have been proposed in literature; for example, Nesterov et al. [30] proposed interior-point polynomial algorithms; Weimin [31] proposed an adaptive seminuclear norm regularization approach; in addition, semidefinite programming (SDP) has also been applied; see [32]. Another approach is to solve rank constraint directly. Grigoriadis and Beran [33] proposed alternating projection algorithms for linear matrix inequalities problems. Mohseni et al. [24] proposed penalty decomposition (PD) method based on penalty function. Gao and Sun [34] recently proposed the majorized penalty approach. Burer et al. [35] proposed nonlinear programming (NLP) reformulation approach.In this paper, we consider problem (3). We introduce a new variable and penalize equality constraints to objective function. By this technique, we can obtain closed-form solution and use it to reformulate the objective function in problem (3) as a difference of convex functions; then (3) is reformulated as DC programming which can be solved via linear approximation method. A function is called DC if it can be represented as the difference of two convex functions. Mathematical programming problems dealing with DC functions are called DC programming problems. Our method is different from original PD [26]; for one thing we penalize equality constraints to objective function and keep other constraints, while PD penalizes all except rank constraint to objective function, including equality and inequality constraints; for another each subproblem of PD is approximately solved by a block coordinate descent method, while we approximate problem (3) to a convex optimization to be solved finally. Compared with PD method, our method uses the closed-form solutions to remove rank constraint, while PD needs to consider rank constraint in each iteration. It is known that rank constraint is a nonconvex and discontinuous constraint, which is not easy to deal. And due to the remaining convex constraints, the formulated programming can be solved by solving a series of convex programming. In our method, we use linear function instead of the last convex function to formulate the programming to convex programming. We test the performance of our method by applying it to affine rank minimization problems and max-cut problems. Numerical results show that the algorithm is feasible and effective.The rest of this paper is organized as follows. In Section1.1, we introduce the notation that is used throughout the paper. Our method is shown in Section 2. In Section 3, we apply it to some practical problems to test its performance.


### 1.1. Notations

In this paper, the symbolR
n denotes the n-dimensional Euclidean space, and the set of all m
×
n matrices with real entries is denoted by R
m
×
n. Given matrices X and Y in R
m
×
n, the standard inner product is defined by 〈
X
,
Y
〉
≔
Tr
⁡
(
X
Y
T
), where Tr
⁡
(
·
) denotes the trace of a matrix. The Frobenius norm of a real matrix X is defined as X
F
≔
Tr
⁡
(
X
X
T
). X
∗ denotes the nuclear norm of X, that is, the sum of singular values of X. The rank of a matrix X is denoted by rank
⁡
(
X
). We use I
q to denote the identity matrix, whose dimension is q. Throughout this paper, we always assume that the singular values are arranged in nonincreasing order, that is, σ
1
≥
σ
2
≥
⋯
≥
σ
r
>
0
=
σ
r
+
1
=
⋯
=
σ
min
⁡
{
m
,
n
}. ∂
f denotes the subdifferential of the function f. We denote the adjoint of A by A
∗. A
A
T
(
k
) is used to denote Ky Fan k-norms. Let P
Ω be the projection onto the closed convex set Ω.


## 1.1. Notations

In this paper, the symbolR
n denotes the n-dimensional Euclidean space, and the set of all m
×
n matrices with real entries is denoted by R
m
×
n. Given matrices X and Y in R
m
×
n, the standard inner product is defined by 〈
X
,
Y
〉
≔
Tr
⁡
(
X
Y
T
), where Tr
⁡
(
·
) denotes the trace of a matrix. The Frobenius norm of a real matrix X is defined as X
F
≔
Tr
⁡
(
X
X
T
). X
∗ denotes the nuclear norm of X, that is, the sum of singular values of X. The rank of a matrix X is denoted by rank
⁡
(
X
). We use I
q to denote the identity matrix, whose dimension is q. Throughout this paper, we always assume that the singular values are arranged in nonincreasing order, that is, σ
1
≥
σ
2
≥
⋯
≥
σ
r
>
0
=
σ
r
+
1
=
⋯
=
σ
min
⁡
{
m
,
n
}. ∂
f denotes the subdifferential of the function f. We denote the adjoint of A by A
∗. A
A
T
(
k
) is used to denote Ky Fan k-norms. Let P
Ω be the projection onto the closed convex set Ω.


## 2. An Efficient Algorithm Based on DC Programming

In this section, we first reformulated problem (3) as a DC programming problem with the convex constrained set, which can be solved via stepwise linear approximation method. Then, the convergence is provided.


### 2.1. Model Transformation

LetX
=
Y. Problem (3) is rewritten as(5)
m
i
n
f
X
s
.
t
.
X
=
Y
,
rank
⁡
Y
≤
k
,
X
∈
Γ
∩
Ω
.Adopting penalty function method to penalize X
=
Y to objective function and choosing proper penalty parameter ρ, then (5) can be reformulated as (6)
m
i
n
f
X
+
ρ
2
X
-
Y
F
2
s
.
t
.
rank
⁡
Y
≤
k
,
X
∈
Γ
∩
Ω
.Solving Y with fixed X
∈
Γ
∩
Ω, (6) can be treated as(7)
m
i
n
X
f
X
+
m
i
n
rank
⁡
Y
≤
k
ρ
2
X
-
Y
F
2
s
.
t
.
X
∈
Γ
∩
Ω
.In fact, given the singular value decomposition, X
=
U
S
V
T, and singular values of X: σ
j
X,  j
=
1
,
…
,
min
⁡
{
m
,
n
}, where U
∈
R
m
×
mand V
∈
R
n
×
n are orthogonal matrices, the minimal optimal problem with respect to Y with X fixed in (7), that is, min
rank
⁡
(
Y
)
≤
k
⁡
(
ρ
/
2
)
X
-
Y
F
2, has closed-form solution Y
=
U
C
V
T with(8)
C
=
D
0
0
0
,where D
=
diag
⁡
(
σ
1
,
σ
2
,
σ
3
,
…
,
σ
k
)
∈
R
k
×
k. Thus, by the closed-form solution and the unitary invariance of F-norm, we have (9)
min
rank
⁡
Y
≤
k
⁡
ρ
2
X
-
Y
F
2
=
min
rank
⁡
Y
≤
k
⁡
ρ
2
U
T
X
V
-
U
T
Y
V
F
2
=
ρ
2
S
-
C
F
2
=
ρ
2
X
F
2
-
ρ
2
∑
i
=
1
k
σ
i
2
X
.Substituting the above into problem (7), problem (3) is reformulated as (10)
m
i
n
f
X
+
ρ
2
X
F
2
-
ρ
2
∑
i
=
1
k
σ
i
2
X
s
.
t
.
X
∈
Γ
∩
Ω
.Clearly,f
(
X
)
+
(
ρ
/
2
)
X
F
2 is a convex function about X. Define function g
:
R
m
×
n
→
R, g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
), 1
≤
k
≤
min
⁡
{
m
,
n
}. Next, we will study the property of g
(
X
). Problem (3) would be reformulated as a DC programming problem with convex constraints if g
(
X
) is a convex function, which can be solved via stepwise linear approximation method. In fact, it does hold, and g
(
X
) equals a trace function.Theorem 1.
Let(11)
W
=
I
k
O
O
O
∈
R
m
×
m
,where O denotes zero matrix whose dimension is determined depending on context; then, for any X
=
U
S
V
T
∈
R
m
×
n, the following hold:(1)
g
(
X
) is a convex function;
(2)
g
(
X
)
=
T
r
(
X
T
U
W
U
T
X
).Proof.
(1) Since{
σ
j
2
(
X
)
j
=
1
,
…
,
min
⁡
(
m
,
n
)
} are eigenvalues of X
X
T, then g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
)
=
X
X
T
(
k
), where ·
(
k
) is Ky Fan k-norms (defined in [36]). Clearly, ·
(
k
) is a convex function. Let φ
(
·
)
=
·
k
,
1
≤
k
≤
min
⁡
{
m
,
n
}. For ∀
A
,
B
∈
R
m
×
n and  α
∈
[
0,1
] we have(12)
g
α
A
+
1
-
α
B
=
φ
α
A
+
1
-
α
B
α
A
+
1
-
α
B
T
=
φ
α
2
A
A
T
+
α
1
-
α
A
B
T
+
α
1
-
α
B
A
T
+
1
-
α
2
B
B
T
≤
α
2
φ
A
A
T
+
α
1
-
α
φ
A
B
T
+
α
1
-
α
φ
B
A
T
+
1
-
α
2
φ
B
B
T
.Note that B
A
T
=
(
A
B
T
)
T; thus they have the same nonzero singular values, and then φ
(
B
A
T
)
=
φ
(
A
B
T
) holds. Together with (12), we have (13)
φ
α
A
+
1
-
α
B
α
A
+
1
-
α
B
T
≤
α
2
φ
A
A
T
+
2
α
1
-
α
φ
A
B
T
+
1
-
α
2
φ
B
B
T
.It is known that for arbitrary matrices A
,
B
∈
R
m
×
m, (14)
σ
j
A
∗
B
≤
1
2
σ
j
A
∗
A
+
B
∗
B
,
1
≤
j
≤
m
,
s
e
e
T
h
e
o
r
e
m
I
X
.
4.2
i
n
36
.Further, (15)
∑
j
=
1
k
σ
j
A
T
B
≤
∑
j
=
1
k
1
2
σ
j
A
∗
A
+
B
∗
B
,
1
≤
k
≤
n
.That is, (16)
φ
A
B
T
≤
1
2
φ
A
A
T
+
B
B
T
,where A
A
T and B
B
T are symmetric matrices. It is also known that, for arbitrary Hermitian matrices A
A
T
,
B
B
T
∈
R
m
×
m, k
=
1,2
,
…
,
m, the following holds (see Exercise II.1.15 in [24]) (17)
∑
j
=
1
k
σ
j
A
A
T
+
B
B
T
≤
∑
j
=
1
k
σ
j
A
A
T
+
∑
j
=
1
k
σ
j
B
B
T
.Thus (18)
φ
A
T
A
+
B
T
B
=
φ
A
A
T
+
B
B
T
≤
φ
A
A
T
+
φ
B
B
T
.Together with (13) and (16), (18) yields (19)
g
α
A
+
1
-
α
B
≤
α
φ
A
A
T
+
1
-
α
φ
B
B
T
=
α
g
A
+
1
-
α
g
B
.Hence, g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
) is a convex function.
(2) LetΘ
=
diag
⁡
(
σ
1
,
σ
2
,
σ
3
,
…
,
σ
r
) and r is the rank of X, k
≤
r
(20)
S
-
=
Θ
O
O
Oand we have (21)
T
r
X
T
U
W
U
T
X
=
T
r
V
S
-
T
U
T
U
W
U
T
U
S
-
V
T
=
T
r
V
S
-
T
W
S
-
V
T
=
T
r
S
-
T
W
S
-
.Moreover(22)
S
-
T
W
S
-
=
Θ
O
O
O
T
I
k
O
O
O
Θ
O
O
O
=
Θ
-
O
O
O
,where Θ
-
=
diag
⁡
(
σ
1
2
,
σ
2
2
,
σ
3
2
,
…
,
σ
k
2
).
Thus,T
r
(
S
-
T
W
S
-
)
=
∑
i
=
1
k
σ
i
2
(
X
)
=
g
(
X
) holds.According to above theorem, (10) can be rewritten as(23)
m
i
n
f
X
+
ρ
2
X
F
2
-
ρ
2
T
r
X
T
U
W
U
T
X
s
.
t
.
X
∈
Γ
∩
Ω
.Since the objective function of (23) is a difference of convex functions, the constrained set Γ
∩
Ω is a closed convex constraint, so (23) is DC programming. Next, we will use stepwise linear approximate method to solve it.


### 2.2. Stepwise Linear Approximative Algorithm

In this subsection, we use stepwise linear approximation method [37] to solve reformulated model (23). Let F
(
X
) represent objective function of model (23); that is, (24)
F
X
=
f
X
+
ρ
2
X
F
2
-
ρ
2
T
r
X
T
U
W
U
T
X
,and the linear approximation function of F
(
X
) can be defined as (25)
H
X
,
X
i
=
f
X
+
ρ
2
X
F
2
-
ρ
2
g
X
i
+
∂
g
X
i
,
X
-
X
i
,where (
∂
g
X
/
∂
X
)
(
∂
T
r
X
T
U
W
U
T
X
/
∂
X
)
=
2
U
W
U
T
X (see 5.2.2 in [36]).Stepwise linear approximative algorithm for solving model (23) can be described as follows.Given the initial matrixX
0
∈
Γ
∩
Ω, penalty parameter ρ
>
0, rank k
>
0, τ
>
1, set i
=
0
.Step 1.
Compute the single value decomposition ofX
i, U
i
Σ
i
V
i
=
X
i
.Step 2.
Update(26)
X
i
+
1
∈
a
r
g
min
X
∈
Γ
∩
Ω
⁡
H
X
,
X
iStep 3.
If stopping criterion is satisfied, stop; else seti
←
i
+
1 and ρ
←
(
1
+
τ
)
ρ and go to Step 1.Remark
(1)
Subproblem (26) can be done by solving a series of convex subproblems. The efficient method for convex programming is multiple and it has perfect theory guarantees.
(a)
WhenΓ
∩
Ω is R
m
×
n, closed-form solutions can be obtained; for example, when solving model (1), X
i
+
1
=
(
A
∗
A
+
ρ
I
)
-
1
(
ρ
U
i
W
(
U
i
)
T
X
i
) are obtained by the first-order optimality conditions (27)
0
∈
∂
f
X
+
ρ
X
-
ρ
U
i
W
U
i
T
X
i
.
(b)
WhenΓ
∩
Ω is a general convex set, the closed-form solution is hard to get, but subproblem (26) can be solved by some convex optimization software programs (such as CVX and CPLEX). According to different closed convex constraint sets, we can pick and choose suitable convex programming algorithm to combine with our method. The toolbox listed in this paper is just an alternative choice.
(2)
The stop rule that we adopted here isX
i
+
1
-
X
i
F
≤
ε with some ε small enough.


### 2.3. Convergence Analysis

In this subsection, we analyze the convergence properties of the sequence{
X
i
} generated by above algorithm. Before we prove the main convergence result, we give the following lemma which is analogous to Lemma 1 in [10].Lemma 2.
For anyY
1 and Y
2
∈
R
m
×
n
(28)
Y
1
¯
-
Y
2
¯
F
≤
Y
1
-
Y
2
F
,where Y
i
¯ is the pre-k singular values approximate of Y
i  (i
=
1,2).Proof.
Without loss of generality, we assumem
≤
n and the singular value decomposition Y
1
=
U
1
S
1
V
1
T
Y
2
=
U
2
S
2
V
2
T with(29)
S
1
=
diag
⁡
σ
0
0
0
∈
R
m
×
n
S
2
=
diag
⁡
γ
0
0
0
∈
R
m
×
n
,where σ
=
(
σ
1
,
σ
2
,
…
,
σ
s
),  σ
1
≥
⋯
≥
σ
s
>
0 and γ
=
(
γ
1
,
γ
2
,
…
,
γ
t
),   γ
1
≥
⋯
≥
γ
t
>
0. Y
i
¯ can be written as Y
i
¯
=
U
i
W
S
i
V
i
T, (i
=
1,2). Note that U
1, V
1, U
2, V
2 are orthogonal matrices; we have(30)
Y
1
-
Y
2
F
2
-
Y
1
¯
-
Y
2
¯
F
2
=
Tr
⁡
Y
1
-
Y
2
T
Y
1
-
Y
2
-
Tr
⁡
Y
1
¯
-
Y
2
¯
T
Y
1
¯
-
Y
2
¯
=
Tr
⁡
Y
1
T
Y
1
-
Y
1
¯
T
Y
1
¯
+
Y
2
T
Y
2
-
Y
2
¯
T
Y
2
¯
-
2
Tr
⁡
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
=
∑
i
=
1
s
σ
i
2
-
∑
i
=
1
k
σ
i
2
+
∑
i
=
1
t
γ
i
2
-
∑
i
=
1
k
γ
i
2
-
2
Tr
⁡
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
.And we note that(31)
T
r
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
=
T
r
Y
1
-
Y
1
¯
T
Y
2
-
Y
2
¯
+
Y
1
-
Y
1
¯
T
Y
2
¯
+
Y
1
¯
T
Y
2
-
Y
2
¯
=
T
r
V
1
S
1
-
W
S
1
T
U
1
T
U
2
S
2
-
W
S
2
V
2
T
+
V
1
S
1
-
W
S
1
T
U
1
T
U
2
W
S
2
V
2
T
+
V
1
W
S
T
U
1
T
U
2
S
2
-
W
S
2
V
2
T
=
T
r
S
1
-
W
S
1
T
U
S
2
-
W
S
2
V
T
+
S
1
-
W
S
1
T
U
W
S
2
V
T
+
W
S
1
T
U
S
2
-
W
S
2
V
T
,where U
=
U
1
T
U
2
,
V
=
V
1
T
V
2 are also orthogonal matrices.
Next we estimate an upper bound forT
r
(
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
). According to Theorem 7.4.9 in [38], we get that an orthogonal matrix U is a maximizing matrix for the problem m
a
x
{
T
r
(
A
U
)
:
U
i
s
o
r
t
h
o
g
o
n
a
l
}, if and only if A
U is positive semidefinite matrix. Thus, T
r
(
(
S
1
-
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T
), T
r
(
(
S
1
-
W
S
1
)
T
U
W
S
2
V
T
), and T
r
(
(
W
S
)
T
U
(
S
2
-
W
S
2
)
V
T
) achieve their maximum, if and only if (
S
1
-
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T,  (
S
1
-
W
S
1
)
T
U
W
S
2
V
T,  and (
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T are all positive semidefinite. It is known that when A
B is positive semidefinite, (32)
T
r
A
B
=
∑
i
σ
i
A
B
≤
∑
i
σ
i
A
σ
i
B
.Applying (32) to above three terms, we have (33)
T
r
S
1
-
W
S
1
T
U
S
2
-
W
S
2
V
T
≤
∑
i
σ
i
S
1
-
W
S
1
σ
i
S
2
-
W
S
2
T
r
S
1
-
W
S
1
T
U
W
S
2
V
T
≤
∑
i
σ
i
S
1
-
W
S
1
σ
i
W
S
2
T
r
W
S
1
T
U
S
2
-
W
S
2
V
T
≤
∑
i
σ
i
W
S
1
σ
i
S
2
-
W
S
2
.Hence, without loss of generality, assuming k
≤
s
≤
t, we have(34)
Y
1
-
Y
2
F
2
-
Y
1
¯
-
Y
2
¯
F
2
≥
∑
i
=
1
s
σ
i
2
-
∑
i
=
1
k
σ
i
2
+
∑
i
=
1
t
γ
i
2
-
∑
i
=
1
k
γ
i
2
-
2
∑
i
=
k
+
1
s
σ
i
γ
i
=
∑
i
=
k
+
1
s
σ
i
-
γ
i
2
+
∑
i
=
s
+
1
t
γ
i
2
≥
0which implies that Y
1
¯
-
Y
2
¯
F
≤
Y
1
-
Y
2
F holds.Next, we claim the value of iterations of stepwise linear approximative algorithm is nonincreasing.Theorem 3.
Let{
X
i
} be the iterative sequence generated by stepwise linear approximative algorithm. Then {
F
(
X
i
)
} is monotonically nonincreasing.Proof.
SinceX
i
+
1
∈
arg
⁡
min
⁡
{
H
(
X
,
X
i
)
∣
X
∈
Γ
∩
Ω
}, we have(35)
H
X
i
+
1
,
X
i
≤
H
X
i
,
X
i
=
f
X
i
+
ρ
2
X
i
F
2
-
ρ
2
g
X
i
=
F
X
i
,and ∂
g
(
X
) is the subdifferential of g
(
X
), which implies (36)
∀
X
Y
,
g
X
≥
g
Y
+
∂
g
Y
,
X
-
Y
.Let X
=
X
i
+
1,  Y
=
X
i in above inequality; we have(37)
g
X
i
+
1
≥
g
X
i
+
∂
g
X
i
,
X
i
+
1
-
X
i
.Therefore,(38)
H
X
i
+
1
,
X
i
=
f
X
i
+
1
+
ρ
2
X
i
+
1
F
2
-
ρ
2
g
X
i
+
∂
g
X
i
,
X
-
X
i
≥
f
X
i
+
1
+
ρ
2
X
i
+
1
F
2
-
ρ
2
g
X
i
+
1
=
F
X
i
+
1
.Consequently, (39)
F
X
i
+
1
≤
H
X
i
+
1
,
X
i
≤
H
X
i
,
X
i
=
F
X
i
.We immediately obtain the conclusion as desired.Next, we show the convergence of iterative sequence generated by stepwise linear approximative algorithm when solving (23). Theorem 4 shows the convergence when solving (23) with Γ
∩
Ω
=
R
m
×
n.Theorem 4.
LetΥ be the set of stable points of problem (23) with no constraints; then the iterative sequence {
X
i
} generated by linear approximative algorithm converges to some X
o
p
t
∈
Υ when ρ
→
∞.Proof.
Assuming thatX
opt
∈
Υ is any optimal solution of problem (23) with Γ
∩
Ω
=
R
m
×
n, according to the first-order optimality conditions, we have (40)
0
=
∂
f
X
opt
+
ρ
X
opt
-
ρ
U
opt
W
U
opt
T
X
opt
.Since X
i
+
1
∈
arg
⁡
min
⁡
f
X
+
(
ρ
/
2
)
X
F
2
-
(
ρ
/
2
)
{
g
(
X
i
)
+
〈
∂
g
(
X
i
)
,
X
-
X
i
〉
}, then(41)
0
=
∂
f
X
i
+
1
+
ρ
X
i
+
1
-
ρ
U
i
W
U
i
T
X
i
.Equation (41) subtracts (40): (42)
X
i
+
1
-
X
opt
=
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
+
U
opt
W
U
opt
T
X
opt
-
U
i
W
U
i
T
X
iin which U
opt
W
(
U
opt
)
T
X
opt
-
U
i
W
(
U
i
)
T
X
i
=
U
opt
W
∑
opt
(
V
opt
)
T
-
U
i
W
∑
i
(
V
opt
)
T
.
By Lemma2, we have (43)
U
opt
W
U
opt
T
X
opt
-
U
i
W
U
i
T
X
i
F
=
U
opt
W
∑
opt
V
opt
T
-
U
i
W
∑
i
V
i
T
F
≤
X
i
-
X
opt
F
.Norming both sides of (42), by the triangle inequalities property of norms, we have (44)
X
i
+
1
-
X
opt
F
=
l
i
m
ρ
→
∞
X
i
+
1
-
X
opt
F
≤
lim
ρ
→
∞
⁡
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
F
+
X
i
-
X
opt
F
=
lim
ρ
→
∞
⁡
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
F
+
lim
ρ
→
∞
⁡
X
i
-
X
opt
F
=
X
i
-
X
opt
F
,which implies that the sequence {
X
i
-
X
opt
F
} is monotonically nonincreasing when ρ
→
∞. We then observe that {
X
i
-
X
opt
F
} is bounded, and hence l
i
m
i
→
∞
X
i
-
X
opt
F exists, denoted by X
¯
-
X
opt
F, where X
¯ is the limited point of {
X
i
}. Let i
→
∞ in (41); by the continuity, we get 0
=
∂
f
(
X
¯
)
+
ρ
X
¯
-
ρ
U
¯
W
U
¯
T
X
¯, which implies that X
¯
∈
Υ. By setting X
¯
=
X
opt, we have (45)
lim
i
→
∞
⁡
X
i
-
X
opt
F
=
X
¯
-
X
opt
F
=
0
,which completes the proof.Now, we are ready to give the convergence of iterative sequence when solving (23) with general closed convex Γ
∩
Ω.Theorem 5.
Letξ be the set of optimal solutions of problem (23); then the iterative sequence {
Y
i
} generated by stepwise linear approximative algorithm converges to some Y
o
p
t
∈
ξ when ρ
→
∞.Proof.
Generally speaking, subproblem (26) can be solved by projecting a solution obtained via solving a problem with no constraints onto a closed convex set; that is, Y
i
+
1
=
P
Γ
∩
Ω
(
X
i
+
1
), where X
i
+
1
∈
arg
⁡
min
⁡
{
H
(
X
,
X
i
)
}. According to Theorem 4, the sequence {
X
i
-
X
opt
F
} is monotonically nonincreasing and convergent when ρ
→
∞; together with the well known fact that P is a nonexpansive operator, we have (46)
Y
i
-
Y
opt
F
=
P
Γ
∩
Ω
X
i
-
P
Γ
∩
Ω
X
opt
F
≤
X
i
-
X
opt
F
,where Y
opt is the projection of X
opt onto Γ
∩
Ω. Clearly, Y
opt
∈
ξ.
Hence,0
≤
lim
i
→
∞
⁡
Y
i
-
Y
opt
F
≤
lim
i
→
∞
⁡
X
i
-
X
opt
F
=
0, which implies that Y
i converge to Y
opt. We immediately obtain the conclusion.


## 2.1. Model Transformation

LetX
=
Y. Problem (3) is rewritten as(5)
m
i
n
f
X
s
.
t
.
X
=
Y
,
rank
⁡
Y
≤
k
,
X
∈
Γ
∩
Ω
.Adopting penalty function method to penalize X
=
Y to objective function and choosing proper penalty parameter ρ, then (5) can be reformulated as (6)
m
i
n
f
X
+
ρ
2
X
-
Y
F
2
s
.
t
.
rank
⁡
Y
≤
k
,
X
∈
Γ
∩
Ω
.Solving Y with fixed X
∈
Γ
∩
Ω, (6) can be treated as(7)
m
i
n
X
f
X
+
m
i
n
rank
⁡
Y
≤
k
ρ
2
X
-
Y
F
2
s
.
t
.
X
∈
Γ
∩
Ω
.In fact, given the singular value decomposition, X
=
U
S
V
T, and singular values of X: σ
j
X,  j
=
1
,
…
,
min
⁡
{
m
,
n
}, where U
∈
R
m
×
mand V
∈
R
n
×
n are orthogonal matrices, the minimal optimal problem with respect to Y with X fixed in (7), that is, min
rank
⁡
(
Y
)
≤
k
⁡
(
ρ
/
2
)
X
-
Y
F
2, has closed-form solution Y
=
U
C
V
T with(8)
C
=
D
0
0
0
,where D
=
diag
⁡
(
σ
1
,
σ
2
,
σ
3
,
…
,
σ
k
)
∈
R
k
×
k. Thus, by the closed-form solution and the unitary invariance of F-norm, we have (9)
min
rank
⁡
Y
≤
k
⁡
ρ
2
X
-
Y
F
2
=
min
rank
⁡
Y
≤
k
⁡
ρ
2
U
T
X
V
-
U
T
Y
V
F
2
=
ρ
2
S
-
C
F
2
=
ρ
2
X
F
2
-
ρ
2
∑
i
=
1
k
σ
i
2
X
.Substituting the above into problem (7), problem (3) is reformulated as (10)
m
i
n
f
X
+
ρ
2
X
F
2
-
ρ
2
∑
i
=
1
k
σ
i
2
X
s
.
t
.
X
∈
Γ
∩
Ω
.Clearly,f
(
X
)
+
(
ρ
/
2
)
X
F
2 is a convex function about X. Define function g
:
R
m
×
n
→
R, g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
), 1
≤
k
≤
min
⁡
{
m
,
n
}. Next, we will study the property of g
(
X
). Problem (3) would be reformulated as a DC programming problem with convex constraints if g
(
X
) is a convex function, which can be solved via stepwise linear approximation method. In fact, it does hold, and g
(
X
) equals a trace function.Theorem 1.
Let(11)
W
=
I
k
O
O
O
∈
R
m
×
m
,where O denotes zero matrix whose dimension is determined depending on context; then, for any X
=
U
S
V
T
∈
R
m
×
n, the following hold:(1)
g
(
X
) is a convex function;
(2)
g
(
X
)
=
T
r
(
X
T
U
W
U
T
X
).Proof.
(1) Since{
σ
j
2
(
X
)
j
=
1
,
…
,
min
⁡
(
m
,
n
)
} are eigenvalues of X
X
T, then g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
)
=
X
X
T
(
k
), where ·
(
k
) is Ky Fan k-norms (defined in [36]). Clearly, ·
(
k
) is a convex function. Let φ
(
·
)
=
·
k
,
1
≤
k
≤
min
⁡
{
m
,
n
}. For ∀
A
,
B
∈
R
m
×
n and  α
∈
[
0,1
] we have(12)
g
α
A
+
1
-
α
B
=
φ
α
A
+
1
-
α
B
α
A
+
1
-
α
B
T
=
φ
α
2
A
A
T
+
α
1
-
α
A
B
T
+
α
1
-
α
B
A
T
+
1
-
α
2
B
B
T
≤
α
2
φ
A
A
T
+
α
1
-
α
φ
A
B
T
+
α
1
-
α
φ
B
A
T
+
1
-
α
2
φ
B
B
T
.Note that B
A
T
=
(
A
B
T
)
T; thus they have the same nonzero singular values, and then φ
(
B
A
T
)
=
φ
(
A
B
T
) holds. Together with (12), we have (13)
φ
α
A
+
1
-
α
B
α
A
+
1
-
α
B
T
≤
α
2
φ
A
A
T
+
2
α
1
-
α
φ
A
B
T
+
1
-
α
2
φ
B
B
T
.It is known that for arbitrary matrices A
,
B
∈
R
m
×
m, (14)
σ
j
A
∗
B
≤
1
2
σ
j
A
∗
A
+
B
∗
B
,
1
≤
j
≤
m
,
s
e
e
T
h
e
o
r
e
m
I
X
.
4.2
i
n
36
.Further, (15)
∑
j
=
1
k
σ
j
A
T
B
≤
∑
j
=
1
k
1
2
σ
j
A
∗
A
+
B
∗
B
,
1
≤
k
≤
n
.That is, (16)
φ
A
B
T
≤
1
2
φ
A
A
T
+
B
B
T
,where A
A
T and B
B
T are symmetric matrices. It is also known that, for arbitrary Hermitian matrices A
A
T
,
B
B
T
∈
R
m
×
m, k
=
1,2
,
…
,
m, the following holds (see Exercise II.1.15 in [24]) (17)
∑
j
=
1
k
σ
j
A
A
T
+
B
B
T
≤
∑
j
=
1
k
σ
j
A
A
T
+
∑
j
=
1
k
σ
j
B
B
T
.Thus (18)
φ
A
T
A
+
B
T
B
=
φ
A
A
T
+
B
B
T
≤
φ
A
A
T
+
φ
B
B
T
.Together with (13) and (16), (18) yields (19)
g
α
A
+
1
-
α
B
≤
α
φ
A
A
T
+
1
-
α
φ
B
B
T
=
α
g
A
+
1
-
α
g
B
.Hence, g
(
X
)
=
∑
i
=
1
k
σ
i
2
(
X
) is a convex function.
(2) LetΘ
=
diag
⁡
(
σ
1
,
σ
2
,
σ
3
,
…
,
σ
r
) and r is the rank of X, k
≤
r
(20)
S
-
=
Θ
O
O
Oand we have (21)
T
r
X
T
U
W
U
T
X
=
T
r
V
S
-
T
U
T
U
W
U
T
U
S
-
V
T
=
T
r
V
S
-
T
W
S
-
V
T
=
T
r
S
-
T
W
S
-
.Moreover(22)
S
-
T
W
S
-
=
Θ
O
O
O
T
I
k
O
O
O
Θ
O
O
O
=
Θ
-
O
O
O
,where Θ
-
=
diag
⁡
(
σ
1
2
,
σ
2
2
,
σ
3
2
,
…
,
σ
k
2
).
Thus,T
r
(
S
-
T
W
S
-
)
=
∑
i
=
1
k
σ
i
2
(
X
)
=
g
(
X
) holds.According to above theorem, (10) can be rewritten as(23)
m
i
n
f
X
+
ρ
2
X
F
2
-
ρ
2
T
r
X
T
U
W
U
T
X
s
.
t
.
X
∈
Γ
∩
Ω
.Since the objective function of (23) is a difference of convex functions, the constrained set Γ
∩
Ω is a closed convex constraint, so (23) is DC programming. Next, we will use stepwise linear approximate method to solve it.


## 2.2. Stepwise Linear Approximative Algorithm

In this subsection, we use stepwise linear approximation method [37] to solve reformulated model (23). Let F
(
X
) represent objective function of model (23); that is, (24)
F
X
=
f
X
+
ρ
2
X
F
2
-
ρ
2
T
r
X
T
U
W
U
T
X
,and the linear approximation function of F
(
X
) can be defined as (25)
H
X
,
X
i
=
f
X
+
ρ
2
X
F
2
-
ρ
2
g
X
i
+
∂
g
X
i
,
X
-
X
i
,where (
∂
g
X
/
∂
X
)
(
∂
T
r
X
T
U
W
U
T
X
/
∂
X
)
=
2
U
W
U
T
X (see 5.2.2 in [36]).Stepwise linear approximative algorithm for solving model (23) can be described as follows.Given the initial matrixX
0
∈
Γ
∩
Ω, penalty parameter ρ
>
0, rank k
>
0, τ
>
1, set i
=
0
.Step 1.
Compute the single value decomposition ofX
i, U
i
Σ
i
V
i
=
X
i
.Step 2.
Update(26)
X
i
+
1
∈
a
r
g
min
X
∈
Γ
∩
Ω
⁡
H
X
,
X
iStep 3.
If stopping criterion is satisfied, stop; else seti
←
i
+
1 and ρ
←
(
1
+
τ
)
ρ and go to Step 1.Remark
(1)
Subproblem (26) can be done by solving a series of convex subproblems. The efficient method for convex programming is multiple and it has perfect theory guarantees.
(a)
WhenΓ
∩
Ω is R
m
×
n, closed-form solutions can be obtained; for example, when solving model (1), X
i
+
1
=
(
A
∗
A
+
ρ
I
)
-
1
(
ρ
U
i
W
(
U
i
)
T
X
i
) are obtained by the first-order optimality conditions (27)
0
∈
∂
f
X
+
ρ
X
-
ρ
U
i
W
U
i
T
X
i
.
(b)
WhenΓ
∩
Ω is a general convex set, the closed-form solution is hard to get, but subproblem (26) can be solved by some convex optimization software programs (such as CVX and CPLEX). According to different closed convex constraint sets, we can pick and choose suitable convex programming algorithm to combine with our method. The toolbox listed in this paper is just an alternative choice.
(2)
The stop rule that we adopted here isX
i
+
1
-
X
i
F
≤
ε with some ε small enough.


## 2.3. Convergence Analysis

In this subsection, we analyze the convergence properties of the sequence{
X
i
} generated by above algorithm. Before we prove the main convergence result, we give the following lemma which is analogous to Lemma 1 in [10].Lemma 2.
For anyY
1 and Y
2
∈
R
m
×
n
(28)
Y
1
¯
-
Y
2
¯
F
≤
Y
1
-
Y
2
F
,where Y
i
¯ is the pre-k singular values approximate of Y
i  (i
=
1,2).Proof.
Without loss of generality, we assumem
≤
n and the singular value decomposition Y
1
=
U
1
S
1
V
1
T
Y
2
=
U
2
S
2
V
2
T with(29)
S
1
=
diag
⁡
σ
0
0
0
∈
R
m
×
n
S
2
=
diag
⁡
γ
0
0
0
∈
R
m
×
n
,where σ
=
(
σ
1
,
σ
2
,
…
,
σ
s
),  σ
1
≥
⋯
≥
σ
s
>
0 and γ
=
(
γ
1
,
γ
2
,
…
,
γ
t
),   γ
1
≥
⋯
≥
γ
t
>
0. Y
i
¯ can be written as Y
i
¯
=
U
i
W
S
i
V
i
T, (i
=
1,2). Note that U
1, V
1, U
2, V
2 are orthogonal matrices; we have(30)
Y
1
-
Y
2
F
2
-
Y
1
¯
-
Y
2
¯
F
2
=
Tr
⁡
Y
1
-
Y
2
T
Y
1
-
Y
2
-
Tr
⁡
Y
1
¯
-
Y
2
¯
T
Y
1
¯
-
Y
2
¯
=
Tr
⁡
Y
1
T
Y
1
-
Y
1
¯
T
Y
1
¯
+
Y
2
T
Y
2
-
Y
2
¯
T
Y
2
¯
-
2
Tr
⁡
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
=
∑
i
=
1
s
σ
i
2
-
∑
i
=
1
k
σ
i
2
+
∑
i
=
1
t
γ
i
2
-
∑
i
=
1
k
γ
i
2
-
2
Tr
⁡
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
.And we note that(31)
T
r
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
=
T
r
Y
1
-
Y
1
¯
T
Y
2
-
Y
2
¯
+
Y
1
-
Y
1
¯
T
Y
2
¯
+
Y
1
¯
T
Y
2
-
Y
2
¯
=
T
r
V
1
S
1
-
W
S
1
T
U
1
T
U
2
S
2
-
W
S
2
V
2
T
+
V
1
S
1
-
W
S
1
T
U
1
T
U
2
W
S
2
V
2
T
+
V
1
W
S
T
U
1
T
U
2
S
2
-
W
S
2
V
2
T
=
T
r
S
1
-
W
S
1
T
U
S
2
-
W
S
2
V
T
+
S
1
-
W
S
1
T
U
W
S
2
V
T
+
W
S
1
T
U
S
2
-
W
S
2
V
T
,where U
=
U
1
T
U
2
,
V
=
V
1
T
V
2 are also orthogonal matrices.
Next we estimate an upper bound forT
r
(
Y
1
T
Y
2
-
Y
1
¯
T
Y
2
¯
). According to Theorem 7.4.9 in [38], we get that an orthogonal matrix U is a maximizing matrix for the problem m
a
x
{
T
r
(
A
U
)
:
U
i
s
o
r
t
h
o
g
o
n
a
l
}, if and only if A
U is positive semidefinite matrix. Thus, T
r
(
(
S
1
-
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T
), T
r
(
(
S
1
-
W
S
1
)
T
U
W
S
2
V
T
), and T
r
(
(
W
S
)
T
U
(
S
2
-
W
S
2
)
V
T
) achieve their maximum, if and only if (
S
1
-
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T,  (
S
1
-
W
S
1
)
T
U
W
S
2
V
T,  and (
W
S
1
)
T
U
(
S
2
-
W
S
2
)
V
T are all positive semidefinite. It is known that when A
B is positive semidefinite, (32)
T
r
A
B
=
∑
i
σ
i
A
B
≤
∑
i
σ
i
A
σ
i
B
.Applying (32) to above three terms, we have (33)
T
r
S
1
-
W
S
1
T
U
S
2
-
W
S
2
V
T
≤
∑
i
σ
i
S
1
-
W
S
1
σ
i
S
2
-
W
S
2
T
r
S
1
-
W
S
1
T
U
W
S
2
V
T
≤
∑
i
σ
i
S
1
-
W
S
1
σ
i
W
S
2
T
r
W
S
1
T
U
S
2
-
W
S
2
V
T
≤
∑
i
σ
i
W
S
1
σ
i
S
2
-
W
S
2
.Hence, without loss of generality, assuming k
≤
s
≤
t, we have(34)
Y
1
-
Y
2
F
2
-
Y
1
¯
-
Y
2
¯
F
2
≥
∑
i
=
1
s
σ
i
2
-
∑
i
=
1
k
σ
i
2
+
∑
i
=
1
t
γ
i
2
-
∑
i
=
1
k
γ
i
2
-
2
∑
i
=
k
+
1
s
σ
i
γ
i
=
∑
i
=
k
+
1
s
σ
i
-
γ
i
2
+
∑
i
=
s
+
1
t
γ
i
2
≥
0which implies that Y
1
¯
-
Y
2
¯
F
≤
Y
1
-
Y
2
F holds.Next, we claim the value of iterations of stepwise linear approximative algorithm is nonincreasing.Theorem 3.
Let{
X
i
} be the iterative sequence generated by stepwise linear approximative algorithm. Then {
F
(
X
i
)
} is monotonically nonincreasing.Proof.
SinceX
i
+
1
∈
arg
⁡
min
⁡
{
H
(
X
,
X
i
)
∣
X
∈
Γ
∩
Ω
}, we have(35)
H
X
i
+
1
,
X
i
≤
H
X
i
,
X
i
=
f
X
i
+
ρ
2
X
i
F
2
-
ρ
2
g
X
i
=
F
X
i
,and ∂
g
(
X
) is the subdifferential of g
(
X
), which implies (36)
∀
X
Y
,
g
X
≥
g
Y
+
∂
g
Y
,
X
-
Y
.Let X
=
X
i
+
1,  Y
=
X
i in above inequality; we have(37)
g
X
i
+
1
≥
g
X
i
+
∂
g
X
i
,
X
i
+
1
-
X
i
.Therefore,(38)
H
X
i
+
1
,
X
i
=
f
X
i
+
1
+
ρ
2
X
i
+
1
F
2
-
ρ
2
g
X
i
+
∂
g
X
i
,
X
-
X
i
≥
f
X
i
+
1
+
ρ
2
X
i
+
1
F
2
-
ρ
2
g
X
i
+
1
=
F
X
i
+
1
.Consequently, (39)
F
X
i
+
1
≤
H
X
i
+
1
,
X
i
≤
H
X
i
,
X
i
=
F
X
i
.We immediately obtain the conclusion as desired.Next, we show the convergence of iterative sequence generated by stepwise linear approximative algorithm when solving (23). Theorem 4 shows the convergence when solving (23) with Γ
∩
Ω
=
R
m
×
n.Theorem 4.
LetΥ be the set of stable points of problem (23) with no constraints; then the iterative sequence {
X
i
} generated by linear approximative algorithm converges to some X
o
p
t
∈
Υ when ρ
→
∞.Proof.
Assuming thatX
opt
∈
Υ is any optimal solution of problem (23) with Γ
∩
Ω
=
R
m
×
n, according to the first-order optimality conditions, we have (40)
0
=
∂
f
X
opt
+
ρ
X
opt
-
ρ
U
opt
W
U
opt
T
X
opt
.Since X
i
+
1
∈
arg
⁡
min
⁡
f
X
+
(
ρ
/
2
)
X
F
2
-
(
ρ
/
2
)
{
g
(
X
i
)
+
〈
∂
g
(
X
i
)
,
X
-
X
i
〉
}, then(41)
0
=
∂
f
X
i
+
1
+
ρ
X
i
+
1
-
ρ
U
i
W
U
i
T
X
i
.Equation (41) subtracts (40): (42)
X
i
+
1
-
X
opt
=
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
+
U
opt
W
U
opt
T
X
opt
-
U
i
W
U
i
T
X
iin which U
opt
W
(
U
opt
)
T
X
opt
-
U
i
W
(
U
i
)
T
X
i
=
U
opt
W
∑
opt
(
V
opt
)
T
-
U
i
W
∑
i
(
V
opt
)
T
.
By Lemma2, we have (43)
U
opt
W
U
opt
T
X
opt
-
U
i
W
U
i
T
X
i
F
=
U
opt
W
∑
opt
V
opt
T
-
U
i
W
∑
i
V
i
T
F
≤
X
i
-
X
opt
F
.Norming both sides of (42), by the triangle inequalities property of norms, we have (44)
X
i
+
1
-
X
opt
F
=
l
i
m
ρ
→
∞
X
i
+
1
-
X
opt
F
≤
lim
ρ
→
∞
⁡
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
F
+
X
i
-
X
opt
F
=
lim
ρ
→
∞
⁡
1
ρ
∂
f
X
i
+
1
-
∂
f
X
opt
F
+
lim
ρ
→
∞
⁡
X
i
-
X
opt
F
=
X
i
-
X
opt
F
,which implies that the sequence {
X
i
-
X
opt
F
} is monotonically nonincreasing when ρ
→
∞. We then observe that {
X
i
-
X
opt
F
} is bounded, and hence l
i
m
i
→
∞
X
i
-
X
opt
F exists, denoted by X
¯
-
X
opt
F, where X
¯ is the limited point of {
X
i
}. Let i
→
∞ in (41); by the continuity, we get 0
=
∂
f
(
X
¯
)
+
ρ
X
¯
-
ρ
U
¯
W
U
¯
T
X
¯, which implies that X
¯
∈
Υ. By setting X
¯
=
X
opt, we have (45)
lim
i
→
∞
⁡
X
i
-
X
opt
F
=
X
¯
-
X
opt
F
=
0
,which completes the proof.Now, we are ready to give the convergence of iterative sequence when solving (23) with general closed convex Γ
∩
Ω.Theorem 5.
Letξ be the set of optimal solutions of problem (23); then the iterative sequence {
Y
i
} generated by stepwise linear approximative algorithm converges to some Y
o
p
t
∈
ξ when ρ
→
∞.Proof.
Generally speaking, subproblem (26) can be solved by projecting a solution obtained via solving a problem with no constraints onto a closed convex set; that is, Y
i
+
1
=
P
Γ
∩
Ω
(
X
i
+
1
), where X
i
+
1
∈
arg
⁡
min
⁡
{
H
(
X
,
X
i
)
}. According to Theorem 4, the sequence {
X
i
-
X
opt
F
} is monotonically nonincreasing and convergent when ρ
→
∞; together with the well known fact that P is a nonexpansive operator, we have (46)
Y
i
-
Y
opt
F
=
P
Γ
∩
Ω
X
i
-
P
Γ
∩
Ω
X
opt
F
≤
X
i
-
X
opt
F
,where Y
opt is the projection of X
opt onto Γ
∩
Ω. Clearly, Y
opt
∈
ξ.
Hence,0
≤
lim
i
→
∞
⁡
Y
i
-
Y
opt
F
≤
lim
i
→
∞
⁡
X
i
-
X
opt
F
=
0, which implies that Y
i converge to Y
opt. We immediately obtain the conclusion.


## 3. Numerical Results

In this section, we demonstrate the performance of our method proposed in Section2 by applying it to solve affine rank minimization problems (image reconstruction and matrix completion) and max-cut problems. The codes implemented in this section are written in MATLAB and all experiments are performed in MATLAB R2013a running Windows 7.


### 3.1. Image Reconstruction Problems

In this subsection, we apply our method to solve one case of affine rank minimization problems. It can be formulated as(47)
m
i
n
1
2
A
X
-
b
F
2
s
.
t
.
rank
⁡
X
≤
k
X
∈
R
m
×
n
,where X
∈
R
m
×
n is the decision variable and the linear map A
:
R
m
×
n
→
R
p and vector b
∈
R
p are known. Clearly, (47) is a special case of problem (3). Hence, our method can be suitably applied to (47).Recht et al. [39] have demonstrated that when we sample linear maps from a class of probability distributions, then they would obey the Restricted Isometry Property. There are two ingredients for a random linear map to be nearly isometric. First, it must be isometric in expectation. Second, the probability of large distortions of length must be exponentially small. In our numerical experiments, we use four nearly isometric measurements. The first one is the ensemble with independent, identically distributed (i.i.d.) Gaussian entries. The second one has entries sampled from an i.i.d. symmetric Bernoulli distribution, which we call Bernoulli. The third one has zeros in two-thirds of the entries and others sampled from an i.i.d. symmetric Bernoulli distribution, which we call Sparse Bernoulli. The last one has an orthonormal basis, which we call Random Projection.Now we conduct numerical experiments to test the performance of our method for solving (47). To illustrate the scaling of low rank recovery for a particular matrix X, consider the MIT logo presented in [39]. The image has 46 rows and 81 columns (the total number of pixels is 46
×
81
=
3726). Since the logo only has 5 distinct rows, it has rank 5. We sample it using Gaussian i.i.d., Bernoulli i.i.d., Sparse Bernoulli i.i.d., and Random Projection measurements matrix with the number of linear constraints ranging between 700 and 2400. We declared MIT logo X to be recovered if X
-
X
opt
F
/
X
F
≤
1
0
-
3; that is, SNR ≥60 dB (S
N
R
=
-
20
log
⁡
(
X
-
X
opt
F
/
X
F
)). Figures 1 and 2 show the results of our numerical experiments.Figure 1
Performance curves under four different measurement matrices: SNR versus the number of linear constraints.Figure 2
Performance curves of nuclear norm heuristic method under four different measurement matrices: SNR versus the number of linear constraints.Figure1 plots the curves of SNR under four different measurement matrices mentioned above with different numbers of linear constraints. The numerical values on the x-axis represent the number of linear constraints p, that is, number of measurements, while those on y-axis represent the Signal to Noise Ration between the true original X and X
opt recovered by our method. Noticeably, Gauss and Sparse Binary can successfully reconstruct MIT logo images from about 1000 constraints, while Random Projection requires about 1200 constraints. Binary requires the minimal constraints about 950. The Gauss and Sparse Binary measurements reach their maximal SNR of 220 dB at around 1500 constraints; the Binary measurements can reach its maximal SNR of 180 dB at around 1400 constraints. Random Projection has more large phase transition point but has the best recovery error. We observe a sharp transition to perfect recovery near 1000 measurements lower than 2
k
(
m
+
n
-
k
). Compared with Figure 2, which shows results of nuclear norm heuristic method in [39] and the method that is a kind of convex relaxations, there is a sharp transition to perfect recovery near 1200 measurements which is approximately equal to 2
k
(
m
+
n
-
k
). The comparison shows our method has better performance that it can recover an image with less constraints; namely, it needs less samples. In addition, we observe that SNR of nuclear norm heuristic method achieves 120 dB when the number of measurements approximately equals 1500, while our method just needs 1200 measurements, which decrease by 20% compared with 1500. If the number of measurements is 1500, accuracy of our method would be higher. Among precisions of four different measurement matrixes, the smallest one is precision of Random Projection. However, the precision value is also higher than 120 dB. The performance illustrates that our method is more accurate. Observing the curve of between 1200 and 1500, we can find that our curve rises smoothly compared with the curve in the same interval in Figure 2; this shows that our method is stable for reconstructing image.Figure3 shows recovered images under four different measurement matrices with 800, 900, and 1000 linear constraints. And Figure 4 shows recovered images by our method and that by SeDuMi. We chose to use this interior-point method because it yields the highest accuracy. Distinctly, the recovered images by SeDuMi have some shadow, while ours is clear. Moreover, as far as time is concerned, our method is more fast and the computational time is not extremely sensitive to p compared to SeDuMi. The reason is that closed-form solutions can be obtained when we solve this problem. Comparisons are shown in Table 1.Table 1
DC versus SeDuMi in time (second).
p
800
900
1000
1100
1200
1500
DC
46.82
50.49
50.78
51.38
51.96
54.09
SeDuMi
93.26
120.16
138.24
159.44
161.76
245.38Figure 3
Recovered images under four measurement matrices with different numbers of linear constraints. From left to right: the original, results under Gauss, Binary, Sparse Binary, and Random Projection measurement. From top to bottom:p
=
800, p
=
900, and p
=
1000.Figure 4
Recovered images under different numbers of linear constraints. From left to right: the original, results under our method, and results by SeDuMi. From top to bottom:p
=
800, p
=
900, and p
=
1000.


### 3.2. Matrix Completion Problems

In this subsection, we apply our method to matrix completion problems. It can be reformulated as(48)
m
i
n
1
2
P
Ω
X
-
P
Ω
M
F
2
s
.
t
.
rank
⁡
X
≤
k
X
∈
R
m
×
n
,where X and M are both m
×
n matrices and Ω is a subset of index pairs (
i
,
j
). Up to now, numerous methods were proposed to solve the matrix completion problem (see, e.g., [9, 10, 29, 31]). It is easy to see that problem (48) is a special case of general rank minimization problem (3) with f
(
X
)
=
(
1
/
2
)
P
Ω
X
-
P
Ω
M
F
2.In our experiments, we first created random matricesM
L
∈
R
m
×
k and M
R
∈
R
m
×
k with i.i.d. Gaussian entries and then set M
=
M
L
M
R
T. We sampled a subset Ω of p entries uniformly at random. We compare against singular value thresholding (SVT) method [9], FPC, FPCA [10], and the OptSpace (OPT) method [29] using (49)
P
Ω
X
-
P
Ω
M
F
P
Ω
M
F
≤
1
0
-
3as a stop criterion. The maximum number of iterations is no more than 1500 in all methods; other parameters are default. We report result averaged over 10 runs.We useS
R
=
p
/
(
m
n
) to denote sampling ratio and F
R
=
k
(
m
+
n
-
k
)
/
p to represent the dimension of the set of rank k matrices divided by the number of measurements. NS and AT denote number of matrices that are recovered successfully and average times (seconds) that are successfully solved, respectively. RE represents the average relative error of the successful recovered matrices. We use the relative error (50)
M
-
X
∗
F
∖
M
Fto estimate the closeness of X
∗ to M, where X
∗ is the optimal solution to (48) produced by our method. We declared M to be recovered if the relative error was less than 1
0
-
3.Table2 shows that our method has the close performance with OptSpace in column of NS. These performances are better than those of SVT, FPC, and FPCA. However, SVT, FPC, and FPCA cannot obtain the optimal matrix within 1500 iterations which shows that our method has advantage in finding the neighborhood of optimal solution. The results do not mean that these algorithms are invalid, but the maximum number of iterations is limited to 1500. When the rank of matrix is large (i.e., k
≥
5), the average time of our algorithm is less than that of OptSpace. In some cases, the relative errors of DC are minimal. When the dimension of matrix becomes larger, DC method has a better performance of NS and similar relative errors with that of OptSpace. In terms of time, our time is longer than that of OptSpace, but shorter than that of FPC and SVT. Tables 3 and 4 show the above results. In conclusion, our method has better performance in “hard” problems; that is, F
R
>
0.34. The problems involved matrices that are not of very low rank, according to [10].Table 2
Numerical results for problems withm
=
n
=
40 and S
R
=
0.5.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0988
10
0.02
1.23
E
-
04
6
0.5550
9
0.12
4.64
E
-
04
SVT
10
0.11
1.59
E
-
04
0
—
—
FPC
9
0.13
6.52
E
-
04
0
—
—
FPCA
10
0.03
3.98
E
-
05
0
—
—
OptSpace
10
0.01
1.31
E
-
04
10
0.21
4.05
E
-
04
DC
2
0.1950
10
0.03
1.50
E
-
04
7
0.6388
9
0.17
6.11
E
-
04
SVT
9
0.18
1.83
E
-
04
0
—
—
FPC
8
0.14
6.52
E
-
04
0
—
—
FPCA
10
0.02
5.83
E
-
05
0
—
—
OptSpace
10
0.01
1.51
E
-
04
9
0.37
5.37
E
-
04
DC
3
0.2888
10
0.05
2.68
E
-
04
8
0.7200
7
0.28
7.90
E
-
04
SVT
2
0.36
2.47
E
-
04
0
—
—
FPC
5
0.25
8.58
E
-
04
0
—
—
FPCA
10
0.06
1.38
E
-
04
0
—
—
OptSpace
10
0.04
1.75
E
-
04
9
0.82
6.68
E
-
04
DC
4
0.3800
10
0.06
2.84
E
-
04
9
0.7987
4
0.40
8.93
E
-
04
SVT
1
0.62
1.38
E
-
04
0
—
—
FPC
1
0.25
8.72
E
-
04
0
—
—
FPCA
0
—
—
0
—
—
OptSpace
10
0.05
2.68
E
-
04
4
1.48
9.28
E
-
04
DC
5
0.4688
10
0.08
3.82
E
-
04
SVT
0
—
—
FPC
0
—
—
FPCA
0
—
—
OptSpace
10
0.09
3.16
E
-
04Table 3
Numerical results for problems withm
=
n
=
100 and SR = 0.2.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0995
10
0.30
1.91
E
-
04
5
0.4875
10
1.50
5.05
E
-
04
SVT
7
1.82
2.04
E
-
04
0
—
—
FPC
10
3.14
8.06
E
-
04
0
—
—
FPCA
10
0.12
1.73
E
-
04
0
—
—
OptSpace
10
0.03
1.71
E
-
04
10
0.58
5.57
E
-
04
DC
2
0.1980
10
0.77
3.34
E
-
04
6
0.5820
6
2.82
7.45
E
-
04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
8
0.21
2.97
E
-
04
0
—
—
OptSpace
10
0.05
2.27
E
-
04
9
0.97
6.57
E
-
04
DC
3
0.2955
10
0.64
3.03
E
-
04
7
0.6755
5
3.24
7.66
E
-
04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
6
0.09
2.20
E
-
04
0
—
—
OptSpace
10
0.12
3.76
E
-
04
5
1.89
7.95
E
-
04
DC
4
0.3920
9
1.16
4.40
E
-
04
SVT
0
—
—
FPC
0
—
—
FPCA
0
—
—
OptSpace
10
0.22
3.85
E
-
04Table 4
Numerical results for problems withm
=
n
=
100 and SR = 0.3.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0663
10
0.19
1.49E − 04
7
0.4530
10
0.61
3.64E − 04
SVT
10
0.68
1.66E − 04
0
—
—
FPC
10
1.33
4.47E − 04
0
—
—
FPCA
10
0.17
6.12E − 05
0
—
—
OptSpace
10
0.05
1.29E − 04
10
0.47
3.39E − 04
DC
2
0.1320
10
0.20
1.69E − 04
8
0.5120
10
0.87
4.70E − 04
SVT
8
1.39
1.66E − 04
0
—
—
FPC
10
2.18
4.84E − 04
0
—
—
FPCA
10
0.16
7.68E − 05
0
—
—
OptSpace
10
0.04
1.48E − 04
10
0.63
3.67E − 04
DC
3
0.1970
10
0.25
1.92E − 04
9
0.5730
9
1.03
4.85E − 04
SVT
4
1.77
1.23E − 04
0
—
—
FPC
9
4.32
5.82E − 04
0
—
—
FPCA
9
0.19
1.21E − 04
0
—
—
OptSpace
10
0.07
1.98E − 04
10
1.40
4.87E − 04
DC
4
0.2613
10
0.32
2.30E − 04
10
0.6333
8
1.36
6.21E − 04
SVT
4
9.82
4.67E − 04
0
—
—
FPC
7
5.02
6.59E − 04
0
—
—
FPCA
10
0.18
6.63E − 04
0
—
—
OptSpace
10
0.09
2.08E − 04
10
1.80
5.97E − 04
DC
5
0.3250
10
0.37
2.72E − 04
11
0.6930
9
2.19
7.29E − 04
SVT
0
—
—
0
—
—
FPC
1
8.18
7.34E − 04
0
—
—
FPCA
10
0.10
2.20E − 04
9
—
—
OptSpace
10
0.19
2.67E − 04
5
4.11
7.05E − 04
DC
6
0.3880
10
0.47
3.00E − 04
12
0.7520
6
2.83
8.79E − 04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
0
0.09
2.20E − 04
0
—
—
OptSpace
10
0.25
3.25E − 04
6
6.44
8.47E − 04Table5 gives the numerical results for problems with large dimensions. In these numerical experiments, the maximum number of iterations is chosen by default of each compared method. We observe that our method can recover problems successfully as other four methods, while our method is slower with the increase of dimension. However, our method is faster than FPC. Though SVT, FPCA, and OptSpace outperform our method in terms of speed, they are specifically developed method for matrix completion problems while our method can be applied to more class of problems. To prove this, we will apply our method to max-cut problems.Table 5
Numerical results for problems with different dimension.
m
(
n
)
k
p
SR
FR
Alg
NS
AT
RA
100
10
5666
0.57
0.34
DC
10
0.14
1.91E − 04
SVT
10
2.81
1.85E − 04
FPC
10
5.23
1.83E − 04
FPCA
10
0.28
4.94E − 05
OptSpace
10
0.9
1.41E − 04
200
10
15665
0.39
0.25
DC
10
0.80
1.65E − 04
SVT
10
2.93
2.53E − 04
FPC
10
4.18
2.73E − 04
FPCA
10
0.08
6.88E − 05
OptSpace
10
0.41
1.72E − 04
300
10
27160
0.30
0.22
DC
10
2.72
1.67E − 04
SVT
10
2.57
1.79E − 04
FPC
10
9.8
1.63E − 04
FPCA
10
0.39
6.71E − 05
OptSpace
10
1.85
1.50E − 04
400
10
41269
0.26
0.19
DC
10
5.73
1.56E − 04
SVT
10
2.78
1.66E − 04
FPC
10
13.28
1.30E − 04
FPCA
10
0.73
6.65E − 05
OptSpace
10
3.47
1.46E − 04
500
10
49471
0.20
0.2
DC
10
13.63
1.73E − 04
SVT
10
2.94
1.74E − 04
FPC
10
25.75
1.44E − 04
FPCA
10
1.12
1.07E − 04
OptSpace
10
5.05
1.60E − 04
1000
10
119406
0.12
0.17
DC
10
356.41
1.70E − 04
SVT
10
4.56
1.66E − 04
FPC
10
110.53
1.29E − 04
FPCA
10
5.16
1.64E − 04
OptSpace
10
14.61
1.50E − 04


### 3.3. Max-Cut Problems

In this subsection, we apply our method to max-cut problems. It can be formulated as follows:(51)
m
i
n
T
r
L
X
s
.
t
.
diag
⁡
X
=
e
rank
⁡
X
=
1
X
⪰
0
,where X
=
x
x
T, x
=
(
x
1
,
…
,
x
n
)
T, L
=
(
1
/
4
)
(
d
i
a
g
(
W
e
)
-
W
), e
=
(
1,1
,
…
,
1
)
T, and W is the weight matrix. Equation (51) is one general case of problem (3) with Γ
∩
Ω
=
{
X
∣
diag
⁡
(
X
)
=
e
X
∈
S
+
}. Hence, our method can be suitably applied to (51).We now address the initialization and termination criteria for our method when applied to (51). We set X
0
=
e
e
T; in addition, we choose the initial penalty parameter ρ to be 0.01 and set the parameter τ to be 1.1. And we use X
i
+
1
-
X
i
F
≤
ε as the termination criteria with ϵ
=
1
0
-
5.In our numerical experiments, we first test several small scale max-cut problems; then we perform our method on the standard “G” set problems written by Professor Ginaldi [40]. In all, we consider 10 problems ranging in size from 800 variables to 1000 variables. The results are shown in Tables 6 and 7, in which “m” and “k” represent corresponding dimension of graph and the rank of solution, respectively. The “obv” denotes objective values of optimal matrix solutions operated by our method.Table 6
Numerical results for famous max-cut problems.
Graph
m
KV
GW-cut
DC
Optimal solutions
UB
k
obv
k
C5
5
4
4.525
5
4
1
(−1; 1; −1; 1; −1)
K5
5
6
6.25
5
6
1
(−1; 1; −1; 1; −1)
AG
5
9.28
9.6040
5
9.28
1
(−1; 1; 1; −1; −1)
BG
12
88
90.3919
12
88
1
(−1; 1; 1; −1; −1; −1; −1; −1; −1; 1; 1; 1)Table 7
Numerical results for G problems.
Graph
m
SUB
GW-cut
DC
DC+
BK
obv
k
obv
k
G01
800
12078
10612
14
11281
1
11518
11624
G02
800
12084
10617
31
11261
1
11550
11620
G03
800
12077
10610
14
11275
1
11540
11622
G11
800
627
551
30
494
2
528
564
G14
800
3187
2800
13
2806
1
2980
3063
G15
800
3169
2785
70
2767
1
2969
3050
G16
800
3172
2782
14
2800
1
2974
3052
G43
1000
7027
6174
—
6321
1
6519
6660
G45
1000
7020
6168
—
6348
1
6564
6654
G52
1000
4009
3520
—
3383
1
3713
3849In Table6, “KV” denotes known optimal objective values and “UB” presents upper bounds for max-cut by the semidefinite relaxation. It is easy to observe that our method is effective for solving these problems and it has computed the global optimal solutions of all problems. Moreover, we obtain lower rank solution compared with semidefinite relaxation method.Table7 shows the numerical results for standard “G” set problems. “SUB” denotes upper bounds by spectral bundle method proposed by Helmberg and Rendl [41]. “GW-cut” denotes the approximative algorithm proposed by Goemans and Williamson [42]. “DC” denotes results generated by our method, D
C
+ is the results of DC with neighborhood approximation algorithm, and the right most column, titled BK, gives best results as we have known according to [43, 44]. We use “—” to indicate that we do not obtain the rank. From Table 3 we can note that our optimal values are close to the results of latest researches, better than values given by GW-cut, and they are all less than the corresponding upper bounds. For the problems, we all get the lower rank solution compared with GW-cut; moreover, most solutions satisfy the constraints r
a
n
k
=
1. The results prove that our algorithm is effective when used to solve max-cut problems.


## 3.1. Image Reconstruction Problems

In this subsection, we apply our method to solve one case of affine rank minimization problems. It can be formulated as(47)
m
i
n
1
2
A
X
-
b
F
2
s
.
t
.
rank
⁡
X
≤
k
X
∈
R
m
×
n
,where X
∈
R
m
×
n is the decision variable and the linear map A
:
R
m
×
n
→
R
p and vector b
∈
R
p are known. Clearly, (47) is a special case of problem (3). Hence, our method can be suitably applied to (47).Recht et al. [39] have demonstrated that when we sample linear maps from a class of probability distributions, then they would obey the Restricted Isometry Property. There are two ingredients for a random linear map to be nearly isometric. First, it must be isometric in expectation. Second, the probability of large distortions of length must be exponentially small. In our numerical experiments, we use four nearly isometric measurements. The first one is the ensemble with independent, identically distributed (i.i.d.) Gaussian entries. The second one has entries sampled from an i.i.d. symmetric Bernoulli distribution, which we call Bernoulli. The third one has zeros in two-thirds of the entries and others sampled from an i.i.d. symmetric Bernoulli distribution, which we call Sparse Bernoulli. The last one has an orthonormal basis, which we call Random Projection.Now we conduct numerical experiments to test the performance of our method for solving (47). To illustrate the scaling of low rank recovery for a particular matrix X, consider the MIT logo presented in [39]. The image has 46 rows and 81 columns (the total number of pixels is 46
×
81
=
3726). Since the logo only has 5 distinct rows, it has rank 5. We sample it using Gaussian i.i.d., Bernoulli i.i.d., Sparse Bernoulli i.i.d., and Random Projection measurements matrix with the number of linear constraints ranging between 700 and 2400. We declared MIT logo X to be recovered if X
-
X
opt
F
/
X
F
≤
1
0
-
3; that is, SNR ≥60 dB (S
N
R
=
-
20
log
⁡
(
X
-
X
opt
F
/
X
F
)). Figures 1 and 2 show the results of our numerical experiments.Figure 1
Performance curves under four different measurement matrices: SNR versus the number of linear constraints.Figure 2
Performance curves of nuclear norm heuristic method under four different measurement matrices: SNR versus the number of linear constraints.Figure1 plots the curves of SNR under four different measurement matrices mentioned above with different numbers of linear constraints. The numerical values on the x-axis represent the number of linear constraints p, that is, number of measurements, while those on y-axis represent the Signal to Noise Ration between the true original X and X
opt recovered by our method. Noticeably, Gauss and Sparse Binary can successfully reconstruct MIT logo images from about 1000 constraints, while Random Projection requires about 1200 constraints. Binary requires the minimal constraints about 950. The Gauss and Sparse Binary measurements reach their maximal SNR of 220 dB at around 1500 constraints; the Binary measurements can reach its maximal SNR of 180 dB at around 1400 constraints. Random Projection has more large phase transition point but has the best recovery error. We observe a sharp transition to perfect recovery near 1000 measurements lower than 2
k
(
m
+
n
-
k
). Compared with Figure 2, which shows results of nuclear norm heuristic method in [39] and the method that is a kind of convex relaxations, there is a sharp transition to perfect recovery near 1200 measurements which is approximately equal to 2
k
(
m
+
n
-
k
). The comparison shows our method has better performance that it can recover an image with less constraints; namely, it needs less samples. In addition, we observe that SNR of nuclear norm heuristic method achieves 120 dB when the number of measurements approximately equals 1500, while our method just needs 1200 measurements, which decrease by 20% compared with 1500. If the number of measurements is 1500, accuracy of our method would be higher. Among precisions of four different measurement matrixes, the smallest one is precision of Random Projection. However, the precision value is also higher than 120 dB. The performance illustrates that our method is more accurate. Observing the curve of between 1200 and 1500, we can find that our curve rises smoothly compared with the curve in the same interval in Figure 2; this shows that our method is stable for reconstructing image.Figure3 shows recovered images under four different measurement matrices with 800, 900, and 1000 linear constraints. And Figure 4 shows recovered images by our method and that by SeDuMi. We chose to use this interior-point method because it yields the highest accuracy. Distinctly, the recovered images by SeDuMi have some shadow, while ours is clear. Moreover, as far as time is concerned, our method is more fast and the computational time is not extremely sensitive to p compared to SeDuMi. The reason is that closed-form solutions can be obtained when we solve this problem. Comparisons are shown in Table 1.Table 1
DC versus SeDuMi in time (second).
p
800
900
1000
1100
1200
1500
DC
46.82
50.49
50.78
51.38
51.96
54.09
SeDuMi
93.26
120.16
138.24
159.44
161.76
245.38Figure 3
Recovered images under four measurement matrices with different numbers of linear constraints. From left to right: the original, results under Gauss, Binary, Sparse Binary, and Random Projection measurement. From top to bottom:p
=
800, p
=
900, and p
=
1000.Figure 4
Recovered images under different numbers of linear constraints. From left to right: the original, results under our method, and results by SeDuMi. From top to bottom:p
=
800, p
=
900, and p
=
1000.


## 3.2. Matrix Completion Problems

In this subsection, we apply our method to matrix completion problems. It can be reformulated as(48)
m
i
n
1
2
P
Ω
X
-
P
Ω
M
F
2
s
.
t
.
rank
⁡
X
≤
k
X
∈
R
m
×
n
,where X and M are both m
×
n matrices and Ω is a subset of index pairs (
i
,
j
). Up to now, numerous methods were proposed to solve the matrix completion problem (see, e.g., [9, 10, 29, 31]). It is easy to see that problem (48) is a special case of general rank minimization problem (3) with f
(
X
)
=
(
1
/
2
)
P
Ω
X
-
P
Ω
M
F
2.In our experiments, we first created random matricesM
L
∈
R
m
×
k and M
R
∈
R
m
×
k with i.i.d. Gaussian entries and then set M
=
M
L
M
R
T. We sampled a subset Ω of p entries uniformly at random. We compare against singular value thresholding (SVT) method [9], FPC, FPCA [10], and the OptSpace (OPT) method [29] using (49)
P
Ω
X
-
P
Ω
M
F
P
Ω
M
F
≤
1
0
-
3as a stop criterion. The maximum number of iterations is no more than 1500 in all methods; other parameters are default. We report result averaged over 10 runs.We useS
R
=
p
/
(
m
n
) to denote sampling ratio and F
R
=
k
(
m
+
n
-
k
)
/
p to represent the dimension of the set of rank k matrices divided by the number of measurements. NS and AT denote number of matrices that are recovered successfully and average times (seconds) that are successfully solved, respectively. RE represents the average relative error of the successful recovered matrices. We use the relative error (50)
M
-
X
∗
F
∖
M
Fto estimate the closeness of X
∗ to M, where X
∗ is the optimal solution to (48) produced by our method. We declared M to be recovered if the relative error was less than 1
0
-
3.Table2 shows that our method has the close performance with OptSpace in column of NS. These performances are better than those of SVT, FPC, and FPCA. However, SVT, FPC, and FPCA cannot obtain the optimal matrix within 1500 iterations which shows that our method has advantage in finding the neighborhood of optimal solution. The results do not mean that these algorithms are invalid, but the maximum number of iterations is limited to 1500. When the rank of matrix is large (i.e., k
≥
5), the average time of our algorithm is less than that of OptSpace. In some cases, the relative errors of DC are minimal. When the dimension of matrix becomes larger, DC method has a better performance of NS and similar relative errors with that of OptSpace. In terms of time, our time is longer than that of OptSpace, but shorter than that of FPC and SVT. Tables 3 and 4 show the above results. In conclusion, our method has better performance in “hard” problems; that is, F
R
>
0.34. The problems involved matrices that are not of very low rank, according to [10].Table 2
Numerical results for problems withm
=
n
=
40 and S
R
=
0.5.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0988
10
0.02
1.23
E
-
04
6
0.5550
9
0.12
4.64
E
-
04
SVT
10
0.11
1.59
E
-
04
0
—
—
FPC
9
0.13
6.52
E
-
04
0
—
—
FPCA
10
0.03
3.98
E
-
05
0
—
—
OptSpace
10
0.01
1.31
E
-
04
10
0.21
4.05
E
-
04
DC
2
0.1950
10
0.03
1.50
E
-
04
7
0.6388
9
0.17
6.11
E
-
04
SVT
9
0.18
1.83
E
-
04
0
—
—
FPC
8
0.14
6.52
E
-
04
0
—
—
FPCA
10
0.02
5.83
E
-
05
0
—
—
OptSpace
10
0.01
1.51
E
-
04
9
0.37
5.37
E
-
04
DC
3
0.2888
10
0.05
2.68
E
-
04
8
0.7200
7
0.28
7.90
E
-
04
SVT
2
0.36
2.47
E
-
04
0
—
—
FPC
5
0.25
8.58
E
-
04
0
—
—
FPCA
10
0.06
1.38
E
-
04
0
—
—
OptSpace
10
0.04
1.75
E
-
04
9
0.82
6.68
E
-
04
DC
4
0.3800
10
0.06
2.84
E
-
04
9
0.7987
4
0.40
8.93
E
-
04
SVT
1
0.62
1.38
E
-
04
0
—
—
FPC
1
0.25
8.72
E
-
04
0
—
—
FPCA
0
—
—
0
—
—
OptSpace
10
0.05
2.68
E
-
04
4
1.48
9.28
E
-
04
DC
5
0.4688
10
0.08
3.82
E
-
04
SVT
0
—
—
FPC
0
—
—
FPCA
0
—
—
OptSpace
10
0.09
3.16
E
-
04Table 3
Numerical results for problems withm
=
n
=
100 and SR = 0.2.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0995
10
0.30
1.91
E
-
04
5
0.4875
10
1.50
5.05
E
-
04
SVT
7
1.82
2.04
E
-
04
0
—
—
FPC
10
3.14
8.06
E
-
04
0
—
—
FPCA
10
0.12
1.73
E
-
04
0
—
—
OptSpace
10
0.03
1.71
E
-
04
10
0.58
5.57
E
-
04
DC
2
0.1980
10
0.77
3.34
E
-
04
6
0.5820
6
2.82
7.45
E
-
04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
8
0.21
2.97
E
-
04
0
—
—
OptSpace
10
0.05
2.27
E
-
04
9
0.97
6.57
E
-
04
DC
3
0.2955
10
0.64
3.03
E
-
04
7
0.6755
5
3.24
7.66
E
-
04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
6
0.09
2.20
E
-
04
0
—
—
OptSpace
10
0.12
3.76
E
-
04
5
1.89
7.95
E
-
04
DC
4
0.3920
9
1.16
4.40
E
-
04
SVT
0
—
—
FPC
0
—
—
FPCA
0
—
—
OptSpace
10
0.22
3.85
E
-
04Table 4
Numerical results for problems withm
=
n
=
100 and SR = 0.3.
Alg
k
FR
NS
AT
RE
k
FR
NS
AT
RA
DC
1
0.0663
10
0.19
1.49E − 04
7
0.4530
10
0.61
3.64E − 04
SVT
10
0.68
1.66E − 04
0
—
—
FPC
10
1.33
4.47E − 04
0
—
—
FPCA
10
0.17
6.12E − 05
0
—
—
OptSpace
10
0.05
1.29E − 04
10
0.47
3.39E − 04
DC
2
0.1320
10
0.20
1.69E − 04
8
0.5120
10
0.87
4.70E − 04
SVT
8
1.39
1.66E − 04
0
—
—
FPC
10
2.18
4.84E − 04
0
—
—
FPCA
10
0.16
7.68E − 05
0
—
—
OptSpace
10
0.04
1.48E − 04
10
0.63
3.67E − 04
DC
3
0.1970
10
0.25
1.92E − 04
9
0.5730
9
1.03
4.85E − 04
SVT
4
1.77
1.23E − 04
0
—
—
FPC
9
4.32
5.82E − 04
0
—
—
FPCA
9
0.19
1.21E − 04
0
—
—
OptSpace
10
0.07
1.98E − 04
10
1.40
4.87E − 04
DC
4
0.2613
10
0.32
2.30E − 04
10
0.6333
8
1.36
6.21E − 04
SVT
4
9.82
4.67E − 04
0
—
—
FPC
7
5.02
6.59E − 04
0
—
—
FPCA
10
0.18
6.63E − 04
0
—
—
OptSpace
10
0.09
2.08E − 04
10
1.80
5.97E − 04
DC
5
0.3250
10
0.37
2.72E − 04
11
0.6930
9
2.19
7.29E − 04
SVT
0
—
—
0
—
—
FPC
1
8.18
7.34E − 04
0
—
—
FPCA
10
0.10
2.20E − 04
9
—
—
OptSpace
10
0.19
2.67E − 04
5
4.11
7.05E − 04
DC
6
0.3880
10
0.47
3.00E − 04
12
0.7520
6
2.83
8.79E − 04
SVT
0
—
—
0
—
—
FPC
0
—
—
0
—
—
FPCA
0
0.09
2.20E − 04
0
—
—
OptSpace
10
0.25
3.25E − 04
6
6.44
8.47E − 04Table5 gives the numerical results for problems with large dimensions. In these numerical experiments, the maximum number of iterations is chosen by default of each compared method. We observe that our method can recover problems successfully as other four methods, while our method is slower with the increase of dimension. However, our method is faster than FPC. Though SVT, FPCA, and OptSpace outperform our method in terms of speed, they are specifically developed method for matrix completion problems while our method can be applied to more class of problems. To prove this, we will apply our method to max-cut problems.Table 5
Numerical results for problems with different dimension.
m
(
n
)
k
p
SR
FR
Alg
NS
AT
RA
100
10
5666
0.57
0.34
DC
10
0.14
1.91E − 04
SVT
10
2.81
1.85E − 04
FPC
10
5.23
1.83E − 04
FPCA
10
0.28
4.94E − 05
OptSpace
10
0.9
1.41E − 04
200
10
15665
0.39
0.25
DC
10
0.80
1.65E − 04
SVT
10
2.93
2.53E − 04
FPC
10
4.18
2.73E − 04
FPCA
10
0.08
6.88E − 05
OptSpace
10
0.41
1.72E − 04
300
10
27160
0.30
0.22
DC
10
2.72
1.67E − 04
SVT
10
2.57
1.79E − 04
FPC
10
9.8
1.63E − 04
FPCA
10
0.39
6.71E − 05
OptSpace
10
1.85
1.50E − 04
400
10
41269
0.26
0.19
DC
10
5.73
1.56E − 04
SVT
10
2.78
1.66E − 04
FPC
10
13.28
1.30E − 04
FPCA
10
0.73
6.65E − 05
OptSpace
10
3.47
1.46E − 04
500
10
49471
0.20
0.2
DC
10
13.63
1.73E − 04
SVT
10
2.94
1.74E − 04
FPC
10
25.75
1.44E − 04
FPCA
10
1.12
1.07E − 04
OptSpace
10
5.05
1.60E − 04
1000
10
119406
0.12
0.17
DC
10
356.41
1.70E − 04
SVT
10
4.56
1.66E − 04
FPC
10
110.53
1.29E − 04
FPCA
10
5.16
1.64E − 04
OptSpace
10
14.61
1.50E − 04


## 3.3. Max-Cut Problems

In this subsection, we apply our method to max-cut problems. It can be formulated as follows:(51)
m
i
n
T
r
L
X
s
.
t
.
diag
⁡
X
=
e
rank
⁡
X
=
1
X
⪰
0
,where X
=
x
x
T, x
=
(
x
1
,
…
,
x
n
)
T, L
=
(
1
/
4
)
(
d
i
a
g
(
W
e
)
-
W
), e
=
(
1,1
,
…
,
1
)
T, and W is the weight matrix. Equation (51) is one general case of problem (3) with Γ
∩
Ω
=
{
X
∣
diag
⁡
(
X
)
=
e
X
∈
S
+
}. Hence, our method can be suitably applied to (51).We now address the initialization and termination criteria for our method when applied to (51). We set X
0
=
e
e
T; in addition, we choose the initial penalty parameter ρ to be 0.01 and set the parameter τ to be 1.1. And we use X
i
+
1
-
X
i
F
≤
ε as the termination criteria with ϵ
=
1
0
-
5.In our numerical experiments, we first test several small scale max-cut problems; then we perform our method on the standard “G” set problems written by Professor Ginaldi [40]. In all, we consider 10 problems ranging in size from 800 variables to 1000 variables. The results are shown in Tables 6 and 7, in which “m” and “k” represent corresponding dimension of graph and the rank of solution, respectively. The “obv” denotes objective values of optimal matrix solutions operated by our method.Table 6
Numerical results for famous max-cut problems.
Graph
m
KV
GW-cut
DC
Optimal solutions
UB
k
obv
k
C5
5
4
4.525
5
4
1
(−1; 1; −1; 1; −1)
K5
5
6
6.25
5
6
1
(−1; 1; −1; 1; −1)
AG
5
9.28
9.6040
5
9.28
1
(−1; 1; 1; −1; −1)
BG
12
88
90.3919
12
88
1
(−1; 1; 1; −1; −1; −1; −1; −1; −1; 1; 1; 1)Table 7
Numerical results for G problems.
Graph
m
SUB
GW-cut
DC
DC+
BK
obv
k
obv
k
G01
800
12078
10612
14
11281
1
11518
11624
G02
800
12084
10617
31
11261
1
11550
11620
G03
800
12077
10610
14
11275
1
11540
11622
G11
800
627
551
30
494
2
528
564
G14
800
3187
2800
13
2806
1
2980
3063
G15
800
3169
2785
70
2767
1
2969
3050
G16
800
3172
2782
14
2800
1
2974
3052
G43
1000
7027
6174
—
6321
1
6519
6660
G45
1000
7020
6168
—
6348
1
6564
6654
G52
1000
4009
3520
—
3383
1
3713
3849In Table6, “KV” denotes known optimal objective values and “UB” presents upper bounds for max-cut by the semidefinite relaxation. It is easy to observe that our method is effective for solving these problems and it has computed the global optimal solutions of all problems. Moreover, we obtain lower rank solution compared with semidefinite relaxation method.Table7 shows the numerical results for standard “G” set problems. “SUB” denotes upper bounds by spectral bundle method proposed by Helmberg and Rendl [41]. “GW-cut” denotes the approximative algorithm proposed by Goemans and Williamson [42]. “DC” denotes results generated by our method, D
C
+ is the results of DC with neighborhood approximation algorithm, and the right most column, titled BK, gives best results as we have known according to [43, 44]. We use “—” to indicate that we do not obtain the rank. From Table 3 we can note that our optimal values are close to the results of latest researches, better than values given by GW-cut, and they are all less than the corresponding upper bounds. For the problems, we all get the lower rank solution compared with GW-cut; moreover, most solutions satisfy the constraints r
a
n
k
=
1. The results prove that our algorithm is effective when used to solve max-cut problems.


## 4. Concluding Remarks

In this paper we used closed-form of penalty function to reformulate convex constrained rank minimization problems as DC programming, which could be solved by stepwise linear approximative method. Under some suitable assumptions, we showed that accumulation point of the sequence generated by the stepwise linear approximative method was a stable point. The numerical results on affine rank minimization and max-cut problems demonstrate that our method is generally comparable or superior to some existing methods in terms of solution quality, and it can be applied to a class of problems. However, the speed of our method needs to be improved.



---

*Source: 7473041-2016-07-27.xml*