arXiv:1512.01748v1 [cs.NA] 6 Dec 2015

Restricted Low-Rank Approximation via ADMM
Ying Zhang
zy013@ie.cuhk.edu.hk
December 8, 2015
Abstract

The matrix low-rank approximation problem with additional convex constraints can find many applications and
has been extensively studied before. However, this problem is shown to be nonconvex and NP-hard; most of the
existing solutions are heuristic and application-dependent.
In this paper, we show that, other than tons of application
in current literature, this problem can be used to recover
a feasible solution for SDP relaxation. By some sophisticated tricks, it can be equivalently posed in an appropriate
form for the Alternating Direction Method of Multipliers
(ADMM) to solve. The two updates of ADMM include
the basic matrix low-rank approximation and projection
onto a convex set. Different from the general non-convex
problems, the sub-problems in each step of ADMM can
be solved exactly and efficiently in spite of their nonconvexity. Moreover, the algorithm will converge exponentially under proper conditions. The simulation results
confirm its superiority over existing solutions. We believe
that the results in this paper provide a useful tool for this
important problem and will help to extend the application
of ADMM to the non-convex regime.

1

Research Proposal

Notations and Operators
Vectors and matrices are denoted by boldface lower and upper case letters respectively. The set of real and natural numbers are represented by R and N.
The set of vectors and matrices with proper sizes are denoted as Rn , Rm×n . The
vectors are by default column vectors and the ith entry of vector x is denoted
as xi . The entry of matrix A in the ith row and j th column is denoted as Ai,j .
The superscript (·)T stands for transpose.
k · kp denotes the p-norm of a vector (p ≥ 1), i.e.,
kxkp =

X
i

|xi |p

! p1

and k · kF denotes the F −norm of a matrix, i.e.,


kAkF = 

X
i,j

 21

A2i,j  .

To be consistent with Matlab operator, diag (A) returns a vector a with
ai = A(i,i if the input A is a matrix and diag (a) returns a matrix A with
ai if i == j,
Ai,j =
. rank(X) returns the rank of matrix X.
0, otherwise

2

1

Research Proposal

1

INTRODUCTION

Introduction

There is a common belief that the complexity of the systems we study is up to
a limited level and the useful information in our observation is usually sparse.
The property of sparsity can be applied to reduce the necessary sampling rate
of signal reconstruction in compressed sensing [6], to avoid over-fitting in regularized machine learning algorithms [24], to increase robustness against noise
in solving linear equations, to train the deep neural network in unsupervised
learning [16] and beyond.
A direct interpretation of sparsity is that there are many zeros entries in
the data. More often, sparsity is reflected not directly by the data we collect,
but by the underlying pattern to be discovered. For example, the sparsity of
the voice will be more prominent in frequency domain [2]; the sparse pattern
of some astronomical or biomedical imaging data will be reviewed by wavelet
transformation [25, 26]. In recent years, we are able to generate or collect a lot
of data samples with many features, and different data samples and features
are usually correlated with each other. If we form the data into a matrix, the
matrix is often of low rank, which is another representation of sparsity and can
be exploited in many applications [8]. To encourage sparsity, we can incorporate
the sparsity interpretation into the objective function as a regularization term,
or put a sparsity upper bound as a hard constraint.
The sparsity usually leads to a non-convex problem and the problem is intractable. For example, we can let the number of nonzero component (0-norm)
in the data to be upper bounded. However, the function of 0-norm is nonconvex. To tackle the non-convexity, the researchers proposed many solutions,
usually belonging to the following two categories. First, we can use its convex envelop to replace the non-convex function, for example, to use 1-norm to
replace 0-norm, and solve the convex problem instead (this method is usually
called convex relaxation). And we can study under what condition the convex
relaxation is exact or how much performance loss we will suffer due to the relaxation [6]. Second, we can directly design some heuristic algorithms to solve the
non-convex problem based on engineering intuition, like the Orthogonal Matching Pursuit (OMP) algorithm in compressed sensing [21]. This kind of method
usually works well in practice but little theoretical performance guarantee can
be obtained.

1.1

Related work

The rank of a matrix, as the dimensionality of the smallest subspace to which the
data belongs to, is essentially important and the research on it can be traced
back to the origin of matrix. In recent years, many problems with low-rank
matrix are of our interests. However, its non-convex nature makes the problems
generally intractable. In this case, some researchers proposed to use the nuclear
norm or log-det function of the matrix to replace the rank function and then
solve the relaxed convex problem [9, 23]. In [1], the matrix being of low rank is
interpreted as the matrix being factorized into two smaller-sized matrices, then
3

1

Research Proposal

INTRODUCTION

the matrix factorization techniques can be applied.
The alternating direction method of multiplier (ADMM) was invented in
the 1970s and is raising its popularity in the era of big data; many large scale
optimization problems that arise in practice can be formulated or equivalently
posed in a form appropriate for ADMM and distributed algorithms can be obtained [3]. Many successful applications including consensus and sharing have
been found. The theory related to the convergence rate of ADMM is an active
ongoing topic and many results for the convex problem are provided in [3] and
the references therein. In practice, we can meet many problems that are nonconvex and ADMM algorithm is proposed to tackle them in a heuristic sense,
such as optimal power flow problem [28], matrix factorization [4], etc. However,
in the non-convex regime, little theoretical guarantee can be obtained. To the
best of our knowledge, [12, 19] are among the first attempts to characterize the
convergence behavior of ADMM for the general non-convex problems and has
limited theoretical performance guarantee.
In this paper, we consider that low-rank approximation problem with additional convex constraints, which is non-convex and NP-hard. By borrowing
some tricks from [3] [28], we can reformulate the problem into the form that can
be solved by ADMM. In each iteration of ADMM, the first update is to solve a
convex problem and the second update is to solve a standard low-rank approximation, the optimal solution of which can be obtained by singular value decomposition. Different from the general non-convex problems, the sub-problems
in ADMM can be solved efficiently. This approach is motivated by [28], in
which the authors also use ADMM to solve a non-convex problem with similar
structures. The different thing is that, due to the special objective function we
consider, i.e., kX − X̂kF , we can show that the first update is a projection onto
a convex set and then establish the convergence of primal variable by assuming
the convergence of dual variable, which is an appealing result for the ADMM
application in non-convex regime.
The remaining part of this paper is organized as follows,
⊲ In Section 2, we introduce the low rank approximation problem with
additional convex constraints and show that it can capture the structured low
rank approximation and feasible solution recovery of SDP relaxation as two
special cases.
⊲ In Section 3, we leverage ADMM to design a generic algorithm and shows
that each step of this algorithm can be solved efficiently. We also provide some
theoretical results to characterize its convergence behavior.
⊲ In Section 4, we shows the performance of our algorithm by extensive
evaluations with synthetic and real-world data.
⊲ Section 5 is for conclusion and future work.

4

Research Proposal

2
2.1

2

RESTRICTED LOW RANK APPROXIMATION

Restricted Low Rank Approximation
Problem formulation

In this paper, we are particularly interested in the data-fitting problem and
restrict our attention to the low-rank solutions. More precisely, we are given a
matrix X̂ and we want to find a low-rank matrix to approximate X̂. The first
version of this problem is formulated as the Low Rank Approximation problem
(LRA) as follows,
LRA

min

kX̂ − XkF

X

s.t.

rank(X) ≤ K;

var.

(1a)

X,

where K is an integer to specify the upper bound of the matrix rank. The
problem is non-convex due to the rank constraint but an optimal solution can
be given by the well known Eckart-Young-Mirsky Theorem,
Theorem 1 (Eckart-Young-Mirsky Theorem [7]) If the matrix X̂ admits
the singular value decomposition X̂ = U ΣV H with Σ = diag([σ1 , σ2 , · · · , σn ])
and σ1 ≥ σ2 ≥ · · · ≥ σn ≥ 0, an optimal solution to problem LRA is given by
X∗ =

K
X

σk uk vkH .

k=1

Furthermore, the minimizer is unique if σK and σK+1 are not equal.
This method is called truncated singular value decomposition (SVD). We
can see that even though LRA is non-convex, it can be solved in polynomial
time because SVD can be computed in polynomial time.
In this paper, we put some additional constraints on LRA and call the
new problem as Restricted Low Rank Approximation problem (RLRA). The
problem is casted as follows,
RLRA

kX̂ − Xk2F

min
X

s.t.

rank(X) ≤ K;
g(X) ≤ 0;

var.

(2a)
(2b)

X,

where g(X) is a convex function, requiring that the approximation X is located
in a convex set.
The difference between LRA and RLRA is the constraint g(X) ≤ 0. The
convexity of g(X) will makes the readers feel that the difference is not significant,
since the convex things are usually treated as easy in most literature. However,
this is not true. The new constraint makes the truncated SVD not applicable
5

Research Proposal

2

RESTRICTED LOW RANK APPROXIMATION

to solve RLRA 1 . More importantly, RLRA is believed to be NP-hard [20].
There is no hope to always achieve the optimal solution in polynomial time
unless P=NP.
There are several proper ways to tackle this non-convex problem. Firstly,
we can replace the rank constraints by some constraints on the nuclear norm
of X [8] and solve the relaxed problem instead; secondly, we can use the kernel
representation, for example,to equivalently transform the rank constraint to the
fact that the matrix can be factorized into two smaller sized matrices and then
apply the matrix factorization technique [20]. In this paper, we propose to use
ADMM, which has found its success in many problems [3], to solve this problem.

2.2

Two specific instances

The generic problem RLRA can be used to solve many problems and we review
two of them in this section to provide more motivations.
2.2.1

Low-rank approximation with linear structures

The data fitting task can be formulated as a structured low-rank approximation problem (SLRA)if the system generating the data is of linear model and
bounded complexity [20]. Many applications can be found in system theory,
signal processing, computer algebra, etc [20].
We denote A a set of matrices with specific affine structures, for example,
Hankel matrix, Toplitiz matrix, etc. The problem is formally given as follows,
kX − X̂k2F

SLRA min
x

s.t.
var.

rank(X) ≤ K,

(3a)

X ∈ A,
X ∈ Cm×n .

(3b)

Different applications lead to different requirements of A [20]. Without
diving into to the details, we list some of them here,
• Hankel matrix: approximate realization, model reduction, output error
identification;
• Sylvester matrix: Pole place by low-order controller, approximate common
divisor;
• Hankel&Toeplitz matrix: Harmonic retrieval;
• Non-negative matrix: image mining, Markov chains.
Some heuristic or local-optimization based algorithms for different problems
are summarized in [20].
1 It

is easy to imagine that the solution by truncated SVD may not respect g(X) ≤ 0.

6

Research Proposal

2.2.2

2

RESTRICTED LOW RANK APPROXIMATION

Feasible solution recovery of SDP relaxation

Many communication problems, like multicast downlink transmit beamforming
problem, can be formulated as a quadratically constrained quadratic program
QCQP, which is non-convex and generally NP-hard,
QCQP

min
x

s.t.

xH Cx
xH Fi x ≥ gi , i = 1, ..., k
xH Hi x = li , i = 1, ..., m
var.x ∈ Rn ,

and it can be shown to be equivalent to problem SDP-QCQP.
QCQP-SDP

min
x

s.t.

var.

trace(C · X)
trace(Fi · X) ≥ gi , i = 1, ..., k

(5a)

trace(Hi · X) = li , i = 1, ..., m
rank(X) = 1,

(5b)
(5c)

X  0, X ∈ Cn×n .

By dropping the rank-1 constraint (5c) we can have a standard SDP problem
(denoted as QCQP-SDPR) and it can be solved by standard solver like CVX
[10]. This technique is called SDP relaxation and more details can be found
in [18].
If the optimal solution of the relaxed problem, denoted as X̂, happens to
respect the rank-1 constraint, then the optimal solution of the problem QCQP,
denoted as x∗ , can be obtained by the fact that X̂ = x∗ x∗T . In this case, the
SDP relaxation is called exact relaxation, and the original QCQP problem can
be solved efficiently and exactly even though it is non-convex. For the optimal
power flow problem, the exact relaxation always happens if some conditions
hold [17] [15].
However, the rank of X̂ is more often larger than 1 2 and the relaxation is
not exact. In this case, another step is needed to recover a good solution if we
do not want to waste pervious effort. The most direct way is to find the rank-1
matrix that is closest to X̂, i.e., solving LRA with K = 1 by truncated SVD.
This approach is suggested in [30] for the state estimation problem of power
system. However, we want to point out that this approach is not guaranteed to
produce a feasible solution of the original problem, because the rank-1 matrix
by truncated SVD may not satisfy the other constraints of SDP-QCQP like
(5a) (5b). We provide a successful case in Fig 1 and an unsuccessful case in
Fig 2 for a more clear illustration; X̃ is the optimal solution of QCQP-SDPR
and X1 is the rank-1 approximation of X̃.
2 Otherwise

SDP relaxation will always solve a non-convex problem exactly, which is not

true

7

3

Research Proposal

ALGORITHM DESIGN

Figure 1: An example that trun- Figure 2: An example that truncated SVD works for FSR-SDPR cated SVD fails to work for FSRSDPR
Some heuristic and problem-dependent algorithms are proposed in [18]. For
example, if it is required that kxk2 ≤ 1 for QCQP, we can first obtain a solution
x̂ from the rank-1 matrix by truncated SVD and then normalize the solution
x̂
if kx̂k2 > 1.
by kx̂k
2
Here we propose a generic method to solve the problem FSR-SDPR (Feasible Solution Recovery of SDP Relaxation) by finding a rank-1 matrix from
which we can recover a feasible solution to QCQP.
FSR-SDPR

min

kX − X̂kF

x

s.t.

g(X) ≤ 0,
rank(X) = 1,
X  0, X ∈ Cn×n ,

var.

where g(X) ≤ 0 represents the original convex constraints like (5a) and (5b).
The rationale behind this approach is that we want to find the point closest
to X̂ in the feasible region. If the objective function, trace(C · X), is Lipschitz
continuous, the performance of the recovered solution is guaranteed to be close
to the optimal value of QCQP-SDPR, hence close to the optimal solution of
QCQP. It is not difficult to see that solving this problem can be viewed as a
special case of RLRA and the solution is far from being trivial.

3

Algorithm Design

Next we will present how to solve the problem RLRA via the ADMM algorithm.

3.1

ADMM algorithm

In this section, we review the basic version of ADMM [3] to bring all the readers
to the same page 3 . We present the standard problem that ADMM can solve
3 The

readers are recommended to read [3] for more details.

8

3

Research Proposal

ALGORITHM DESIGN

in SP, in which the objective function is separable and two variables x, y are
coupled with each other by a linear constraint.
SP

min

f (x) + g(y)

s.t.
var.

Ax + By = c,
x, y.

x,y

(7a)

The augmented Lagrange multiplier function of the above problem is given
by

ρ
Lρ (x, y, λ) = f (x) + g(y) + λ(Ax + By − c) + kAx + By − ck22 ,
2
the scaled form of which is
ρ
Lρ (x, y, u) = f (x) + g(y) + kAx + By − c + uk22 ,
2

(8)

and we will use the scaled form in the sequel unless specified.
In each iteration, ADMM algorithm consists of the following three steps
• x update:
xk+1 = argminx Lρ (x, y k , uk ).
• y update:
y k+1 = argminy Lρ (xk+1 , y, uk ).
• u update:
uk+1 = uk + Axk+1 + By k+1 − c.
The ADMM algorithm is guaranteed to converge under some conditions, as
shown in Theorem 2.
Theorem 2 ( [3]) If f (x) and g(y) are closed, proper , convex and the unaugmented Lagrange function L0 has a saddle point, ADMM algorithm is guaranteed
to converge in the following sense,
• The residual Ax + By − c will converge to 0, i.e., the solution xk , y k will
approach feasibility.
• The objective value f (xk ) + g(y k ) will converge to the optimality.
• The dual variable u will converge to the dual optimal point.
We highlight that under the conditions provided in Theorem 2, the primal
variables x, y are not guaranteed to converge. But since the objective value is
guaranteed to converge, all the solutions after enough iterations will produce
the same results and are equivalently good.

9

3

Research Proposal

3.1.1

ALGORITHM DESIGN

Some remarks

We remark some properties of ADMM here. First, the ADMM algorithm can
produce a reasonably accurate solution very fast, but takes more time to generate a high accurate solution. For many machine learning or statistical tasks,
the overall performance depends on both problem formulation (including feature engineering and model selection) and parameter estimation (solving the
optimization problem). Usually the performance bottleneck comes from the
first one and a more accurate solution will not lead to significant performance
improvement. In this way, a reasonably accurate solution is totally acceptable
and this property of ADMM is appealing. Second, the optimization variable
of Problem SP is (x, y). In conventional algorithm such as gradient descent,
x, y are updated simultaneously, but in ADMM, we update x first and use the
newly updated x to update y. The intuition is that we want the newly updated information to take effect as soon as possible, which is similar to the
logic of Gauss-Sidel algorithm. Last but not least, the ADMM algorithm can
be implemented in a distributed manner with a parameter server [3].

3.2

Problem reformulation and applying ADMM

The original problem formulation in RLRA does not have the form of SP to
be readily solved by ADMM. We need to reformulate the problem. Firstly we
define two indicator functions as follows,
(
0, if rank(X) ≤ K,
I(X) =
(9)
+∞, otherwise,
(
0, if g(X) ≤ 0,
J (X) =
(10)
+∞, otherwise.
And we reformulate Problem RLRA equivalent to Problem RLRA-ADMM
as,
RLRA-ADMM

min
X,Y

kX − X̂kF + J (X) + I(Y )
|
{z
}
| {z }
X involved

s.t.
var.

X −Y =0

Y involved

(11a)

X, Y.

The equivalence of the two problems is formally established in Lemma 1.
Lemma 1 If Problem RLRA is feasible, then the optimal solutions and optimal
objective values of RLRA and RLRA-ADMM are the same.
And more importantly, RLRA-ADMM is readily solved by AMDD algorithm.
We denote f (X) = kX − X̂kF +I(X), g(Y ) = J(Y ) and the augmented function
with dual variable U as
ρ
Lρ (X, Y, U ) = f (X) + g(Y ) + kX − Y + U k2F .
2
10

3

Research Proposal

ALGORITHM DESIGN

Following the ADMM procedure, we will have the following three updates in
each iteration.
• X update:
X k+1 = argminX Lρ (X, Y k , U k ).
• Y update:
Y k+1 = argminY Lρ (X k+1 , Y, U k ).
• U update:
U k+1 = U k + X k+1 − Y k+1 .
In the next part, we carefully study the details of the algorithm and show
that each update can be carried out efficiently even though the problems can
be non-convex.
3.2.1

The subproblems

The U update is simple and direct. We revisit the other two in this part.
3.2.1.1

X update

The optimization problem in X update is equivalent to
X-MIN min
X

s.t.

ρ
kX − X̂k2F + kX − Y k + U k k2F
2
g(X) ≤ 0.

The problem X-MIN is convex and can be solved efficiently. We wit can be
viewed as a projection of Xk+1 + Uk onto the convex set S = {Y|g(Y) ≤ 0}.
It can be solved efficiently or even have closed-form solutions.
Lemma 2 The optimal solution of X-MIN can be obtained by solving
min
X

kX −

s.t.

1
1+

ρ
2




ρ
X̂ + (Y k − U k ) k2F
2

g(X) ≤ 0.


which is a projection of 1+1 ρ X̂ + ρ2 (Y k − U k ) onto the convex set {X|g(X) ≤
2

0}.

Proof We prove the equivalence between the two problems by showing the
linear relationship between the two objective functions. We will use the equation
that kAk2F = trace(AAH ) in the proof.

11

3

Research Proposal

ALGORITHM DESIGN

ρ
kX − X̂k2F + kX − Y k + U k k2F
2
ρ
=trace((X − X̂)(X − X̂)H ) + trace((X − Y k + U k )(X − Y k + U k )H )


 2 ρ

ρ
H
=trace 1 +
XX − 2 X̂ + (Y k − U k ) X H + C1
2
2
!



 
H


1
ρ
1
ρ k
ρ
trace
X−
= 1+
X−
X̂ + (Y − U k )
X̂ + (Y k − U k )
2
1 + ρ2
2
1 + ρ2
2
+ C2
ρ
1
= 1+
kX −
2
1+


ρ
2




ρ
X̂ + (Y k − U k ) k2F + C2 .
2

The proof is completed.
3.2.1.2

Y update

The optimization problem in Y update is equivalent to
Y-MIN

min
Y

s.t.

ρ k+1
kX
− Y + Uk k2F
2
rank(Y) ≤ 0.

The above optimization problem is non-convex because the rank constraint
is non-convex. Thus the problem is challenging on the first sight. The following
lemma shows that its optimal solution can be obtained by truncated SVD. We
summarize the algorithm into Algorithm 1 to end this part.
Algorithm 1 ADMM-RLRA: ADMM algorithm for RLRA
Require: pm ,pg ,pe (t),ek (t)
Ensure: uk (t),v k (t)
1: initialization
2: while not terminate do
3:
Obtain Xk+1 by solving X-MIN
4:
Obtain Yk+1 by solving Y-MIN
5:
Uk+1 = Uk + Xk+1 − Yk+1
6:
k =k+1
7: end while

3.3

On the convergence of the algorithm

The convergence of ADMM for the non-convex problem is an open problem [3].
Some positive results are obtained with some assumptions [12, 19]. We provide
some preliminary results in this section, regarding to the convergence and feasibility of Algorithm 1. An assumption that the dual variable Uk converges is
12

3

Research Proposal

ALGORITHM DESIGN

used in the theoretical analysis, which is also the assumption for the convergence analysis of ADMM for polynomial optimization [14], non-negative matrix
factorization [29] and non-negative matrix factorization [27].
Lemma 3 If the dual variable Uk converges, the solution Xk in the X update
of Algorithm 1 will approach feasibility of RLRA.
Proof Since Uk converges, Xk − Yk converges to 0. Since Yk satisfies the
constraint g(X) ≤ 0 for all k, Xk will satisfy the same condition when k is large
enough. Meanwhile, X k satisfies the rank constraint for all k, then Xk is a
feasible solution for RLRA.
The proof is completed.
When we characterize the convergence of ADMM in Theorem 2, the primal
variables are not guaranteed to converge, but can oscillate in the optimal region.
For the non-convex problem, because we do not know whether the objective
function will converge or not, the convergence of the primal variable is more
important, which we will discuss next.
Since Uk converges Ū, Xk − Yk converges to 0, then the X update can be
denoted as

1 
ρ k
k
2
Xk+1 = argmin kX −
ρ X̂ + (X − U ) kF .
1
+
2
rank(X)≤K
2


Let D = X̂ − ρ2 Ū , it can be further simplified as
Xk+1 = argmin kX −
rank(X)≤K

1
1+

ρ
2

ρ
2


Xk + D k2F


= argminkX − αXk + (1 − α)D k2F ,

α=

g(X)≤0

ρ
2+ρ

= C(Xk ).

Based on this understanding, we have the following theorem to characterize
its convergence.
Theorem 3 If the dual variable converges, the primal variable will converge,
and when k is large enough, we can have
kXk+2 − Xk+1 kF ≤

ρ
kXk+1 − Xk kF .
ρ+2

Proof X k+2 is the projection of αX k+1 + (1 − α)D and X k+1 is the projection
of αX k + (1 − α)D, with the fact that the feasible region of the projections is
convex, and projection onto a convex set is non-expansive [11] we can have
kXk+2 − Xk+1 kF ≤ kαX k+1 + (1 − α)D − αX k − (1 − α)DkF
= αkX k+1 − X k kF
The proof is completed.
13

3

Research Proposal

ALGORITHM DESIGN

We provide some remarks regarding this theorem. Firstly, the convergence of X k
will guarantee that we can have a local optimal or stationary point. Secondly,
the inequality indicates that the primal variable will converge exponentially and
a smaller ρ will lead to a faster convergence.

3.4

Another ADMM

3.4.1

Algorithm design

In this part, we provide another ADMM algorithm, which considers different
constraints compared with the algorithm we previously proposed, i.e., in X
update we consider the rank constraint while in Y update we consider the
convex constraint.
We denote f˜(X) = kX − X̂kF + I(X), g̃(Y ) = J(Y ) and the augmented
function with dual variable U as
ρ
L̃ρ (X, Y, U ) = f˜(X) + g̃(Y ) + kX − Y + U k2F .
2
Following the ADMM procedure, we will have the following three updates in
each iteration.
• X update:
X k+1 = argminX L̃ρ (X, Y k , U k ).
• Y update:
Y k+1 = argminY L̃ρ (X k+1 , Y, U k ).
• U update:
U k+1 = U k + X k+1 − Y k+1 .
3.4.1.1

The Subproblems

By similar arguments, the optimization problem in X update is equivalent to

1 
ρ k
k
min
kX −
(Y
−
U
)
k2F
X̂
+
X
1 + ρ2
2
s.t.
rank(X) ≤ K.


which is to find a low-rank matrix to approximate 1+1 ρ X̂ + ρ2 (Y k − U k ) and
2
can be solved by truncated SVD.
The optimization problem in Y update is equivalent to
min
Y

s.t.

ρ k+1
kX
− Y + Uk k2F
2
g(Y) ≤ 0.

The problem is convex and it can be viewed as a projection of Xk+1 + Uk
onto the convex set S = {Y|g(Y) ≤ 0}. It can be solved efficiently or even have
closed-form solutions.
14

3

Research Proposal

ALGORITHM DESIGN

Figure 3: An illustration of X update
3.4.2

Convergence analysis

The theoretical analysis is also based on the assumption that the dual variable
will converge. With the same argument in Section 3.3, the value of primal
variable is iteratively updated by the function Y = H(X), where
H(X) = argmin kY − (αX + (1 − α)D) k2F ,
rank(Y)≤K

α=

ρ 4
,
2+ρ

i.e., X k+1 = H(X k ) An interpretation of the update is provide in Fig 3.
If the feasible approximation region is convex, the update rule will converge
with kXk+2 −Xk+1 kF ≤ αkXk+1 −Xk kF , but unfortunately, the set {X|rank ≤
K} is non-convex and the update is not guaranteed to converge. In the following,
we provide some simple necessary conditions for the convergence of the primal
variable. Firstly, if Xk converges to X∗ , then X∗ is a fixed point of H(X) 5 , as
shown in Proposition 1.
Proposition 1 . For any sequence Xk converging to X∗ we will have
X∗ = H(X∗ ).
We can see that if Xk converges to X∗ , then X∗ is a fixed point of Y = H(X).
We define a set F as
F = {X|X = H(X), g(X) ≤ 0},
4 To make function H(X) well-defined, we assume that the minimizer of the X update
is unique, which means that, in each
iteration, σK 6= σK+1 holds for the singular value

decomposition of αXk + (1 − α)D . This assumption is used in the sequel.
5 This result depends on the continuity of H(X) at X∗ , which requires a formal proof. We
conjecture that H(X) is continuous at the points where it is well-defined.

15

3

Research Proposal

ALGORITHM DESIGN

and provide a necessary condition for the convergence of Xk in Lemma 4.
Lemma 4 Under the condition that the dual variable Uk converges to Ū, if Xk
converges to X∗ , then X∗ ∈ F .
Lemma 4 is useful in the sense that, with the knowledge of D, we can restrict
our attention to F. On one hand, if we observe that the dual variable converges,
we do not need to keep updating since the only possible limiting points are in
F ; on the other hand, if we find F is empty, then we can say that the primal
variable will not converge and there is no need to spend more effort on updating.
3.4.2.1

the Fixed Points of H(X)

It turns out that the fixed points of Y = H(X) are highly related to the singular
value decomposition of D, which we will show in this section.
To warm up, we firstly provide one fixed point in Corollary 1.
Corollary 1 If X̃ is the optimal solution of

min

rank(X)≤K

kX − DkF by truncated

SVD, then X̃ = H(X̃).
Proof
Let D = U ΣV T , we can obtain that the optimal solution is X̃ =

Σ1:K 0
U
V T . Then,
0
0


Σ1:K 0
αX̃ + (1 − α)D = αU
V T + (1 − α)U ΣV T
0
0


Σ1:K
0
V T,
=U
0
(1 − α)ΣK+1:N
which is the singular value
 decomposition of αX̃ + (1 − α)D. Then the value of
Σ1:K 0
∗
H(X ) is U
V T.
0
0
The proof is completed.
We note that even though X̃ is a fixed point of Y = H(X), it may not be a
stationary point of X update because X̃ may not satisfy g(X) ≤ 0.
Next, we characterize all possible fixed points of H(X) in Corollary 2.
Corollary 2 Suppose D admits the singular value decomposition of U ΣV T with
σ1 ≥ σ2 ≥ ... ≥ σn ≥ 0, and I = {i1 , i2 , ..., iK } is a subset of {1, 2, ..., n}
with size K. If for any i ∈ I and j ∈
/ I we can have σi ≥ (1 − α)σj , then
P
X̃ = i∈I σi ui viT is a fixed point of Y = H(X). On the other hand, if X̃ is a
P
fixed point, then there must exists such a subset I, such that X̃ = i∈I σi ui viT .

Proof Suppose X̃ is a fixed point of H(X) and αX̃ + (1 − α)D admits
the

Σ
0
1:K
singular value decomposition of UΣVT . Then X̃ = U
VT . By
0
0
16

4 SIMULATIONS

Research Proposal

substituting into (1 − α)D + αX̃ = UΣVT we can have


Σ1:K 0
(1 − α)D + αU
VT = UΣVT .
0
0
Thus
D=U



Σ1:K
0

0
ΣK+1:n
1−α



VT .

However, to have the standard
 SVD of D, we need to permutate the diagonal
ΣK
0
elements of
to make them decreasing and the corresponding
ΣK+1:n
0
1−α
columns of U, V.


Σ1:K
0
are decreasing, which
However, the diagonal elements of
0
ΣK+1:n
establishes the results in this corollary.
We provide some remarks regarding this corollary.
• By Corollary 2, it is easy to see that Corollary 1 is true.

n
• The number of fixed points is upper bounded by K
because we need to
choose K singular values to form one fixed point of C(X).
• A larger value of α (meaning a larger ρ) will make the condition σi ≥
(1 − α)σj more difficult to satisfy and the possible fixed points will be
fewer.

4

Simulations

In this section, we conduct some simulations to show the performance of the
proposed algorithm. Our purpose is to investigate (i) the performance comparison with other algorithms under different scenarios, (ii) the impact of values of
ρ, and (iii) the performance on a real application.

4.1

Non-negative low-rank approximation

Firstly, we consider a special case with the convex constraints being Xi,j ≥
0, ∀i, j, i.e., the non-negative low-rank approximation problem.
Other than the algorithm we propose, there are two alternatives to solve
this problem: Alternating Direction Projection (ADP) [5], and Non-negative
Matrix Factorization (NMF) [20]. We briefly discuss them here.
In each iteration of ADP, we first make a projection onto the rank-K set and
then another projection onto the convex set {X|g(X) ≤ 0}, both of which can
be solved efficiently. One advantage of ADP is that the algorithm is guaranteed
to converge, but this is also its disadvantage, because it is easily to be trapped
by a local optimal. In NMF, we firstly equivalently transform the rank-K
constraints to the fact that X = AB, A ∈ RM×K , B ∈ RK×N and impose
17

4 SIMULATIONS

Research Proposal

93

90

92

89

96

ADMM−RLRA
ADP
NNMF

94

90

Objective values

ADMM−RLRA
ADP
NNMF

89
88

88
87

Objective values

91
Objective values

98

ADMM−RLRA
ADP
NNMF

86

92
90
88
86

85

87

84

86

84

85
0

83
0

5

10
15
20
Number of iterations

25

30

82
5

10
15
20
Number of iterations

25

30

80
0

5

10
15
20
Number of iterations

25

30

Figure 4: Non-negative Figure 5: Non-negative Figure 6: Non-negative
low-rank approximation low-rank approximation low-rank approximation
with K = 3
with K = 6
with K = 10
constraints Ai,j ≥ 0, Bi,j ≥ 0 to ensure that Xi,j ≥ 0, which is obviously
more strict. Then we apply the ADMM algorithm propose in [29] to obtain the
factorization.
In the experiments, the original data X̂ ∈ R100×80 is randomly generate and
ρ = 5 for the algorithm we propose. The objective values of different algorithms
in different scenarios (different rank constraints) are shown in from Fig 4 to
Fig 6.
It can be seen that the performance of our proposed algorithm is always the
best in three cases. NMF fails to converge when K = 10. Another observation
is that the objective value of AD will remain the same after several iterations,
the reason of which is that the optimization variable will not change once it
is feasible. The reason that our algorithm outperforms NMF is that NMF
requires A and B are both non-negative, which is more restrictive than the
original problem.
4.1.1

The impact of ρ

We have little theoretical results on how the parameter ρ will impact the performance and how to choose a proper value of ρ for different problems. We provide
some simulation results regarding to the convergence and performance here. In
this simulation, we fix the data matrix X̂, keep K = 5 and examine the residual
value kX − YkF and objective value for ρ = 1, 5, 9, 15. The result is shown in
Fig 7 8.
From Fig 7, we can see that the residual value converges to 0 when ρ = 5, 9, 15
and converges faster with a larger ρ, but fails to converge when ρ = 1; on the
other hand, as shown in Fig 8, for the scenarios when the algorithm converges,
the objective value decreases more slowly if ρ is larger. As we can see, a larger
ρ will make the algorithm more easily to converge, but will result in a worse
performance, indicating that there is a underlying tradeoff, which we leave for
future work.

18

4 SIMULATIONS

Research Proposal

2

90
ρ=1
ρ=5
ρ=9
ρ=15

1

10

0

10

−1

10

0

10
20
Number of iterations

30

Objective values

)QRUPRIUHVLGXDO

10

ρ=1
ρ=5
ρ=9
ρ=15

88

86

84

82
0

10
20
Number of iterations

30

Figure 7: Values of residual with Figure 8: Values of objective with
different ρ
different ρ

4.2

Image denoising by partial prior knowledge

Next, we consider an application of image denoising. More precisely, we have
a noisy image X̂ and two kinds of prior knowledge: (a) the original image is of
low rank and the rank upper bound is K; (b) the original values of some pixels.
And we want to obtain a clear image with less noise.
Singular value decomposition is shown useful in image denoising [22] [13].
But direct truncated SVD cannot use the second kind of prior knowledge in the
scenario we consider. Here we formulate an optimization problem in the form
of RLRA with the convex constraint g(X) ≤ 0 being Xi,j = xi,j , (i, j) ∈ P,
where P is the set of pixels with known values. We do the experiment on MIT
logo (the rank of the original image is 5) with ρ = 5 and the values of 5% of the
pixels are known as a prior. The truncated SVD algorithm (TSVD) is used as
a comparison. The original image, noised image, denoised image by ADMMRLRA and denoised image by TSVD are shown from Fig 9 to Fig 12. The
image qualities measured by PSNR and SNR are shown in Table 1.
Table 1: Image quality, PSNR, SNR
PSNR
SNR
Noisy Image -3.504
-4.895
By TSVD
8.675
7.284
By RLRA
12.895 11.504
The simulation results shows that the second kind of prior knowledge is useful
and ADMM-RLRA can use that to almost double the PSNR and improve the
image quality significantly, which verifies the effectiveness of our algorithm.

19

Research Proposal

5

CONCLUSION AND FUTURE WORK

Figure 9: Original image: MIT logo Figure 10: Noisy image with Gaussian noise

Figure 11:
RLRA

5

Recovered image by Figure 12:
TSVD

Recovered image by

Conclusion and Future Work

This paper proposes to use ADMM algorithm to solve the restricted low-rank
approximation problem. We show that all subproblems in ADMM can be solved
efficiently. And more interestingly, under the condition that the dual variable
converges, the update of primal variable is carried out by a function and the
possible limiting points are its stationary points,. We provide some preliminary
theoretical results based on this understanding and show its performance by
experiments.
The theoretical part in this paper is limited because most of them is based on
the assumption that the dual variable converges, which is not guaranteed. The
most important future work is to study under what condition the dual variable
will converge. From Corollary 2 and experiments in Section 4.1.1, we can see
that the parameter ρ plays an important role and an interesting question is to
ask how the value of ρ will affect the convergence of this algorithm and how to
choose a proper ρ.
There can be several different ways to reformulate the original problem and
apply ADMM. As we mention in Section 3.4, we can change the update order.
And also, we can let I(X) be the indicator function for the convex constraints
and J (Y ) be the indicator function for the rank constraint. It is interesting
to study how the different approaches will affect the solution theoretically and

20

REFERENCES

Research Proposal

empirically.

References
[1] A. Berman and N. Shaked-Monderer. Completely positive matrices. World
Scientific, 2003.
[2] P. Bofill and M. Zibulevsky. Underdetermined blind source separation using
sparse representations. Signal processing, 81(11):2353–2362, 2001.
[3] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed
optimization and statistical learning via the alternating direction method
of multipliers. Foundations and Trends R in Machine Learning, 3(1):1–122,
2011.
[4] R. Chartrand. Nonconvex splitting for regularized low-rank+ sparse decomposition. Signal Processing, IEEE Transactions on, 60(11):5810–5819,
2012.
[5] M. T. Chu, R. E. Funderlic, and R. J. Plemmons. Structured low rank
approximation. Linear algebra and its applications, 366:157–172, 2003.
[6] D. L. Donoho. Compressed sensing. Information Theory, IEEE Transactions on, 52(4):1289–1306, 2006.
[7] C. Eckart and G. Young. The approximation of one matrix by another of
lower rank. Psychometrika, 1(3):211–218, 1936.
[8] M. Fazel. Matrix rank minimization with applications. PhD thesis, PhD
thesis, Stanford University, 2002.
[9] M. Fazel, H. Hindi, and S. P. Boyd. Log-det heuristic for matrix rank
minimization with applications to hankel and euclidean distance matrices.
In American Control Conference, 2003. Proceedings of the 2003, volume 3,
pages 2156–2162. IEEE, 2003.
[10] M. Grant and S. Boyd. Cvx: Matlab software for disciplined convex programming.
[11] J.-B. Hiriart-Urruty and C. Lemaréchal. Convex analysis and minimization algorithms I: fundamentals, volume 305. Springer Science & Business
Media, 2013.
[12] M. Hong, Z.-Q. Luo, and M. Razaviyayn. Convergence analysis of admm
for a family of nonconvex problems.
[13] Z. Hou. Adaptive singular value decomposition in wavelet domain for image
denoising. Pattern Recognition, 36(8):1747–1763, 2003.

21

REFERENCES

Research Proposal

[14] B. Jiang, S. Ma, and S. Zhang. Alternating direction method of multipliers for real and complex polynomial optimization models. Optimization,
63(6):883–898, 2014.
[15] J. Lavaei and S. H. Low. Zero duality gap in optimal power flow problem.
Power Systems, IEEE Transactions on, 27(1):92–107, 2012.
[16] H. Lee, R. Grosse, R. Ranganath, and A. Y. Ng. Convolutional deep belief
networks for scalable unsupervised learning of hierarchical representations.
In Proceedings of the 26th Annual International Conference on Machine
Learning, pages 609–616. ACM, 2009.
[17] S. H. Low. Convex relaxation of optimal power flow: a tutorial. In Bulk
Power System Dynamics and Control-IX Optimization, Security and Control of the Emerging Power Grid (IREP), 2013 IREP Symposium, pages
1–15. IEEE, 2013.
[18] Z.-Q. Luo, W.-K. Ma, A.-C. So, Y. Ye, and S. Zhang. Semidefinite relaxation of quadratic optimization problems. Signal Processing Magazine,
IEEE, 27(3):20–34, 2010.
[19] S. Magnússon, P. C. Weeraddana, M. G. Rabbat, and C. Fischione. On
the convergence of alternating direction lagrangian methods for nonconvex
structured optimization problems. arXiv preprint arXiv:1409.8033, 2014.
[20] I. Markovsky. Structured low-rank approximation and its applications.
Automatica, 44(4):891–909, 2008.
[21] D. Needell and J. A. Tropp. Cosamp: Iterative signal recovery from incomplete and inaccurate samples. Applied and Computational Harmonic
Analysis, 26(3):301–321, 2009.
[22] A. Rajwade, A. Rangarajan, and A. Banerjee. Image denoising using the
higher order singular value decomposition. Pattern Analysis and Machine
Intelligence, IEEE Transactions on, 35(4):849–862, 2013.
[23] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions
of linear matrix equations via nuclear norm minimization. SIAM review,
52(3):471–501, 2010.
[24] B. Schölkopf and A. J. Smola. Learning with kernels: Support vector machines, regularization, optimization, and beyond. MIT press, 2002.
[25] J.-L. Starck and J. Bobin. Astronomical data analysis and sparsity: from
wavelets to compressed sensing. Proceedings of the IEEE, 98(6):1021–1030,
2010.
[26] M. Unser. Wavelets, sparsity and biomedical image reconstruction.

22

REFERENCES

Research Proposal

[27] Y. Xu, W. Yin, Z. Wen, and Y. Zhang. An alternating direction algorithm
for matrix completion with nonnegative factors. Frontiers of Mathematics
in China, 7(2):365–384, 2012.
[28] S. You and Q. Peng. A non-convex alternating direction method of multipliers heuristic for optimal power flow. In Smart Grid Communications
(SmartGridComm), 2014 IEEE International Conference on, pages 788–
793. IEEE, 2014.
[29] Y. Zhang. An alternating direction algorithm for nonnegative matrix factorization. preprint, 2010.
[30] H. Zhu and G. Giannakis. Estimating the state of ac power systems using
semidefinite programming. In North American Power Symposium (NAPS),
2011, pages 1–7. IEEE, 2011.

23

