Download Report

On how to solve large-scale log-determinant
optimization problems
Chengjing Wang
∗
January 25, 2014
Abstract
We propose a proximal augmented Lagrangian method and a hybrid method, i.e.,
employing the proximal augmented Lagrangian method to generate a good initial
point and then employing the Newton-CG augmented Lagrangian method to get
a highly accurate solution, to solve large-scale nonlinear semideﬁnite programming
problems whose objective functions are a sum of a convex quadratic function and
a log-determinant term. We demonstrate that the algorithms can supply a high
quality solution eﬃciently even for some ill-conditioned problems.
Keywords: Quadratic programming, Log-determinant optimization problem, Proximal augmented Lagrangian method, Augmented Lagrangian method, Newton-CG method
1
Introduction
In this paper, by deﬁning log 0 := −∞, we consider the following standard primal and
dual nonlinear semideﬁnite programming problems whose objective functions are a sum
of a convex quadratic function and a log-determinant term (QP-Logdet) :
(P )
1
min{ ⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X : A(X) = b, X ≽ 0},
X
2
(D)
1
max {− ⟨X, Q(X)⟩ + bT y + µ log det Z + nµ(1 − log µ) :
X,y,Z
2
−Q(X) + A∗ y + Z = C, Z ≽ 0},
∗
School of Mathematics and Institute of Applied Mathematics, Southwest Jiaotong University, No.999,
Xian Road, Chengdu 611756, China ([email protected]). The author’s research was supported by
the National Natural Science Foundation of China under grant 11201382, the Youth Fund of Humanities
and Social Sciences of the Ministry of Education under grant 12YJC910008, the project of the science
and technology department of Sichuan province under grant 2012ZR0154, and the Fundamental Research
Funds for the Central Universities under grants SWJTU12CX055 and SWJTU12ZT15.
1
where Q : S n → S n is a given self-adjoint positive semideﬁnite linear operator, C ∈ S n ,
b ∈ Rm , µ ≥ 0 is a given parameter, A : S n → Rm is a linear map and A∗ : Rm → S n is
the adjoint of A. Note that the linear maps A and A∗ can be expressed, respectively, as
[
A(X) = ⟨A1 , X⟩, . . . , ⟨Am , X⟩
]T
,
∗
A (y) =
m
∑
yk Ak ,
(1)
k=1
where Ak , k = 1, · · · , m are given matrices in S n . As for the explanation of all other
main notations one may see Subsection 1.1.
The perturbed QP-Logdet problem has the form
{1
}
(P ε )
min
⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X : A(X) = b, X ≽ εI .
X
2
(Dε )
1
max {− ⟨X, Q(X)⟩ + bT y + ⟨Z, εI⟩ + µ log det(Q(X) + C − A∗ y − Z)
X,y,Z
2
+nµ(1 − log µ) : −Q(X) + A∗ y + Z + µX −1 = C, Z ≽ 0}.
The QP-Logdet problem (P ) is itself a classical convex model problem in optimization
theory. It can be regarded as an extension of the qudratic semideﬁnite programming
problem (QSDP) and the log-determinant (Logdet) problem, so it shares the structures of
both problems, and it goes without saying that the QP-Logdet problem is considerable.
For the QSDP, it is certainly a heart problem in nonlinear semideﬁnite programming
problems, which has been considered by Toh [35], Toh, T¨
ut¨
unc¨
u and Todd [36, 37], Zhao
[45], Jiang, Sun and Toh [14], etc.. For the Logdet problem, which has a very important
application in covariance selection [5] and has been intensively studied over the past several
years, including the work of Dahl, Vandenberghe and Roychowdhury [4], d’Aspremont,
Banerjee and El Ghaoui [6], Li and Toh [15], Lu [16, 17], Lu and Zhang [18], Olsen,
Oztoprak, Nocedal and Rennie [24], Scheinberg, Ma and Goldfarb [30], Scheinberg and
Rish [31], Toh [34], Wang, Sun and Toh [40], Yang, Sun and Toh [41], Yang, Shen, Wonka,
Lu and Ye [43], Yuan [44], etc.. As far as the QP-Logdet problem be concerned, it also
arises in many practical applications such as robust simulation of global warming policies
[13], speech recognition [39], and so on. Thus the algorithms developed to solve this kind
of problems can potentially ﬁnd wide applications.
For the QP-Logdet problems of small and medium size, the interior-point method
(IPM) with a direct solver is certainly an eﬃcient and robust approach; however, for these
large-scale problems, the IPM lacks the ability due to the need of computing, storing, and
factorizing the Schur matrices that are typically dense. In view of this, we need to design
new approaches.
The proximal augmented Lagrangian (PAL) method proposed in [7] is a fast primaldual approach. The key idea of this approach is to apply the proximal technique to
the augmented Lagrangian function of the primal problem at every iteration so that the
2
simpliﬁed inner problem has an analytical or at least a semi-analytical solution which can
be solved very fast, and then to update the dual variables. The convergence of the PAL
method has been shown in [7]. The biggest advantage of the PAL method is that there
is no need to solve any linear system of equations to obtain step directions. It is a good
approach for generating a good initial point. Furthermore, for some ill-conditioned inner
problems whose Hessian matrices are of near low ranks, the proximal technique of the
PAL method is actually better than the Newton-CG method since it can be regarded as a
regularization treatment in some sense. Allowing for the advantages of the PAL method,
we may apply it to the QP-Logdet problem, especially for generating an initial point.
Although the PAL method is a good choice for the QP-Logdet problem, it is after all a
gradient method, once the iteration point is near the optimal solution, we may consider an
approach with locally faster convergence rate. The augmented Lagrangian (AL) method
[29] is just in the position to play this role. It is a classical method that can be viewed as a
special form of the proximal point algorithm (PPA). (As for the PPA, it was proposed by
Martinet [19], and further studied by Rockafellar [28, 29].) Although the AL method for
convex programs is a gradient ascent method applied to the corresponding dual problems
[27], Sun, Sun and Zhang revealed that under the constraint nondegenerate conditions
[1], it could be locally regarded as an approximate generalized Newton method applied
to a semismooth equation. It is just this reason that inspired us to apply the AL method
to solve the QP-Logdet problem. Great successes of the applications of the AL method
to large-scale semideﬁnite programming problems can also be seen in [46, 40, 41] and the
references therein.
If applying the AL method to problem (P), we have to solve the inner problems
which have no closed-form solutions. The well-tested quadratically convergent semismooth
Newton method introduced by Qi and Sun [25] is an ideal approach to fulﬁll this task.
As for the semismooth Newton direction, the preconditioned conjugate gradient (PCG)
method is a good choice because the Newton system of equations is large and a direct
solver is not proper.
Mainly based on Rockafellar’s theoretical framework the global convergence and local
convergence rate of the sequence generated by the Newton-CG AL (NAL) method can be
established.
The numerical results demonstrate that the proposed approach can be very eﬃcient and robust to supply highly accurate solutions for large-scale problems not only for
synthetic problems but also for real data problems. In this paper, we solve synthetic
QP-Logdet problems with n up to 2, 000 and m up to 1, 186, 173 in about 2 hours.
The remaining parts of this paper are organized as follows. In Section 2, we give a
brief introduction on some preliminary contents. In Section 3, we describe the algorithms
in details. In Section 4, we establish the convergence theory of the proposed approaches.
In Section 5 we describe some numerical issues on the semismooth Newton-CG method.
In Section 6, we present the numerical results. Finally, we give some concluding remarks.
3
1.1
Additional Notations
In this paper, all vectors are assumed to be ﬁnite dimensional. The symbol Rn denotes
the n-dimensional Euclidean space, and On denotes the set of n × n orthogonal matrices.
The set of all m × n matrices with real entries is denoted by Rm×n , ⟨·, ·⟩ stands for
the standard trace inner product in S n , ∥ · ∥ denotes the induced Frobenius norm, S+n
n
and S++
denote the sets of positive semideﬁnite matrices and positive deﬁnite matrices,
respectively. The symbol ◦ denotes the Hadamard product, and ⊗ denotes the Kronecker
product. The symbol FP ε := {X ∈ S n : A(X) = b, X ≽ εI} denotes the feasible set of
the problem (P ε ), and ΠK (·) denotes the metric projector onto the closed convex set K.
Let vec : Rm×n → Rmn be the vectorization operator on matrices deﬁned by stacking the
columns of a matrix to a long vector one by one.
2
Preliminaries
In order to be able to present our ideas more clearly, we ﬁrst introduce some concepts
related to the AL method based on the classical papers by Rockafellar [28, 29].
2.1
Maximal monotone operators
Let H be a real Hilbert space with an inner product ⟨·, ·⟩. A multifunction T : H ⇒ H is
said to be a monotone operator if
⟨z − z ′ , w − w′ ⟩ ≥ 0, whenever w ∈ T (z), w′ ∈ T (z ′ ).
(2)
It is said to be maximal monotone if, in addition, the graph
G(T ) = {(z, w) ∈ H × H| w ∈ T (z)}
is not properly contained in the graph of any other monotone operator T ′ : H ⇒ H.
For example, if T is the subdiﬀerential ∂f of a lower semicontinuous convex function
f : H → (−∞, +∞], f ̸≡ +∞, then T is maximal monotone (see Minty [21] or Moreau
[22]), and the relation 0 ∈ T (z) means that f (z) = min f .
2.2
Representations in terms of maximal monotone operators
The following Karush-Kuhn-Tucker (KKT) conditions are necessary and suﬃcient for the
optimality of (P ) and (D):
A(X) − b = 0,
Q(X) + C − A∗ y − Z = 0,
XZ = µI, X ≽ 0, Z ≽ 0.
4
The following KKT conditions are also necessary and suﬃcient for the optimality of (P ε )
and (Dε ):
A(X) − b = 0,
Q(X) + C − A∗ y − Z − µX −1 = 0,
⟨X − εI, Z⟩ = 0, X − εI ≽ 0, Z ≽ 0.
Throughout this paper, the following conditions for (P ε ) and (Dε ) are assumed to
hold.
Assumption 2.1. Problem (P ε ) satisﬁes the condition
n
∃ X0 ∈ S++
such that A(X0 ) = b, X0 ≻ εI.
Assumption 2.2. For problem (P ε ), there exists an α ∈ R such that the level set
{x| f0 (X) ≤ α, A(X) = b, X ≽ εI} is nonempty and bounded, where f0 (x) stands for
the objective function.
Let l(X; y, Z) : S n × Rm × S n → R
the extended form:


 L0 (X; y, Z)
−∞
l(X; y, Z) =


+∞
be the ordinary Lagrangian function for (P ε ) in
if X ∈ FP ε and (y, Z) ∈ Rm × S+n ,
if X ∈ FP ε and (y, Z) ̸∈ Rm × S+n ,
if X ̸∈ FP ε ,
where
L0 (X; y, Z) =
1
⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X − ⟨y, A(X) − b⟩ − ⟨Z, X − εI⟩.
2
The essential objective function in (P ε ) is given by
{ 1
⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X if X ∈ FP ε ,
2
f (X) =
sup
l(X; y, Z) =
+∞
otherwise,
y∈Rm ,Z∈S n
while the essential objective function in (Dε ) is given by
{
g0 (y, Z) if (y, Z) ∈ Rm × S+n ,
g(y, Z) = inf n l(X; y, Z) =
X∈S
−∞
if (y, Z) ̸∈ Rm × S+n ,
5
As in [29], we can deﬁne the following operators
Tl (X, y, Z) = {(V, u, W ) ∈ S n × Rm × S n | (V, −u, −W ) ∈ ∂l(X; y, Z)},
(X, y, Z) ∈ S n × Rm × S n ,
Tf (X) = {V ∈ S n | V ∈ ∂f (X)}, X ∈ S n ,
Tg (y, Z) = {(u, W ) ∈ Rm | (−u, −W ) ∈ ∂g(y, Z)}, (y, Z) ∈ Rm × S n .
We can observe that l is a closed proper saddle function in the sense of [26, p.363], and the
mapping Tl is a maximal monotone operator in S n ×Rm ×S n [26, Cor.37.5.2]. Meanwhile,
f is a closed proper convex function, so Tf is a maximal monotone operator in S n [21, 22].
Similarly, g is a closed proper concave function, so Tg is a maximal monotone operator in
Rm × S n .
To discuss the rate of convergence, we introduce some related concepts and conclusions.
Deﬁnition 2.1. (cf. [28]) For a maximal monotone operator T , we say that its inverse
T −1 is Lipschitz continuous at the origin (with modulus a ≥ 0) if there is a unique solution
z¯ to z = T −1 (0), and for some τ > 0 we have
∥z − z¯∥ ≤ a∥w∥, where z ∈ T −1 (w) and ∥w∥ ≤ τ.
Assumption 2.3. Let X be the optimal solution to problem (P ε ) and (y, Z) is the corresponding Lagrangian multiplier. We say that the primal constraint nondegeneracy [1]
holds at X to problem (P ε ) if
(
)
(
) ( m )
{0}
A
R
n
S +
=
,
n
lin(TS+ (X − εI))
I
Sn
where lin(TS+n (X − εI)) denotes the lineality space of TS+n (X − εI), and TS+n (X − εI)
denotes the tangent cone of S+n at X − εI.
Assumption 2.4. Let X be the optimal solution to problem (P ε ) and (y, Z) is the corresponding Lagrangian multiplier. We say that the strong second order suﬃcient condition
[3] holds at X if
⟨H, ∇2XX L0 (X; y, Z)(H)⟩ + Υ(X−εI) (Z, H) > 0,
∀ H ∈ aﬀ(C(X)) \ {0},
where the linear-quadratic function ΥB : S n × S n → R is deﬁned by
ΥB (Γ, A) := 2⟨Γ, AB + A⟩, (Γ, A) ∈ S n × S n ,
where B + is the Moore-Penrose pseudo-inverse of B.
6
Next we present the direct condition for the Lipschitz continuity of Tg−1 due to Rockafellar [29, Prop.2].
Proposition 2.1. If the strong second-order suﬃcient condition and the primal constraint
nondegeneracy condition are satisﬁed, then Tl−1 is actually single-valued and continuously
diﬀerentiable on a neighborhood of the origin, and so are Tf−1 and Tg−1 . Thus in particular,
these mappings are all Lipschitz continuous at the origin.
Remark 2.1. We must emphasize that the so called strong second-order condition in [29]
actually includes the primal constraint nondegeneracy condition.
3
3.1
The algorithms
The PAL method
In this subsection, we describe the PAL method in details. The augmented Lagrangian
function of (P) is
Lσ (X; y) =
1
σ
⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X − ⟨y, A(X) − b⟩ + ∥A(X) − b∥2 , X ≽ 0.
2
2
At the kth iteration, we solve the following subproblem
}
{
1
argmin Lσk (X; y k ) + ∥X − X k ∥2Tk | X ≽ 0
2
{1
}
= argmin ⟨X, (Q + σk A∗ A + Tk )(X)⟩ − ⟨M k , X⟩ − µ log det X | X ≽ 0 , (3)
2
where M k = A∗ y k + σk A∗ b + Tk (X k ) − C. Pick Tk ≽ 0 such that Q + σk A∗ A + Tk = αk I,
with αk = λmax (Q + σk A∗ A), then M k = −Q(X k ) − C + αk X k + A∗ (y k + σk (b − A(X k ))),
and (3) simpliﬁes to
{1
}
1 k 2
argmin αk ∥X − M ∥ − µ log det X | X ≽ 0 .
(4)
2
αk
Based on Lemma 2.1 in [40], problem (4) has a closed-form solution
X k+1 = ϕ+
γk (
1 k
k
k T
M ) := P k diag(ϕ+
γk (d ))(P ) ,
αk
(5)
Z k+1 = ϕ−
γk (
1 k
k
k T
M ) := P k diag(ϕ−
γk (d ))(P ) ,
αk
(6)
and further
where α1k M k = P k Dk (P k )T is the eigenvalue decomposition of α1k M k with P k ∈ On and
1
k
k
− 1
k
+
Dk = diag(dk ), dk ∈ Rn , ϕ+
γk ( αk M ) and ϕγk ( αk M ) are matrix valued functions, ϕγk (d )
7
√
k
ϕ−
γk (d )
ϕ+
γk (x)
x2 +4γ
k
are vector valued functions, whose scalar counterparts are
=
2
√
x2 −4γk
and ϕ−
, correspondingly, and γk = αµk . From (5) and (6), we can see
γk (x) =
2
that X k+1 and Z k+1 are positive deﬁnite automatically. Once X k+1 is obtained, one may
update the dual variable to get y k+1 .
As for λmax (Q + σk A∗ A), we may apply the power method to compute it.
Now we may summarize the PAL method as below
and
n
Algorithm 1: The PAL method. Input X 0 ∈ S++
, y 0 ∈ Rm , σ0 = 10. Iterate the
following steps:
Step 1. Compute αk = λmax (Q + σk A∗ A).
Step 2. Set Tk = αk I − Q − σk A∗ A, and M k = A∗ y k + σk A∗ b + Tk (X k ) − C. Compute
X k+1 = ϕ+
γk (
1 k
1 k
M ) and Z k+1 = ϕ−
M ).
γk (
αk
αk
Step 3. Set
y k+1 = y k − τ σk (A(X k+1 ) − b),
where τ ∈ [1, 2) is a given parameter.
Step 4. Compute
RPk+1 =
∥b − A(X k+1 )∥
,
max{1, ∥b∥}
k+1
RD
=
∥Q(X k+1 ) + C − A∗ y k+1 − Z k+1 ∥
,
max{1, ∥C∥}
k+1
If max{RPk+1 , RD
} < Tol1, stop; else, choose σk+1 ; end.
3.2
The NAL method
In this subsection, we apply the NAL method to the perturbed problem (P ε ) since this
method can not guarantee the positive deﬁniteness of the variable X automatically.
Given a penalty parameter σ > 0, the augmented Lagrangian function for problem
8
(P ε ) is deﬁned as
Lσ (X; y, Z) =
1
σ
⟨X, Q(X)⟩ + ⟨C, X⟩ − µ log det X − ⟨y, A(X) − b⟩ + ∥A(X) − b∥2 +
2
2
1
1
∥ΠS+n (Z − σ(X − εI))∥2 −
∥Z∥2 .
2σ
2σ
At the kth iteration, we solve the following subproblem
{
}
argmin Lσk (X; y k , Z k ) | X ≽ εI .
(7)
Since the subproblem (7) has no closed-form solution, it needs elaborating on. We can
calculate the ﬁrst-order derivative of Θk (X) := Lσk (X; y k , Z k ) with respect to X
∇Θk (X) = Q(X) + C − A∗ y k − µX −1 + σk A∗ (A(X) − b) − ΠS+n (Z k − σk (X − εI)).
Furthermore, base on [26, Thm. 23.8] and [32, Lem. 2.1] we can calculate the Clarke
generalized Jacobian ∂ 2 Θk (X) := ∂∇Θk (X) as below
∂ 2 Θk (X)[H] = (Q + σk A∗ A)[H] + µX −1 HX −1 + σk ∂ΠS+n (Z k − σk (X − εI)), ∀ H ∈ S n .
Actually, from [20, Prop.1] every element in ∂ΠS+n (·) is positive semideﬁnite, so every
element in ∂ 2 Θk (X) is positive deﬁnite. Therefore we can apply the semismooth Newton
method developed in [25] with the line search technique which guarantees the global
convergence. And we may obtain the Newton direction by solving the linear system of
equations
Vbσ0k [H] = −∇Θk (X k ),
where
Vbσ0k := Q + σk A∗ A + µ(X k )−1 ⊗ (X k )−1 + Vσ0k ,
here Vσ0k ∈ ∂ΠS+n (Z k − σk (X k − εI)) may be chosen in the same way as that in [46]
Vσ0k [H] = σk P k (Ωk ◦ ((P k )T (H)P k ))(P k )T ,
where P k ∈ On and
Z k − σk (X k − εI) = P k Λk (P k )T ,
where Λk is the diagonal matrix with diagonal entries consisting of the eigenvalues λk1 ≥
λk2 ≥ · · · ≥ λkn of Z k − σk (X k − εI); if deﬁning three index sets
αk := {i | λki > 0}, βk := {i | λki = 0}, and γk := {i | λki < 0},
[
then
k
Ω =
Eγ¯k γ¯k νγ¯k γk
0
νγ¯Tk γk
]
, νij :=
9
λki
, i ∈ γ¯k , j ∈ γk ,
λki − λkj
here γ¯k = {1, · · · , n} \ γk , and Eγ¯k γ¯k ∈ S |¯γk | is the matrix of ones.
Now we summarize the algorithm as below
Algorithm 2: The NAL method. Input X 0 ∈ FP ε , y 0 ∈ Rm , Z 0 ∈ S+n , σ0 = 10.
Iterate the following steps:
Step 1. Find an approximate minimizer
X k+1 ≈ arg minn Θk (X).
X∈S++
(8)
Step 2. Compute
y k+1 = y k − σk (A(X k+1 ) − b),
Z k+1 = ΠS+n (Z k − σk (X k+1 − εI)).
Step 3. Compute
RPk+1 =
∥b − A(X k+1 )∥
,
max{1, ∥b∥}
k+1
RD
=
∥Q(X k+1 ) + C − A∗ y k+1 − Z k+1 − µ(X k+1 )−1 ∥
,
max{1, ∥C∥}
k+1
If max{RPk+1 , RD
} < Tol2, stop; else, σk+1 = 2σk ; end.
In Algorithm 2, the principal computing costs lie in Step 1, i.e., computing the approximate minimizer of the inner problem (8). So we state the algorithm about how to
compute the approximate minimizer in details below.
10
Algorithm 3: The semismooth Newton-CG Method.
Step 1. Given ζ ∈ (0, 12 ), τ1 , τ2 ∈ (0, 1), and δ ∈ (0, 1), choose X 0 (≽ εI).
Step 2. For j = 0, 1, 2, . . . ,
Step 1.1. Apply the PCG method to ﬁnd an approximate solution H j to
(Vbσ0k + ϵj I)[H] = −∇Θk (X j ),
(9)
where ϵj := τ1 min{τ2 , ∥∇Θk (X j )∥}.
Step 1.2. Set αj = δ mj , where mj is the ﬁrst nonnegative integer m for which
Θk (X j + δ m H j ) ≤ Θk (X j ) + ζδ m ⟨∇Θk (X j ), H j ⟩.
Step 1.3. Set
X j+1 = X j + αj H j .
From Algorithm 3, we can see that the main computing costs lie in solving the linear
operator equation (9). The linear operator corresponds to a fourth-order n dimemsional
tensor or a n2 × n2 matrix. For large-scale problems, n is very large, so direct solvers are
not suitable and PCG solvers are ideal approaches for (9). Among various PCG solvers,
we adopt the symmetric QMR algorithm [9]. But we must note that comparing with
the algorithm in [9], we change the matrix to the linear operator in (9), and change the
vectors to the corresponding matrices. As these changes are only a trivial generalization
of the symmetric QMR algorithm, we do not go to details.
4
Convergence analyses
The convergence analyses of these two methods can be derived from Fazel, Pong, Sun and
Tseng’s paper [7] and Rockafellar’s paper [28, 29] without many diﬃculties, respectively.
For the sake of completeness, we also present these results below.
4.1
Convergence analysis for the PAL method
Theorem 4.1. Assume that the solution set of (P) is nonempty and Assumption 2.1
holds. Assume that Q + σk A∗ A + √Tk is positive deﬁnite. Let {X k , y k , Z k } be generated
from Algorithm 1, and if τ ∈ (0, 1+2 5 ), then the sequence {X k } converges to the optimal
solution to (P) and {y k , Z k } converges to the optimal solution to the dual problem (D).
Proof. Since Q + σk A∗ A + Tk = λmax (Q + σk A∗ A)I ≻ 0, based on Theorem B.1 in [7]
and the KKT condition, the conclusion of the theorem is obvious.
11
4.2
Convergence analysis for the NAL method
Since we can not solve the inner problems exactly, we will use the following stopping
criteria considered by Rockafellar [28, 29] for terminating Algorithms 2 and 3:
(A) Θk (X
k+1
) − inf Θk (X) ≤
ϵ2k /2σk ,
ϵk ≥ 0,
∞
∑
ϵk < ∞;
k=0
(B) Θk (X
k+1
) − inf Θk (X) ≤
δk2 /2σk ∥(y k+1 , Z k+1 )
− (y , Z )∥ , δk ≥ 0,
k
k
2
∞
∑
δk < ∞;
k=0
(B ′ ) ∥∇Θk (X k+1 )∥ ≤ δk′ /σk ∥(y k+1 , Z k+1 ) − (y k , Z k )∥, 0 ≤ δk′ → 0.
We directly obtain from [28, 29] the following convergence results.
Theorem 4.2. Let Algorithm 2 be executed with stopping criterion (A). If Assumption
2.1 is satisﬁed, then the generated sequence {(y k , Z k )} is bounded and {(y k , Z k )} converges
to (y, Z), where (y, Z) is the unique optimal solution to (Dε ), and {X k } is asymptotically
minimizing for (P ε ) with max(Dε ) = inf(P ε ).
The boundedness of {y k , Z k } under (A) is actually equivalent to the existence of an
optimal solution to (Dε ).
If {y k , Z k } is bounded and Assumption 2.2 is satisﬁed, then the sequence {X k } is also
bounded, and the accumulation point of the sequence {X k } is the unique optimal solution
to (P ε ).
Theorem 4.3. Let Algorithm 2 be executed with stopping criteria (A) and (B). If Assumption 2.1 holds and Tg−1 is Lipschitz continuous at the origin with modulus ag , then
{(y k , Z k )} converges to the unique optimal solution (y, Z) with max(Dε ) = inf(P ε ), and
for all k suﬃciently large,
∥(y k+1 , Z k+1 ) − (y, Z)∥ ≤ θk ∥(y k , Z k ) − (y, Z)∥,
where
2 −1/2
θk = [ag (a2g + σk2 )−1/2 + δk ](1 − δk )−1 → θ∞ = ag (a2g + σ∞
)
< 1, σk → σ∞ ,
Moreover, the conclusions of Theorem 4.2 about {(y k , Z k )} are valid.
If in addition to (B) and the condition on Tg−1 one has (B’) and the stronger condition
that Tl−1 is Lipschitz continuous at the origin with modulus al (≥ ag ), then X k → X,
where X is the unique optimal solution to (P ε ), and one has that for all k suﬃciently
large,
∥X k+1 − X∥ ≤ θk′ ∥(y k+1 , Z k+1 ) − (y k , Z k )∥,
where
′
= al /σ∞ .
θk′ = al (1 + δk′ )/σk → θ∞
12
Remark 4.1. In Algorithm 2 we can also add the term 2σ1k ∥X −X k ∥2 to Θk (X). Actually,
in our Matlab code, one can optionally add this term. This actually corresponds to the
proximal method of multipliers considered in [29, Section 5]. Convergence analysis for
this improvement can be conducted in a parallel way as for Algorithm 2.
Note that in the stopping criteria (A) and (B), inf Θk (X) is an unknown value. Since
b
Θk (X) := Θk (X) + 2σ1k ∥X − X k ∥2 is a strongly convex function with modulus σ1k , then
one has the estimation
b k (X k+1 ) − inf Θ
b k (X) ≤ 1 ∥∇Θ
b k (y k+1 )∥2 ,
Θ
2σk
thus criteria (A) and (B) can be practically modiﬁed as follows:
b k (y k+1 )∥ ≤ ϵk , ϵk ≥ 0,
∥∇Θ
b k (X k+1 )∥ ≤ δk ∥(y k+1 , Z k+1 ) − (y k , Z k )∥, δk ≥ 0,
∥∇Θ
∞
∑
k=0
∞
∑
ϵk < ∞;
δk < ∞.
k=0
Remark 4.2. The condition that Tg−1 is Lipschitz continuous at the origin is not easy to
be realized, however, if Assumptions 2.3 and 2.4 are satisﬁed, then Tl−1 , Tf−1 and Tg−1 are
all single-valued and continuously diﬀerentiable on a neighborhood of the origin, and in
particular they are all Lipschitz continuous at the origin.
5
Numerical issues in the associated semismooth NewtonCG method
When applying Algorithm 3 to solve the inner problem (7), the most expensive step lies
in solving the linear system of equations
where
(M + εI)vec(H) = vec(−∇Θk (X)),
(10)
M := Q + σAT A + µX −1 ⊗ X −1 + σP ⊗ P diag(vec(Ω))P ⊗ P,
(11)
Q and A denote the matrix representations of Q and A, to obtain the approximate Newton
direction. In order to solve (10) as eﬃciently as possible, it is necessary to analyze the
condition number of the coeﬃcient matrix and design an eﬃcient preconditioner.
5.1
Conditioning of the coeﬃcient matrix
Let
S = X − εI.
13
For simplicity, we suppose that strict complementarity holds for Z, S, i.e., Z + S ≻ 0.
From the fact that ZS = 0, we have
[ Z
]
Λ
0
PT,
(12)
Z − σS = P
0 −σΛS
where ΛZ = diag(λZ ) ∈ Rr×r and ΛS = diag(λS ) ∈ R(n−r)×(n−r) are diagonal matrices
of positive eigenvalues of Z and S, respectively. Deﬁne the index sets γ¯ := {1, . . . , r},
γ := {r + 1, . . . , n}. Let
]
[
λZi
Eγ¯γ¯ νγ¯γ
,
ν
:=
Ω=
, i ∈ γ¯ , j ∈ γ,
(13)
ij
νγ¯Tγ 0
λZi + σλSj−r
and
c1 =
min(λZ )
max(λZ )
,
c
=
< σ.
2
min(λZ )/σ + max(λS )
max(λZ )/σ + min(λS )
Then we have
Proposition 5.1. Suppose that the strict complementarity holds for Z and S, then we
have the following bound on the condition number of M
κ(M ) ≤
λmax (Q) + σλmax (AT A) + µ/λ2min (X) + σλmax (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )
,
λmin (Q) + σλmin (AT A) + µ/λ2max (X) + c1 λmin (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )
(14)
where Pe1 = Pγ¯ ⊗ Pγ¯ , Pe2 = Pγ ⊗ Pγ¯ , Pe3 = Pγ¯ ⊗ Pγ , Pγ¯ and Pγ are submatrices of
P = [Pγ¯ Pγ ].
Proof. From (11), (12) and (13), we obtain
M = Q + σAT A + µX −1 ⊗ X −1 + σ(Pe1 Pe1T + Pe2 D2 Pe2T + Pe3 D3 Pe3T ),
where D2 = diag(vec(νγ¯γ )), and D3 = diag(vec(νγ¯Tγ )). Since
c1 ≤ σνij ≤ c2 , i ∈ γ¯ , j ∈ γ,
then we can derive
c1 (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T ) ≼ σ(Pe1 Pe1T + Pe2 D2 Pe2T + Pe3 D3 Pe3T ) ≼ σ(Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T ).
Furthermore, from [11, Thm 4.2.12],
µ
λ2max (X)
I ≼ µX −1 ⊗ X −1 ≼
14
µ
λ2min (X)
I.
Hence from [10, Thm 4.3.7],
µ
+ c1 λmin (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )]I ≼ M ≼
2
λmax (X)
µ
+ σλmax (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )]I.
[λmax (Q) + σλmax (AT A) + 2
λmin (X)
[λmin (Q) + σλmin (AT A) +
So we can get the following bound on the condition number of M :
κ(M ) ≤
λmax (Q) + σλmax (AT A) + µ/λ2min (X) + σλmax (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )
.
λmin (Q) + σλmin (AT A) + µ/λ2max (X) + c1 λmin (Pe1 Pe1T + Pe2 Pe2T + Pe3 Pe3T )
The upper bound in (14) suggests that: (i) with σ and 1/c1 increasing, the condition
number κ(M ) may increase. (ii) With Q, AT A and P ⊗ P diag(vec(Ω))P ⊗ P being of
(near) low ranks, µX −1 ⊗ X −1 can potentially improve the condition number, for the
term µ/λ2max (X)I can “lift” the minimal eigenvalue of M ; in other words, if Q, AT A and
P ⊗ P diag(vec(Ω))P ⊗ P are of (near) low ranks, then with µ decreasing, the condition
number κ(M ) may probably increase.
5.2
A diagonal preconditioner
To achieve faster convergence for the PCG method to solve (10), one may select a proper
preconditioner. In our implementation, we devise an easy-to-compute diagonal preconditioner by using an idea ﬁrst developed in [8]. The preconditioner has the following
form
MD := diag(Q) + σdiag(AT A) + µdiag(vec(Ψ)) + σdiag(d),
where Ψ ∈ S n , d ∈ Rn , and
2
(Ψ)ij = (X −1 )ii (X −1 )jj + (X −1 )2ij ,
d(ij) = ((P ◦ P )Ω(P ◦ P )T )ij , 1 ≤ i, j ≤ n.
The biggest advantage of the preconditioner MD is that only O(n3 ) ﬂops are needed to
compute it.
6
Numerical experiments
In this section, we present some numerical results to demonstrate the performances of
the approaches with both synthetic and real data. We implemented the approaches in
Matlab R2013a. All runs were performed on a PC (Intel Core 2 Duo 2.60 GHz with 12
GB RAM).
15
We measure the infeasibilities and optimality for the primal and dual problems by
RP , RD and RG which are described in Algorithms 1 and 2. In general cases, we use
the PAL method to generate an initial point and then switch to the NAL method for
accelerating the local convergence. We stop the PAL method when
max{RP , RD } < Tol1,
and stop the NAL method when
max{RP , RD } < Tol2,
where Tol1 and Tol2 are pre-speciﬁed accuracy tolerances with Tol1 = 5 × 10−3 and
Tol2 = 10−6 as the default. In some hard cases, we use the PAL method only and
terminate it with Tol1 = 10−6 . Furthermore, we also use the relative gap
RG :=
|pobj − dobj|
1 + |pobj| + |dobj|
to measure the quality of the solution, where pobj and dobj denote the primal and dual
objective function values, respectively. For the PAL method, we cap the iteration number
to be 1000; for the Newton-CG method, we set the maximal outer iteration number to
be 100 and the maximal inner iteration number to be 15. We choose the initial iterate
X0 = I, and σ0 = 10. If the NAL method is applied to the above problems, the inequality
constraint X ≽ 0 is replaced by X ≽ εI with ε = 10−16 .
6.1
Synthetic experiments I
In this subsection, we focus our numerical experiments on the following special problems
}
{1
⟨X, H ◦ X⟩ + ⟨C, X⟩ − µ log det X | Xii = 1, i = 1, 2, · · · , n, X ≽ 0
(A1) min
X
2
{
}
Xij = Xji = 0, ∀ (i, j) ∈ E
1
(A2) min
⟨X, H ◦ X⟩ + ⟨C, X⟩ − µ log det X X
2
Xii = 1, i = 1, 2, · · · , n, X ≽ 0
{1
}
(A3) min
⟨X, H ◦ X⟩ + ⟨C, X⟩ − µ log det X | Tr(X) = 1, X ≽ 0
X
2
}
{
Xij = Xji = 0, ∀ (i, j) ∈ E
1
(A4) min
⟨X, H ◦ X⟩ + ⟨C, X⟩ − µ log det X .
X
2
Tr(X) = 1, X ≽ 0
The matrices H and C are generated randomly, and E is a random subset of {(i, j) | 1 ≤
i < j, j = 2, . . . , n}.
The performances of the algorithms are presented in the following tables. For each
instance, we report the matrix dimension (n) and the number of the equality constraints
16
(m); the number of outer iterations (it) and the total number of inner iterations (itsub);
the primal (pobj) and dual (dobj) objective values; the primal infeasibility (RP ), the dual
infeasibility (RD ), the relative gap (RG ); the time (in the format of hours:minutes:seconds)
taken.
Firstly, we present the performances of the pure PAL method and the hybrid method
(i.e., the PAL method and the NAL method) on the problems (A1) and (A2) with µ = 1
in Tables 1 and 2, respectively. But for the hybrid method, we only report the details of
the Newton-CG method. From the two tables, we can see that although the PAL method
can solve the problems rapidly, the hybrid method is still superior to it for all instances.
Especially for instances with relatively few constraints, the hybrid method outperforms
the PAL method obviously. Actually, the superiority of the hybrid method will be even
obvious if we want to get a highly accurate solution.
Table 1: Performances of the PAL method on problems (A1)-(A4)
with µ = 1.
Problem
n|m
it
rand-A1-0500
rand-A1-1000
rand-A1-1500
rand-A1-2000
rand-A2-0500
rand-A2-1000
rand-A2-1500
rand-A2-2000
rand-A3-0500
rand-A3-1000
rand-A3-1500
rand-A3-2000
rand-A4-0500
rand-A4-1000
rand-A4-1500
rand-A4-2000
500 | 500;
1000 | 1000;
1500 | 1500;
2000 | 2000;
500 | 12322;
1000 | 25086;
1500 | 38389;
2000 | 51322;
500 | 1;
1000 | 1;
1500 | 1;
2000 | 1;
500 | 74080;
1000 | 296243;
1500 | 666966;
2000 | 1186173;
93|
88|
88|
95|
106|
102|
99|
95|
570|
1000|
1000|
1000|
1000|
1000|
1000|
1000|
pobj
1.51611644
2.88889067
4.02954208
5.12194289
1.48549608
2.92216880
4.07660888
5.14445107
3.10830391
6.90874980
1.09705221
1.52001337
3.10830458
6.90848660
1.09609128
1.51822498
17
dobj
3
3
3
3
3
3
3
3
3
3
4
4
3
3
4
4
1.51612384
2.88891640
4.02959372
5.12194882
1.48549478
2.92216820
4.07660986
5.14444160
3.10830441
6.90875463
1.09708284
1.52027947
3.10830554
6.90875346
1.09707760
1.52026740
3
3
3
3
3
3
3
3
3
3
4
4
3
3
4
4
RP /RD /RG
Time
4.2-7| 9.4-7| 2.4-6
6.8-7| 9.9-7| 4.5-6
8.9-7| 9.9-7| 6.4-6
9.2-7| 2.5-7| 5.8-7
5.1-7| 9.2-7| 4.4-7
2.7-7| 9.6-7| 1.0-7
6.6-8| 9.9-7| 1.2-7
3.3-7| 9.1-7| 9.2-7
9.8-7| 9.3-7| 7.9-8
4.8-6| 4.6-6| 3.5-7
2.0-4| 1.9-4| 1.4-5
1.3-3| 1.3-3| 8.8-5
1.6-4| 3.5-4| 1.5-7
2.8-4| 4.7-4| 1.9-5
6.6-3| 7.0-3| 4.5-4
1.0-2| 7.5-3| 6.7-4
18
1:20
3:41
8:13
20
1:47
4:28
8:53
1:08
10:36
31:24
1:16:31
3:29
17:55
46:26
1:41:03
Table 2: Performances of the hybrid method on problems (A1)(A4) with µ = 1.
Problem
n|m
it|itsub
rand-A1-0500
rand-A1-1000
rand-A1-1500
rand-A1-2000
rand-A2-0500
rand-A2-1000
rand-A2-1500
rand-A2-2000
rand-A3-0500
rand-A3-1000
rand-A3-1500
rand-A3-2000
rand-A4-0500
rand-A4-1000
rand-A4-1500
rand-A4-2000
500 | 500;
1000 | 1000;
1500 | 1500;
2000 | 2000;
500 | 12322;
1000 | 25086;
1500 | 38389;
2000 | 51322;
500 | 1;
1000 | 1;
1500 | 1;
2000 | 1;
500 | 74080;
1000 | 296243;
1500 | 666966;
2000 | 1186173;
5| 8
6| 9
6| 8
6| 7
10| 11
11| 12
10| 11
10| 11
13| 15
15| 17
20| 22
16| 19
32| 55
36| 63
52| 92
79| 140
pobj
1.51610927
2.88889217
4.02952847
5.12195286
1.48549372
2.92217094
4.07661066
5.14445263
3.10830441
6.90875466
1.09708296
1.52028038
3.10830709
6.90875652
1.09708312
1.52028052
dobj
3
3
3
3
3
3
3
3
3
3
4
4
3
3
4
4
1.51593972
2.88900007
4.02940576
5.12224022
1.48549527
2.92217815
4.07661215
5.14445533
3.10830442
6.90875467
1.09708296
1.52028038
3.10830709
6.90875652
1.09708312
1.52028052
3
3
3
3
3
3
3
3
3
3
4
4
3
3
4
4
RP /RD /RG
Time
5.3-7| 1.7-7| 5.6-5
1.3-7| 3.9-8| 1.9-5
3.2-7| 8.2-8| 1.5-5
5.1-7| 1.4-7| 2.8-5
7.7-7| 6.4-7| 5.2-7
3.8-7| 2.0-7| 1.2-6
3.8-7| 2.3-7| 1.8-7
4.5-7| 2.4-7| 2.6-7
3.6-9| 5.8-7| 2.4-9
1.3-9| 4.1-7| 7.4-10
2.8-9| 9.0-7| 9.6-10
1.2-9| 7.9-7| 6.6-10
1.4-10| 7.6-7| 1.5-12
5.5-11| 7.4-7| 4.4-13
8.7-11| 6.9-7| 3.6-13
1.9-10| 9.1-7| 5.9-13
08
41
1:47
3:44
15
1:20
3:26
7:50
35
5:25
22:55
1:07:14
1:37
13:39
59:34
2:09:36
Secondly, we present the performances of the PAL method and the hybrid method on
problem (A1) with various µ. From Tables 3 and 4, we can see that with µ decreasing,
the computing time increases very mildly. In contrast, with µ decreasing, the computing
time increases very obviously. The rationale is analyzed behind Proposition 5.1, that is,
in this problem, Q, AT A and P ⊗ P diag(vec(Ω))P ⊗ P are of near low ranks and a small
µ leads to the ill-conditioning of the inner problem, so the NAL method needs many PCG
iterations to solve the linear system of equations, however, the PAL method directly gives
a closed-form solution of a regularized problem.
18
Table 3: Performances of the PAL method on problem (A1) with
various µ.
Problem
rand1000-1
rand1000-2
rand1000-3
rand1000-4
rand1000-5
rand1000-6
rand1000-7
rand1000-8
rand1000-9
rand1000-10
n|m
µ
it
|
|
|
|
|
|
|
|
|
|
1.0000;
0.7500;
0.5000;
0.2500;
0.1250;
0.0625;
0.0313;
0.0157;
0.0078;
0.0039;
88|
113|
112|
126|
171|
169|
179|
179|
180|
180|
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
pobj
2.88889067
2.84053903
2.77900646
2.69409906
2.63413549
2.59461152
2.56991656
2.55504884
2.54622819
2.54123460
RP /RD /RG
dobj
3
3
3
3
3
3
3
3
3
3
2.88891640
2.84053619
2.77901061
2.69412261
2.63413210
2.59461261
2.56991919
2.55505345
2.54623730
2.54124523
3
3
3
3
3
3
3
3
3
3
6.8-7|
8.2-7|
8.8-7|
1.6-7|
8.6-7|
8.8-7|
7.7-7|
9.4-7|
7.7-7|
8.2-7|
9.9-7|
1.9-7|
5.5-7|
9.9-7|
2.2-7|
4.7-7|
4.9-7|
6.4-7|
7.2-7|
8.0-7|
Time
4.5-6
5.0-7
7.5-7
4.4-6
6.4-7
2.1-7
5.1-7
9.0-7
1.8-6
2.1-6
1:20
1:42
1:41
1:57
2:45
2:41
2:44
2:47
2:45
2:45
Table 4: Performances of the hybrid method on problem (A1) with
various µ.
Problem
rand1000-1
rand1000-2
rand1000-3
rand1000-4
rand1000-5
rand1000-6
rand1000-7
rand1000-8
rand1000-9
rand1000-10
n|m
µ
it|itsub
|
|
|
|
|
|
|
|
|
|
1.0000;
0.7500;
0.5000;
0.2500;
0.1250;
0.0625;
0.0313;
0.0157;
0.0078;
0.0039;
6| 9
6| 9
6| 11
12| 24
12| 24
11| 22
11| 22
11| 22
11| 26
11| 30
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
pobj
2.88889217
2.84053768
2.77900549
2.69409749
2.63412896
2.59460453
2.56991042
2.55504168
2.54622239
2.54122850
19
RP /RD /RG
dobj
3
3
3
3
3
3
3
3
3
3
2.88900007
2.84064543
2.77914546
2.69409705
2.63412821
2.59460367
2.56990961
2.55504091
2.54622168
2.54122780
3
3
3
3
3
3
3
3
3
3
1.3-7|
1.4-7|
1.9-7|
9.5-8|
1.6-7|
1.9-7|
1.8-7|
1.7-7|
1.6-7|
1.5-7|
3.9-8|
3.6-8|
5.6-8|
4.7-7|
8.3-7|
9.8-7|
9.3-7|
9.1-7|
8.5-7|
8.4-7|
1.9-5
1.9-5
2.5-5
8.2-8
1.4-7
1.7-7
1.6-7
1.5-7
1.4-7
1.4-7
Time
43
45
59
3:18
3:51
5:12
6:59
9:08
12:17
15:50
1000
800
800
time (seconds)
time (seconds)
1000
600
400
200
0 0
10
600
400
200
−1
10
0 0
10
−2
µ
10
−1
10
−2
µ
10
Figure 1: Comparison of the performances of two methods on problem (A1) with various
µ.
6.2
Synthetic experiments II
In this subsection, we focus the numerical experiments on the following problem:
{1
}
(B1) min
∥H ◦ (X − G)∥2 − µ log det X | Xii = 1, i = 1, 2, · · · , n, X ≽ 0
X
2
}
{
Xij = Xji = 0, ∀ (i, j) ∈ E
1
∥H ◦ (X − G)∥2 − µ log det X .
(B2) min
X
2
Xii = 1, i = 1, 2, · · · , n, X ≽ 0
The matrices H and G are generated randomly, where H is a weighted matrix whose
entries are between 0 and 1, and many entries are nearly zero, and E is a random subset
of {(i, j) | 1 ≤ i < j, j = 2, . . . , n}.
The objective function of (B1) or (B2) is
1
1
⟨X, H ◦ H ◦ X⟩ − ⟨H ◦ H ◦ G, X⟩ + ⟨H ◦ G, H ◦ G⟩ − µ log det X,
2
2
(15)
which is a special case of (P ), so the algorithms we developed can be applied to solve (B1)(B2). We list the performances of the PAL method and the hybrid method on problems
(B1) and (B2) in Tables 5 and 6 as below. We can see that for these two problems, the
PAL method is much faster than the hybrid method. The reason lies in that the entries
of the weighted matrix H in (15) are small, then the entries of H ◦ H are even smaller
which leads to a more ill-conditioned Hessian of the inner problem and in this situation
the NAL method has no superiority.
20
Table 5: Performances of the PAL method on problems (B1)-(B2)
with µ = 1.
Problem
n|m
it
rand-B1-0500
rand-B1-1000
rand-B1-1500
rand-B1-2000
rand-B2-0500
rand-B2-1000
rand-B2-1500
rand-B2-2000
500 | 500;
1000 | 1000;
1500 | 1500;
2000 | 2000;
500 | 12322;
1000 | 25086;
1500 | 38389;
2000 | 51322;
147|
114|
106|
107|
175|
143|
135|
129|
pobj
3.90611146
1.76566465
4.38503029
8.41256796
5.30756232
1.91004704
4.53956279
8.56585712
RP /RD /RG
dobj
3
4
4
4
3
4
4
4
3.90613831
1.76570729
4.38510787
8.41269558
5.30756434
1.91005117
4.53957194
8.56586838
3
4
4
4
3
4
4
4
8.9-7|
7.5-7|
2.5-7|
1.6-7|
8.6-7|
4.7-7|
3.4-8|
2.3-8|
3.1-7|
9.7-7|
9.4-7|
9.6-7|
2.1-7|
9.8-7|
9.9-7|
9.6-7|
3.4-6
1.2-5
8.8-6
7.6-6
1.9-7
1.1-6
1.0-6
6.6-7
Time
29
1:49
4:47
9:47
35
2:21
6:05
11:53
Table 6: Performances of the hybrid method on problems (B1)(B2) with µ = 1.
Problem
n|m
rand-B1-0500
rand-B1-1000
rand-B1-1500
rand-B1-2000
rand-B2-0500
rand-B2-1000
rand-B2-1500
rand-B2-2000
500 | 500;
1000 | 1000;
1500 | 1500;
2000 | 2000;
500 | 12322;
1000 | 25086;
1500 | 38389;
2000 | 51322;
it
9|
7|
7|
7|
6|
7|
7|
7|
25
20
19
20
17
19
19
19
pobj
3.90611894
1.76566628
4.38503287
8.41257271
5.30756255
1.91004712
4.53956281
8.56585840
RP /RD /RG
dobj
3
4
4
4
3
4
4
4
3.90609819
1.76566528
4.38502532
8.41255864
5.30757122
1.91004680
4.53956184
8.56585294
3
4
4
4
3
4
4
4
6.9-7|
2.7-7|
3.7-7|
4.8-7|
3.5-7|
1.4-7|
1.7-7|
2.7-7|
7.6-7|
5.2-7|
6.5-7|
7.8-7|
9.8-7|
2.8-7|
3.1-7|
4.4-7|
2.7-6
2.9-7
8.6-7
8.4-7
8.2-7
8.2-8
1.1-7
3.2-7
Time
1:24
7:22
22:01
58:27
1:01
6:51
22:13
51:15
Based on the synthetic experiments in the above two subsections, we conclude that
if the Hessian of the inner problem is good-conditioned, we prefer to adopt the hybrid
method, but if the Hessian is ill-conditioned, the PAL method only is a good choice.
6.3
Real data experiments
In this section, we consider the following model problem which is a special case of problem
(P ) with Q ≡ 0:
{
}
(C) min ⟨C, X⟩ − log det X | Xij = Xji = 0, ∀ (i, j) ∈ E, X ≽ 0 ,
X
21
where the matrix C is real data coming from the two gene expression data sets, and E is
a predetermined index set. The model (C) ﬁnds wide applications in covariance selection.
One gene set is the Rosetta Inpharmatics Compendium of gene expression proﬁles
described by Hughes et al. [12]. The data set contains 253 samples with n = 6136
variables. We aim to estimate the covariance matrix of a Gaussian graphic model whose
conditional independence is unknown. Another gene set is the Iconix microarray data
obtained from 255 drug-treated rat livers; see Natsoulis et al. [23] for details. For both
data sets, although our method can handle problems with larger matrix dimensions, we
test only on a subset of the data. We create 5 subsets by taking 500, 1000, 1500, 2000 and
2500 variables with the highest variances, respectively. And as the variances vary widely,
we normalize the sample covariance matrices to have unit variances on the diagonal.
As the model (C) is a special case of (P ) with Q ≡ 0, so the Hessian of the inner
problem may probably be ill-conditioned and we only apply the PAL method to solve it.
Furthermore, we also compare the performance of the PAL method with that of the PPA
proposed in [40].
Table 7: Performances of the PAL method on problem (C).
Problem
Rosetta-0500
Rosetta-1000
Rosetta-1500
Rosetta-2000
Rosetta-2500
Iconix-0500
Iconix-1000
Iconix-1500
Iconix-2000
Iconix-2500
n|m
it
|
|
|
|
|
|
|
|
|
|
13|
12|
12|
12|
12|
28|
25|
21|
21|
16|
500
1000
1500
2000
2500
500
1000
1500
2000
2500
999;
1999;
2999;
3999;
4999;
999;
1999;
2999;
3999;
4999;
pobj
2.29472323
3.06492978
3.56384267
3.94407754
4.19969175
1.30457645
2.28035658
2.92514554
3.42933318
3.84427349
22
RP /RD /RG
dobj
1
1
1
1
1
2
2
2
2
2
2.29473435
3.06495129
3.56386620
3.94410344
4.19971940
1.30457705
2.28035773
2.92514672
3.42933507
3.84427544
1
1
1
1
1
2
2
2
2
2
4.3-7|
4.7-7|
3.7-7|
3.3-7|
2.9-7|
3.8-7|
1.3-7|
1.6-7|
6.6-8|
7.5-8|
5.1-7|
8.7-7|
6.7-7|
5.9-7|
5.4-7|
7.9-7|
9.0-7|
8.2-7|
9.7-7|
9.1-7|
2.4-6
3.5-6
3.3-6
3.2-6
3.3-6
2.3-7
2.5-7
2.0-7
2.8-7
2.5-7
Time
03
16
46
1:49
3:27
08
28
58
2:14
3:05
Table 8: Comparison of the PAL method and the PPA on problem
(C) with accuracy 10−6 .
Problem
Rosetta-0500
Rosetta-1000
Rosetta-1500
Rosetta-2000
Rosetta-2500
Iconix-0500
Iconix-1000
Iconix-1500
Iconix-2000
Iconix-2500
n|m
500
1000
1500
2000
2500
500
1000
1500
2000
2500
|
|
|
|
|
|
|
|
|
|
it
PAL
PPA
13
12
12
12
12
28
25
21
21
16
2
2
2
2
3
3
3
4
4
4
999
1999
2999
3999
4999
999
1999
2999
3999
4999
Time
PAL PPA
pobj
PAL
2.29472323
3.06492978
3.56384267
3.94407754
4.19969175
1.30457645
2.28035658
2.92514554
3.42933318
3.84427349
PPA
1
1
1
1
1
2
2
2
2
2
2.29473468
3.06495844
3.56387453
3.94411131
4.19979828
1.30457676
2.28035671
2.92514506
3.42933324
3.84427356
1
1
1
1
1
2
2
2
2
2
03
16
46
1:49
3:27
08
28
58
2:14
3:05
08
43
2:01
4:51
13:36
14
1:20
5:59
13:27
25:45
Tables 7 and 8 not only show that the PAL method is very eﬃcient to solve the real
data problems, but also show that the PAL method is about 2 to 7 times faster than the
PPA.
Concluding remarks
In this paper, we designed a PAL method and a hybrid method to solve log-determinant
optimization problems. We established the rigorous convergence results based on the
fundamental theoretical frameworks. Extensive numerical experiments conducted on both
synthetic problems and real data problems demonstrated that our methods are very robust
and eﬃcient. Based on the results of the numerical experiments, we conclude that for
general problems, we prefer to adopt the hybrid method; but for some ill-conditioned
problems, the PAL method may be a better choice.
Acknowledgements
I sincerely appreciate the Institute for Mathematical Sciences, National University of
Singapore for supporting me to visit the institute and attend the workshop “Optimization:
Computation, Theory and Modeling” in 2012 so that I can have a good opportunity to
have fruitful discussions with Professors Defeng Sun and Kim-Chuan Toh on the PAL
method. I also appreciate Dr. Xinyuan Zhao in Beijing University of Technology for
many discussions on this topic.
23
References
[1] F. Alizadeh, J. P. A. Haeberly, and O. L. Overton, Complementarity and nondegeneracy in semideﬁnite programming, Math. Programming, 77 (1997), 111–128.
[2] O. Banerjee, L. El Ghaoui, A. d’Aspremont, Model selection through sparse maximum
likelihood estimation, J. Mach. Learn. Res., 9 (2008), 485–516.
[3] J. F. Bonnans and A. Shapiro, Perturbation Analysis of Optimization Problems,
Springer, New York, 2000.
[4] J. Dahl, L. Vandenberghe, and V. Roychowdhury, Covariance selection for nonchordal graphs via chordal embedding, Optim. Methods Softw., 23 (2008), 501–520.
[5] A. Dempster, Covariance selection, Biometrics, 28 (1972), 157–175.
[6] A. d’Aspremont, O. Banerjee, and L. El Ghaoui, First-order methods for sparse
covariance selection, SIAM J. Matrix Anal. Appl., 30 (2008), 56–66.
[7] M. Fazel, T.-K. Pong, D. Sun, and P. Tseng, Hankel matrix rank minimization with
applications to system identiﬁcation and realization, SIAM J. Matrix Anal. Appl., 34
(2013), 946–977.
[8] Y. Gao and D. Sun, Calibrating least squares semideﬁnite programming with equality
and inequality constraints, SIAM J. Matrix Anal. Appl., 31 (2009), 1432–1457.
[9] R. W. Freund and N. M. Nachtigal, A new Krylov subspace method for symmetric
indeﬁnite linear systems, ORNL/TM-12754, 1994.
[10] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, UK, 1985.
[11] R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge University
Press, Cambridge, UK, 1991.
[12] T. R. Hughes, M. J. Marton, A. R. Jones, C. J. Roberts, R. Stoughton, C. D.
Armour, H. A. Bennett, E. Coﬀey, H. Dai, Y. D. He, M. J. Kidd, A. M. King, M.
R. Meyer, D. Slade, P. Y. Lum, S. B. Stepaniants, D. D. Shoemaker, D. Gachotte,
K. Chakraburtty, J. Simon, M. Bard, and S. H. Friend, Functional discovery via a
compendium of expression proﬁles, Cell, 102 (2000), 109–126.
[13] Z. Hu, J. Cao, and L. J. Hong, Robust simulation of global warming policies using
the DICE model, Manage. Sci., 58 (2012), 1–17.
[14] K. F. Jiang, D. F. Sun, and K.-C. Toh, An inexact accelerated proximal gradient
method for large scale linearly constrained convex SDP, SIAM J. Optim., 22 (2012),
1042–1064.
24
[15] L. Li, and K.-C. Toh,, An inexact interior point method for L1-regularized sparse
covariance selection, Math. Program. Comput., 2 (2010), 291–315.
[16] Z. Lu, Smooth optimization approach for sparse covariance selection, SIAM J. Optim.,
19 (2009), 1807–1827.
[17] Z. Lu, Adaptive ﬁrst-order methods for general sparse inverse covariance selection,
SIAM J. Matrix Anal. Appl., 31 (2010), 2000–2016.
[18] Z. Lu and Y. Zhang, Penalty decomposition methods for L0-norm minimization,
Proceedings of Neural Information Processing Systems (NIPS), 2011: 46–54.
[19] B. Martinet, Regularisation d’in´equations variationelles par approximations successives, Rev. Fran¸caise d’Informat. Recherche Op´erationnelle, (1970), 154–159.
[20] F. Meng, D. Sun, and G. Zhao, Semismoothness of solutions to generalized equations
and the Moreau–Yosida regularization, Math. Programming, 104 (2005), 561–581.
[21] G. J. Minty, On the monotonicity of the gradient of a convex function, Paciﬁc J.
Math., 14 (1964), 243–247.
[22] J. J. Moreau, Proximit´e et dualit´e dans un espace Hilbertien, Bull. Soc. Math. France,
93 (1965), 273–299.
[23] G. Natsoulis, C. I. Pearson, J. Gollub, B. P. Eynon, J. Ferng, R. Nair, R. Idury,
M. D. Lee, M. R. Fielden, R. J. Brennan, A. H. Roter and K. Jarnagin, The liver
pharmacological and xenobiotic gene response repertoire, Mol. Syst. Biol., 175 (2008),
1–12.
[24] P. Olsen, F. Oztoprak, J. Nocedal, and S. Rennie, Newton-like methods
for sparse inverse covariance estimation, available at http://www.optimizationonline.org/DB HTML/2012/06/3506.html.
[25] H. Qi and D. Sun, A quadratically convergent Newton method for computing the
nearest correlation matrix, SIAM J. Matrix Anal. Appl., 28 (2006), 360–385.
[26] R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, 1970.
[27] R. T. Rockafellar, A dual approach to solving nonlinear programming problems by
unconstrained optimization, Math. Programming 5 (1973), 354–373.
[28] R. T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM J.
Control Optim., 14 (1976), 877–898.
[29] R. T. Rockafellar, Augmented Lagrangains and applications of the proximal point
algorithm in convex programming, Math. Oper. Res., 1 (1976), 97–116.
25
[30] K. Scheinberg, S. Ma, and D. Goldfarb, Sparse inverse covariance selection via alternating linearization methods, Twenty-Fourth Annual Conference on Neural Information Processing Systems (NIPS), 2010: 2101–2109.
[31] K. Scheinberg and I. Rish, Learning sparse Gaussian Markov networks using a
greedy coordinate ascent approach, in Machine Learning and Knowledge Discovery
in Databases, Lecture Notes in Comput. Sci. 6323, J. L. Balcazar, F. Bonchi, A.
Gionis, and M. Sebag, eds., Springer-Verlag, Berlin, 2010, 196–212.
[32] D. Sun, The strong second order suﬃcient condition and constraint nondegeneracy
in nonlinear semideﬁnite programming and their implications, Math. Oper. Res., 31
(2006), 761–776.
[33] D. Sun, J. Sun, and L. Zhang, The rate of convergence of the augmented Lagrangian
method for nonlinear semideﬁnite programming, Math. Programming, 114 (2008),
349–391.
[34] K.-C. Toh, Primal-dual path-following algorithms for determinant maximization problems with linear matrix inequalities, Comput. Optim. Appl., 14 (1999), 309–330.
[35] K.-C. Toh, An inexact primal-dual path following algorithm for convex quadratic
SDP, Math. Programming, 112 (2008), 221–254.
[36] R. H. T¨
ut¨
unc¨
u, K.-C. Toh, and M. J. Todd, Solving semideﬁnite-quadratic-linear
programs using SDPT3, Math. Programming, 95 (2003), 189–217.
[37] K.-C. Toh, R. H. T¨
ut¨
unc¨
u, and M. J. Todd, Inexact primal-dual path-following algorithms for a special class of convex quadratic SDP and related problems, Pac. J.
Optim., 3 (2007), 135–164.
[38] N.-K. Tsing, M. K. H. Fan, and E. I. Verriest, On analyticity of functions involving
eigenvalues, Linear Algebra Appl. 207 (1994), 159–180.
[39] B. Varadarajan, D. Povey, and S. M. Chu, Quick fmllr for speaker adaptation in
speech recognition, Acoustics, Speech and Signal Processing, 2008. ICASSP 2008.
IEEE International Conference on.
[40] C. Wang, D. Sun, and K.-C. Toh, Solving log-determinant optimization problems by a
Newton-CG primal proximal point algorithm, SIAM J. Optim., 20 (2010), 2994–3013.
[41] J. Yang, D. Sun, and K.-C. Toh, A proximal point algorithm for log-determinant
optimization with group Lasso regularization, SIAM J. Optim., 23 (2013), 857–893.
[42] K. Yosida, Functional Analysis, Springer Verlag, Berlin, 1964.
26
[43] S. Yang, X. Shen, P. Wonka, Z. Lu, and J. Ye, Fused multiple graphical Lasso,
available at http://people.math.sfu.ca/∼zhaosong/ResearchPapers/FMGL.pdf.
[44] X. Yuan, Alternating direction methods for sparse covariance selection, J. Sci. Comput., 51 (2012), 261–273.
[45] X.-Y. Zhao, A Semismooth Newton-CG Augmented Lagrangian Method for Large
Scale Linear and Convex Quadratic SDPs, PhD thesis, National University of Singapore, 2009.
[46] X.-Y. Zhao, D. Sun, and K.-C. Toh, A Newton-CG augmented Lagrangian method
for semideﬁnite programming, SIAM J. Optim., 20 (2010), 1737–1765.
27