Biconvex Landscape In SDP-Related Learning

En-Liang Hu

^{1}

Bo Wang

^{2}

^{1}

Department of Mathematics, Yunnan Normal University, {ynel.hu@gmail.com}.

^{2}

Xiaomi Inc., {wangbo9@xiaomi.com}.

Abstract

Many machine learning problems can be reduced to learning a low-rank positive semidefinite matrix (denoted as $Z$ ), which encounters semidefinite program (SDP). Existing SDP solvers are often expensive for large-scale learning. To avoid directly solving SDP, some works convert SDP into a nonconvex program by factorizing $Z$ as $X X^{⊤}$ . However, this would bring higher-order nonlinearity, resulting in scarcity of structure in subsequent optimization. In this paper, we propose a novel surrogate for SDP-related learning, in which the structure of subproblem is exploited. More specifically, we surrogate unconstrained SDP by a biconvex problem, through factorizing $Z$ as $X Y^{⊤}$ and using a Courant penalty to penalize the difference of $X$ and $Y$ , in which the resultant subproblems are convex. Furthermore, we provide a theoretical bound for the associated penalty parameter under the assumption that the objective function is Lipschitz-smooth, such that the proposed surrogate will solve the original SDP when the penalty parameter is larger than this bound. Experiments on two SDP-related machine learning applications demonstrate that the proposed algorithm is as accurate as the state-of-the-art, but is faster on large-scale learning.

Introduction

In this paper, we consider the minimization of a smooth convex function over symmetric positive semidefinite (PSD) matrices, that is

min Z \in S_{+}^{n} ~ f (Z),

(1)

where $~ f (\cdot)$ is differentiable convex, and $S_{+}^{n} = {Z | Z = Z^{⊤}, Z ⪰ 0, Z \in R^{n \times n}}$ is the cone of symmetric PSD matrices, in which $⪰$ denotes the Löwner partial order. Many machine learning problems can be reduced as (1). Prominent examples include adaptive clustering [\citeauthoryearRoyer2017], sparse PCA [\citeauthoryearD’aspremont et al.2007], connected subgraph detection [\citeauthoryearAksoylar et al.2017], nonlinear dimensionality reduction [\citeauthoryearWeinberger et al.2004], kernel learning [\citeauthoryearZhuang et al.2011], distance metric learning [\citeauthoryearXing et al.2002], multitask learning [\citeauthoryearObozinski et al.2009], streaming model [\citeauthoryearTropp et al.2017], matrix completion [\citeauthoryearSrebro et al.2005], and inference in graphical model [\citeauthoryearErdogdu et al.2017], etc.

Generally speaking, PSD constraint is the most challenging aspect of solving (1). In order to overcome this difficulty, some works [\citeauthoryearBurer and Monteiro2003, \citeauthoryearJourneé et al.2010, \citeauthoryearLi and Tang2016] proposed to factorize $Z$ as $X X^{⊤}$ , and surrogate (1) by a nonconvex program as

min X \in R^{n \times r} ~ f (X X^{⊤}),

(2)

where $X$ is a real matrix and $r ≪ n$ if it is low-rank. Compared with (1), problem (2) does not have PSD constraint, but the objective becomes a nonconvex function.

To solve problem (2), we refer to the existing results in [\citeauthoryearBurer and Monteiro2003, \citeauthoryearJourneé et al.2010, \citeauthoryearLi and Tang2016], and then from the perspectives therein we can summarize the connections between (1) and (2) as follows:

A local minimizer $¯ X$ of problem (2) provides a stationary point $¯ Z = ¯ X {¯ X}^{T}$ of problem (1) if $¯ X$ is rank deficient. Furthermore, the stationary point $¯ Z$ is thus optimal solution of problem (1).
In the case $r = n$ , any local minimizer $¯ X$ of problem (2) provides the stationary point $¯ Z = ¯ X {¯ X}^{T}$ of problem (1), $¯ Z$ is thus optimal solution of (1).

The above results show that a local minimizer $¯ X$ of (2) will solve (1) if we set $r > r a n k (¯ X)$ (i.e., $¯ X$ is rank deficient, refer to [\citeauthoryearLi and Tang2016] for more details). The recent work [\citeauthoryearGe et al.2017] further claimed that: for (2), all local minimizers are also globally optimal, and no high-order saddle points exist. Nevertheless, the factorized formulation (2) would bring higher-order nonlinearity, and then result in scarcity of structure in subsequent optimization.

Different from (2), in this paper we propose a novel surrogate of (1), in which the subproblem structure is convex. Specifically, we convert SDP (1) into a biconvex problem using Courant penalty [\citeauthoryearNocedal and Wright2006], in which the resultant subproblems are convex. We will show that the proposed surrogate contributes to faster computation to optimize (1) especially when $~ f (\cdot)$ is a quadratic function.

Empirical results on two applications demonstrate that our proposal can bring computational gain, which is substantial without loss of accuracy.

Related Work

The interior point method [\citeauthoryearNesterov and Nemirovski1994] can produce highly accurate solutions for SDP problem. However, it is not scalable because at least $O (n^{3})$ time and $O (n^{2})$ space are needed in each iteration. Moreover, the interior point method is difficult to utilize additional information about the problem structure, such as that the target matrix is low-rank. To exploit low-rank structure, a popular approach for solving SDP is the Frank-Wolfe (FW) algorithm [\citeauthoryearJaggi2011], which is based on sparse greedy approximation. For smooth problems, the complexity of FW is $O (N / ϵ^{1.5})$ arithmetic operations [\citeauthoryearJaggi2011], that is: $O (1 / ϵ)$ iterations are demanded to converge to an $ϵ$ -accurate solution, and in each iteration, $O (N log (n) / \sqrt{ϵ})$ arithmetic operations are needed, where $N$ is the number of nonzero entries in gradient. When the gradient of objective is dense, FW can still be expensive.

Alternatively, as $Z$ in (1) is symmetric and positive semidefinite, it can be factorized as $X X^{⊤}$ . Thus, problem (1) can then be rewritten as a nonconvex program as (2) [\citeauthoryearBurer and Monteiro2003]. It is known that a local minimizer $¯ X$ of (2) corresponds to a minimizer of (1) if $¯ X$ is rank-deficient [\citeauthoryearBurer and Monteiro2003, \citeauthoryearJourneé et al.2010, \citeauthoryearLi and Tang2016]. For linear objective and linear constraints, the nonconvex program (2) has been solved with L-BFGS [\citeauthoryearBurer and Monteiro2003, \citeauthoryearNocedal and Wright2006]. However, the convergence properties are unclear when solving a general nonconvex program by L-BFGS. Block-cyclic coordinate minimization has been used in [\citeauthoryearHu et al.2011] to solve a special nonconvex program of SDP, but a closed-form solution is preferred in each block coordinate update, which might be overly restrictive.

Notation

The transpose of vector or matrix is denoted by the superscript $^{T}$ . The identity matrix is denoted by $I$ with appropriate size. For a matrix $A = [A_{i j}]$ , $tr (A)$ is its trace, $∥ A ∥ = \sqrt{\sum_{i j} A_{i j}^{2}}$ is its Frobenius norm, and $v e c (A)$ unrolls $A$ into a column vector. For two matrices $A, B$ , $⟨ A, B ⟩ = tr (A^{⊤} B)$ , and $A \otimes B$ is their Kronecker product. For a function $g (x)$ , $\nabla_{x} g (\cdot)$ is the derivative w.r.t. variable $x$ , and $g$ is $L$ -smooth and $σ$ -strongly convex if for any $x_{1}, x_{2}$ , there exists a Lipschitz constant $0 \leq σ \leq L < \infty$ such that $σ ∥ x_{1} - x_{2} ∥ \leq ∥ \nabla_{x} g (x_{1}) - \nabla_{x} g (x_{2}) ∥ \leq L ∥ x_{1} - x_{2} ∥$ . If $g (x, y)$ is a function w.r.t. $x$ and $y$ , then $g (\cdot, y)$ and $g (x, \cdot)$ are functions w.r.t. $x$ with constant $y$ and w.r.t. $y$ with constant $x$ respectively.

Problem Formulation

Motivation

For current large-scale optimization in machine learning, the first-order optimization methods, such as (accelerated/stochastic) gradient descent, are popular because of their computational simplicity. In the class of these methods, the most expensive operation usually lies in repeatedly searching a stepsize for objective descent at each iteration, which involves computation of the objective many times. If a good stepsize is possible in an analytical and simple form, the computation complexity would be decreased greatly.

However, a simple stepsize is impossible in (2) even if $~ f$ is a simple quadratic function. This is because that when $~ f$ in (1) is quadratic w.r.t. $Z$ , the objective in (2) rises into a quartic function w.r.t. $X$ , such that the stepsize search in (2) involves some complex computations including constructing the coefficients of a quartic polynomial, solving a cubic equation and comparing the resulted different solutions [\citeauthoryearBurer and Choi2006]. Again, this needs compute the objective many times.

The above difficulties motivate us to propose biconvex formulation (3) instead of (2). In each subproblem of (3), the stepsize search is usually easier than in (2). For example, when $~ f$ is a quadratic function, the subobjective $F (\cdot, Y; γ)$ or $F (X, \cdot; γ)$ in (3) is still quadratic w.r.t. $X$ or $Y$ respectively. Hence, the optimal stepsize to descend objective w.r.t $X$ or $Y$ is analytical, simple and unique for (3) ((16) is an example).

Problem Formulation

Instead of factorizing $Z$ as $X X^{⊤}$ , we factorize $Z$ as $X Y^{⊤}$ , and penalize the difference between $X$ and $Y$ , which can be formulated as

min X, Y \in R^{n \times r} F (X, Y; γ) \equiv f_{s} (X, Y) + \frac{γ}{2} {∥ X - Y ∥}^{2},

(3)

where $γ > 0$ is a penalty parameter and $f_{s} (X, Y)$ is the symmetrization of $~ f (X Y^{⊤})$ in terms of $X$ and $Y$ , namely $f_{s} (X, Y) = f_{s} (Y, X)$ is satisfied. For example, $f_{s}$ can be defined as

(4)

and we will see (in proof of Theorem 1) that the symmetry of $f_{s}$ is key to bound $γ$ . Note that the penalty term in (3) is the classic quadratic (or Courant) penalty for the constraint $X = Y$ . In general, problem (3) approaches problem (2) only when $γ$ goes to infinity. However, we will show (in Theorem 1) that $γ$ can be bounded if $f_{s} (\cdot, Y)$ and $f_{s} (X, \cdot)$ are Lipschitz-smooth, such that when $γ$ is larger than that bound, (3) is equivalent to (2) in that a stationary point $(¯ X, ¯ Y)$ of (3) provides $¯ X$ as a stationary point of (2). Hence, any local minimizer of (3) will produce optimal solution of (1) if $r$ is set as large as enough. Moreover, it will be shown that the objective in (3) is biconvex, and thus we can choose a convex optimization algorithm to decrease the objective function $F$ w.r.t. $X$ and $Y$ alternately. Consequently, this opens a door to surrogate an unconstrained SDP by biconvex optimization.

Contributions

Our main contributions can be summarized as follows:

We propose a biconvex penalty formulation (3) to surrogate unconstrained SDP. This opens a door to solve many SDP-based learning problems using biconvex optimization.
We propose a theoretical bound of the penalty parameter $γ$ . This bridges the stationary point of (3) to the stationary point of (2) when $γ$ is set as larger than this bound.
We show that problem (3) is evidently easier to be solved than (2) if $~ f (\cdot)$ is especially a quadratic function. This is validated on two SDP-related applications in the experimental section.

Proposed Algorithm

Biconvex Penalty

We first give the definition of biconvexity as follows

Definition 1.

$g (x, y)$ is biconvex if for any $x_{i}, y_{i}$ , and any $0 \leq α_{i}, β_{i} \leq 1$ ( $i = 1, 2$ ) such that $α_{1} + α_{2} = β_{1} + β_{2} = 1$ , we have $g (α_{1} x_{1} + α_{1} x_{2}, β_{1} y_{1} + β_{2} y_{2}) \leq α_{1} β_{1} g (x_{1}, y_{1}) + α_{1} β_{2} g (x_{1}, y_{2}) + α_{2} β_{1} g (x_{2}, y_{1}) + α_{2} β_{2} g (x_{2}, y_{2})$ . Namely, $g (x, y)$ is convex in $x$ with fixed $y$ and convex in $y$ with fixed $x$ respectively.

The following Proposition shows that $f_{s}$ in (4) and the objective in (3) are biconvex if $~ f$ in (1) is convex.

Proposition 1.

If $~ f$ in (1) is convex w.r.t. $Z$ , then $f_{s}$ in (4) and $F$ in (3) are biconvex w.r.t. $X$ and $Y$ .

We assume that $f_{s}$ in (3) is coercive and satisfies the following assumption.

Assumption 1.

$f_{s} (\cdot, Y)$ is $L^{X}$ -smooth and $σ^{X}$ -strongly convex, and $f_{s} (X, \cdot)$ is $L^{Y}$ -smooth and $σ^{Y}$ -strongly convex.

The Theoretical Bound

In the following Theorem, we provide a theoretical bound of penalty parameter $γ$ to connect (3) to (2), and thus to (1).

Theorem 1.

If $(¯ X, ¯ Y)$ is a stationary point of (3), then $¯ X = ¯ Y$ when $γ > \frac{1}{4} min {L^{X} - σ^{X}, L^{Y} - σ^{Y}}$ .

Corollary 1.

When the condition of Theorem 1 is satisfied, we have: (i) $¯ X$ is a stationary point of (2) if $(¯ X, ¯ Y)$ is a stationary point of (3); (ii) $¯ X$ is a local minimizer of (2) if $(¯ X, ¯ Y)$ is a local minimizer of (3).

Remark 1.

Following Corollary 1, a local minimizer $(¯ X, ¯ Y)$ of (3) provides $¯ X$ as a local minimizer of (2) if $γ > \frac{1}{4} min {L^{X} - σ^{X}, L^{Y} - σ^{Y}}$ . Furthermore, the produced $¯ Z (= ¯ X {¯ X}^{⊤})$ will solve (1) if the condition of rank deficiency (i.e., $r > r a n k (¯ X)$ ) is satisfied [\citeauthoryearLi and Tang2016].

Our Algorithm

Now we focus on solving (3) under the assumption that the penalty parameter $γ$ satisfies Theorem 1. Depending on the convexity of subproblems, we can solve (3) in alternate convex search. Specifically, we can solve the convex subproblem w.r.t. $X$ with fixed $Y$ and the convex subproblem w.r.t. $Y$ with fixed $X$ alternately, that is

	$X_{k}$	$=$	${a r g m i n}_{X \in R^{n \times r}} F (X, Y_{k - 1}; γ);$		(5)
	$Y_{k}$	$=$	${a r g m i n}_{Y \in R^{n \times r}} F (X_{k}, Y; γ) .$		(6)

In fact, it makes sense to solve (5) and (6) only inexactly. An example is multiconvex optimization [\citeauthoryearXu and Yin2013], in which an alternating proximal gradient descent was proposed to solve multiconvex objective.

We first give a brief introduction to accelerated proximal gradient (APG) descent method, and then fuse it into our algorithm. A well-known APG method is FISTA [\citeauthoryearBeck and Teboulle2009], which was proposed to solve convex composite objective as like as

{min}_{x} ℓ (x) + h (x),

where $ℓ (x)$ is differentiable loss and $h (x)$ is nondifferentiable regularizer. When $h (x) = 0$ (i.e., ’proximal’ is absent, as used in our case later), FISTA is degenerated into computing the sequence ${x_{s}}$ via the iteration

(7)

where $τ_{s}$ is stepsize for descent.

We can apply (7) (as a degenerated FISTA) to (5) and (6). Since $F (\cdot, Y_{k - 1}; γ)$ is convex and differentiable, regarding it as $ℓ (\cdot)$ we can inexactly solve (5) by using (7) as the process of inner iterations. Analogously, we can inexactly solve (6) by using (7) as inner procedure. The smaller the number of inner iterations, the faster the interaction between the updates of $X_{k}$ and $Y_{k}$ each other. Specially, if we set $1$ as the number of inner iterations to inexactly solve (5) and (6) respectively, optimization in them can be realized as Algorithm 1 (except steps $4$ , $7$ and $11$ ), that is our alternating accelerated gradient descent (AAGD) algorithm, where steps $5$ and $6$ inexactly solve (5) and steps $8$ and $9$ inexactly solve (6).

1: Initialization:

X_{0}

Y_{0}

X_{- 1} = X_{0}

Y_{- 1} = Y_{0}

t_{0} = 1

2: for

k = 1, \dots,

t_{k} = \frac{1}{2} (1 + \sqrt{1 + 4 t_{k - 1}^{2}}), ω_{k} = \frac{t_{k - 1} - 1}{t_{k}}

4: estimate the penalty parameter

γ_{\frac{k}{2}}

;

{^X}_{k - 1} = X_{k - 1} + ω_{k} (X_{k - 1} - X_{k - 2})

;

X_{k} = {^X}_{k - 1} - τ_{k}^{X} \nabla_{X} F ({^X}_{k - 1}, Y_{k - 1}; γ_{\frac{k}{2}})

;

7: estimate the penalty parameter

γ_{k}

;

{^Y}_{k - 1} = Y_{k - 1} + ω_{k} (Y_{k - 1} - Y_{k - 2})

;

Y_{k} = {^Y}_{k - 1} - τ_{k}^{Y} \nabla_{Y} F (X_{k}, {^Y}_{k - 1}; γ_{k})

;

10: if a stopping criterion is satisfied then

11: return

¯ Z = Y_{k} Y_{k}^{⊤}

(or

X_{k} X_{k}^{⊤}

12: end if

13: end for

Algorithm 1 Alternating Accelerated Gradient Descent (AAGD).

From another perspective, Algorithm 1 (except steps $4$ and $7$ and $11$ ) can be regarded as a degenerated realization (with a slight difference in rectified $ω_{k}$ ) of alternating proximal gradient method in multiconvex optimization [\citeauthoryearXu and Yin2013]. Hence, the convergence of Algorithm 1 is solid following [\citeauthoryearXu and Yin2013].

Applications

NPKL: nonparametric kernel learning

Given $n$ patterns, let $M$ be the must-link set containing pairs that should belong to the same class, and $C$ be the cannot-link set containing pairs that should not belong to the same class. We use one of formulations in [\citeauthoryearZhuang et al.2011], which learns a kernel matrix using the following SDP problem

min Z \in S_{+}^{n} {~ f}^{(1)} (Z) = \frac{1}{2} \sum {j, i} \in T (Z_{j i} - T_{j i})^{2} + λ tr (Z Δ)

(8)

where $Z = [Z_{j i}] \in R^{n \times n}$ is the target kernel matrix, $Δ$ is the graph Laplacian matrix of the data, $T = M \cup C \cup {j = i}$ , $T = [T_{j i}]$ with $T_{j i} = 1$ if ${j, i} \in M$ or $j = i$ , and 0 if ${j, i} \in C$ , and $λ$ is a tradeoff parameter. The first term of objective in (8) measures the difference between $Z_{j i}$ and $T_{j i}$ , and the second term $tr (Z Δ)$ encourages smoothness on the data manifold by aligning $Z$ with $Δ$ .

Biconvex reformulation

Let $f_{s}^{(1)} (X, Y) = {~ f}^{(1)} (X Y^{⊤})$ , we have

$f_{s}^{(1)} (X, Y) = \frac{1}{2} \sum {j, i} \in T [(tr (X^{⊤} S_{j i} Y))^{2} - 2 T_{j i} tr (X^{⊤} S_{j i} Y)]$		(9)
	$+ λ tr (X^{⊤} Δ Y) + C$
	$= \frac{1}{2} v e c (X)^{⊤}^A (Y) v e c (X) - v e c (X)^{⊤}^b (Y) + C$

where

$^A (Y)$	$=$	$(Y^{⊤} \otimes I) P P^{⊤} (Y^{⊤} \otimes I)^{⊤},$
$^b (Y)$	$=$	$(Y^{⊤} \otimes I) q,$
$P$	$=$	$\sum {j, i} \in T v e c (S_{j i}) v e c (S_{j i})^{⊤},$
$q$	$=$	$\sum {j, i} \in T T_{j i} v e c (S_{j i}) - λ v e c (Δ),$
$S_{j i}$	$=$	$I (:, j) (I (i, :))^{⊤},$

and $I$ is denoted as the identity matrix with appropriate size, $I (:, i)$ is the $i$ -th column of $I$ , $C$ is a constant.

Clearly, $f_{s}^{(1)}$ is symmetric in terms of $X$ and $Y$ with $T = T^{⊤}$ and $Δ = Δ^{⊤}$ . Plugging $f_{s} = f_{s}^{(1)}$ into (3), we have

$F (X, Y; γ) = f_{s}^{(1)} (X, Y) + \frac{γ}{2} {∥ X - Y ∥}^{2}$
	$=$	$\frac{1}{2} v e c (X)^{⊤} A (Y, γ) v e c (X) - v e c (X)^{⊤} b (Y, γ)$
		$- \frac{γ}{2} v e c (Y)^{⊤} v e c (Y) + C$

where

	$A (Y, γ)$	$=$	$γ I +^A (Y),$
	$b (Y, γ)$	$=$	$(Y^{⊤} \otimes I) (q - γ v e c (I)) .$

So, the subproblem about $X$ (corresponding to (5)) is

	$min X \in R^{n \times r}$	$F (X, Y_{k - 1}; γ_{\frac{k}{2}})$
	$\Leftrightarrow min X \in R^{n \times r}$	$\frac{1}{2} v e c (X)^{⊤} A_{k}^{Y} v e c (X) - v e c (X)^{⊤} b_{k}^{Y}$			(10)

where $A_{k}^{Y} = A (Y_{k - 1}, γ_{\frac{γ}{2}})$ , $b_{k}^{Y} = b (Y_{k - 1}, γ_{\frac{γ}{2}})$ .

Since $F$ is symmetric in terms of $X$ and $Y$ , similar to (10) the subproblem about $Y$ (corresponding to (6)) is

	$min Y \in R^{n \times r}$	$F (X_{k}, Y; γ_{k})$
	$\Leftrightarrow min Y \in R^{n \times r}$	$\frac{1}{2} v e c (Y)^{⊤} A_{k}^{X} v e c (Y) - v e c (Y)^{⊤} b_{k}^{X}$			(11)

where $A_{k}^{X} = A (X_{k}, γ_{k})$ , $b_{k}^{X} = b (X_{k}, γ_{k})$ .

Estimating the bound of $γ$

From (9), the Hessian of $f_{s}^{(1)} (X, Y)$ in terms of variable $X$ is

H (Y) =^A (Y) = (Y^{⊤} \otimes I) P P^{⊤} (Y^{⊤} \otimes I)^{⊤} .

So the Lipschitz smooth constant $L^{X}$ (in Assumption 1) can be estimated as the maximal eigenvalue of $H (Y)$ , that is

	$L^{X}$	$=$	${∥ H (Y) ∥}_{2} \leq ({∥ P ∥}_{2} {∥ Y ∥}_{2})^{2} \leq ({∥ P ∥}_{2} ∥ Y ∥)^{2}$		(12)
		$=$	${∥ P ∥}_{2}^{2} tr (Y^{⊤} Y)$		(12)

where ${∥ \cdot ∥}_{2}$ is the spectral norm, and the last inequality uses ${∥ \cdot ∥}_{2} \leq ∥ \cdot ∥$ , and the last equality uses ${∥ Y ∥}^{2} = tr (Y^{⊤} Y)$ . In practice, $L^{X}$ can be estimated dynamically, such that at $k$ -th iteration, $L^{X}$ can be approximated as

L^{X} \approx L_{k}^{X} ≜ {∥ P ∥}_{2}^{2} tr% (Y_{k - 1}^{⊤} Y_{k - 1}) .

(13)

From the symmetry of $f_{s}^{(1)}$ , similar to (13), $L^{Y}$ can be approximated as

L^{Y} \approx L_{k}^{Y} ≜ {∥ P ∥}_{2}^{2} tr% (X_{k}^{⊤} X_{k}) .

(14)

Thus, in Algorithm 1 we can dynamically set the penalty parameter by Theorem 1 as

γ_{\frac{k}{2}} > \frac{1}{4} L_{k}^{X}, a n d γ_{k} > \frac{}{1} 4 L_{k}^{Y} .

(15)

Additionally, instead of (13) and (14), we can obtain Lipschitz smooth constants $L_{k}^{X}$ and $L_{k}^{Y}$ by line search or backtracking.

Computing the optimal stepsizes

The objective functions of (10) and (11) are quadratic, so the optimal stepsize to descend subobjective w.r.t $X$ or $Y$ is analytical, simple and unique, which are provided by the following Proposition.

Proposition 2.

For subproblems (10) and (11), the optimal stepsizes in Algorithm 1 have closed forms as

τ_{k}^{X} = \frac{{∥ ∥ d_{k}^{X} ∥ ∥}^{2}}{{d_{k}^{X}}^{⊤} A_{k}^{Y} d_{k}^{X}}, τ_{k}^{Y} = \frac{{∥ ∥ d_{k}^{Y} ∥ ∥}^{2}}{{d_{k}^{Y}}^{⊤} A_{k}^{X} d_{k}^{Y}}

(16)

where $d_{k}^{X} = v e c (\nabla_{X} F ({^X}_{k - 1}, Y_{k - 1}; γ_{\frac{k}{2}}))$ and $d_{k}^{Y} = v e c (\nabla_{Y} F (X_{k}, {^Y}_{k - 1}; γ_{k}))$ .

CMVU:colored maximum variance unfolding

The colored maximum variance unfolding (CMVU) is a âcoloredâ variants of maximum variance unfolding [\citeauthoryearWeinberger et al.2004], subjected to class labels information. The slack CMVU [\citeauthoryearSong et al.2008] was proposed as a SDP as

min Z \in S_{+}^{n} {~ f}^{(2)} (Z) = \frac{1}{2} \sum {i, j} \in N (Z_{i i} + Z_{j j} - 2 Z_{i j} - d_{i j})^{2} - λ % t r (H Z H L) .

(17)

where $d_{i j}$ denotes the square Euclidean distance between the $i$ -th and $j$ -th objects, $N$ denotes the set of neighbor pairs, whose distances are to be preserved in the unfolding space, $L$ is a kernel matrix of the labels, $H = [H_{i j}] = [δ_{i j} - \frac{1}{n}]$ centers the data and the labels in the unfolding space, and $λ$ controls the tradeoff between dependence maximization and distance preservation. Let $E_{i j} = I (:, i) - I (:, j)$ , (17) can be reformulated as

min Z \in S_{n}^{+} {~ f}^{(2)} (Z) = \frac{1}{2} \sum {i, j} \in N (E_{i j}^{⊤} Z E_{i j} - d_{i j})^{2} - λ % tr (Z H L H) .

Let $f_{s}^{(2)} (X, Y) = {~ f}^{(2)} (X Y^{⊤})$ , it is easy to know that $f_{s}^{(2)}$ satisfies the symmetry. Let $S_{j i} = E_{i j} E_{i j}^{⊤}$ and viewing $d_{i j}$ as $T_{j i}$ in (8), the subproblems of (17) w.r.t. $X$ and $Y$ can be derived and they are similar to (10) and (11) respectively.

Complexity Analysis

The iterations in Algorithms 1 are inexpensive. It mainly costs $O (n r^{2})$ operations for estimating penalty parameter in $(???)$ and $(???)$ since ${∥ P ∥}_{2}$ is a constant needed to be computed only once. As in Proposition 2, it is still $O (n r^{2})$ operations for calculating the optimal stepsizes. The space complexity is scalable because only $O (n r)$ space is needed indeed.

Experiments

In this part, we perform experiments on two SDP-related applications in the previous section: NPKL, CMVU and SD.

The proposed algorithm, AAGD in Algorithm 1, will be compared with the following state-of-the-art baselines in solving SDP in (1).

FW: accelerated Frank-Wolf algorithm [\citeauthoryearLaue2012] to solve (1);
AGD: transform (1) to (2), and solve (2) using the accelerated gradient descent [\citeauthoryearGhadimi and Lan2016, \citeauthoryearLi and Lin2015]. The stepsize search is followed as in [\citeauthoryearBurer and Choi2006].
L-BFGS: transform (1) to (2), and solve (2) using limited memory Broyden-Fletcher-Goldfarb-Shanno method (as similar as done in [\citeauthoryearBurer and Monteiro2003]).
SQLP: semidefinite-quadratic-linear program [\citeauthoryearWu et al.2009], an interior point method, is specified to quadratic SDP. Naturally, SQLP would be faster than a general purpose SDP-solver to solve a quadratic SDP problem like (8) or (17).

All these algorithms are implemented in Matlab. The stopping condition of FW is reached when the relative change in objective value is smaller than $10^{- 6}$ , and then the stopping conditions of the others are reached when their objective values are not larger than the final objective value of FW or the number of iterations exceeds $2000$ . All experiments are run on a PC with a 3.07GHz CPU and 24GB RAM.

The $r$ (rank of the solution) is automatically chosen by the FW algorithm, and is then used by all the other algorithms except SQLP. How to estimate the rank $r$ in (2) and (3) is beyond the scope of this paper. Interested readers can refer to [\citeauthoryearJourneé et al.2010, \citeauthoryearBurer and Monteiro2003] for details.

Results on NPKL

Experiments are performed on the adult data sets¹ (including $9$ sets: a1a $\sim$ a9a, see Table 1) which have been commonly used for benchmarking NPKL algorithms.

As in [\citeauthoryearZhuang et al.2011], the learned kernel matrix is then used for kernelized $k$ -means clustering with the number of clusters $k$ equals to the number of classes. We set $λ = 10$ and repeat $30$ times clustering on each data set with random start point $(X_{0}, Y_{0})$ and random pair constraints ${M$ , $C}$ , and then the average results ( $\pm$ standard deviations) are reported.

	a1a	a2a	a3a	a4a	a5a	a6a	a7a	a8a	a9a
number of samples	1,605	2,265	3,185	4,781	6,414	11,220	16,100	22,696	32,561
number of must-link pairs	1204	2,269	2,389	3,586	4,811	8,415	12,075	17,022	24,421
number of cannot-link pairs	1204	2,269	2,389	3,586	4,811	8,415	12,075	17,022	24,421

Table 1: Adult data sets a1a

\sim

a9a used in the nonparametric kernel learning experiment. Each data sample has 123 features.

	FW	AGD	AAGD	L-BFGS	SQLP
a1a	98.40 $\pm$ 0.48	98.39 $\pm$ 0.48	98.38 $\pm$ 0.48	98.38 $\pm$ 0.53	98.42 $\pm$ 0.50
a2a	98.07 $\pm$ 0.39	98.07 $\pm$ 0.39	98.07 $\pm$ 0.38	98.07 $\pm$ 0.44	98.13 $\pm$ 0.47
a3a	97.70 $\pm$ 0.43	97.69 $\pm$ 0.43	97.69 $\pm$ 0.42	97.72 $\pm$ 0.45	97.71 $\pm$ 0.41
a4a	97.58 $\pm$ 0.31	97.59 $\pm$ 0.29	97.58 $\pm$ 0.30	97.60 $\pm$ 0.23	97.63 $\pm$ 0.36
a5a	97.63 $\pm$ 0.24	97.62 $\pm$ 0.24	97.63 $\pm$ 0.24	97.58 $\pm$ 0.27	–
a6a	97.64 $\pm$ 0.17	97.64 $\pm$ 0.18	97.64 $\pm$ 0.18	97.64 $\pm$ 0.20	–
a7a	97.67 $\pm$ 0.13	97.66 $\pm$ 0.14	97.66 $\pm$ 0.13	–	–
a8a	97.64 $\pm$ 0.11	97.64 $\pm$ 0.12	97.64 $\pm$ 0.11	–	–
a9a	97.55 $\pm$ 0.13	97.55 $\pm$ 0.13	97.55 $\pm$ 0.13	–	–

Table 3: CPU time (in sec) on the adult data sets a1a

\sim

a9a in Table 1. The best and comparable results (according to the pairwise t-test with 95% confidence) are highlighted.

	FW	AGD	AAGD	L-BFGS	SQLP
a1a	59.84 $\pm$ 20.07	52.43 $\pm$ 39.44	37.62 $\pm$ 31.14	131.29 $\pm$ 11.49	339.01 $\pm$ 15.11
a2a	102.94 $\pm$ 28.08	78.32 $\pm$ 44.18	53.15 $\pm$ 30.46	237.75 $\pm$ 16.34	858.34 $\pm$ 26.95
a3a	219.08 $\pm$ 47.97	93.47 $\pm$ 56.90	83.10 $\pm$ 61.77	464.53 $\pm$ 34.75	2049.20 $\pm$ 105.34
a4a	318.17 $\pm$ 75.12	111.64 $\pm$ 44.97	104.29 $\pm$ 41.90	1008.17 $\pm$ 49.09	7342.42 $\pm$ 467.26
a5a	402.82 $\pm$ 94.79	126.36 $\pm$ 73.75	106.57 $\pm$ 41.87	1806.99 $\pm$ 107.91	–
a6a	747.39 $\pm$ 151.68	224.87 $\pm$ 144.58	165.63 $\pm$ 46.01	4961.67 $\pm$ 389.22	–
a7a	1030.07 $\pm$ 226.39	275.92 $\pm$ 88.91	207.83 $\pm$ 60.28	–	–
a8a	1518.96 $\pm$ 270.53	396.67 $\pm$ 195.11	290.96 $\pm$ 67.80	–	–
a9a	2134.26 $\pm$ 544.91	529.23 $\pm$ 131.37	389.74 $\pm$ 91.00	–	–

Table 2: Rand indices (%) on the adult data sets a1a

\sim

a9a in Table 1. The best and comparable results (according to the pairwise t-test with 95% confidence) are highlighted.

As in [\citeauthoryearRand1971], clustering accuracy is measured by the Rand index $\frac{a + b}{0.5 n (n - 1)}$ , where $a$ is the number of pattern pairs belonging to the same class and are placed in the same cluster by $k$ -means, $b$ is the number of pattern pairs belonging to different classes and are placed in different clusters, and the denominator $0.5 n (n - 1)$ is the total number of all different pairs. The higher the Rand index, the better the accuracy.

Results on the clustering accuracy and running CPU time are shown in Tables 3 and 3, respectively. As can be seen, all algorithms obtain almost the same accuracy, and AAGD is the fastest. Figure 1 shows convergence of the objective on the smallest a1a data set and the largest a9a data set. It can be observed that, to arrive the final objective-value (a broken mauve line), AAGD converges significantly much faster than the others. The optimal stepsizes exist as simple and closed forms (i.e., as (16)) for AAGD, while this is unavailable for AGD, which is one of reasons that AAGD is faster than AGD.

Figure 2 shows the progress of ${∥ X_{k} - Y_{k} ∥}^{2}$ with iterations. It converges to zero clearly, implying that the penalty term will vanish after convergence.

Figure 2: The progress of residual value ${∥ X_{k} - Y_{k} ∥}^{2}$ on a9a.

Results on CMVU

Two benchmark data sets², USPS Digits and Newsgroups 20 are used in our experiments, and their information can be referred to [\citeauthoryearSong et al.2008]. As in [\citeauthoryearSong et al.2008], we construct the set $N$ by considering the $1 %$ nearest neighbors of each point. The tradeoff parameter $λ$ is set to $1$ as a default. The problem dimensionality $n$ is straightly set as the number of all data points.

Figure 3 shows the convergences of the objective on the USPS Digits and Newsgroups 20 data sets respectively. It can be observed that, to arrive the final objective-value (a broken mauve line), AAGD converges significantly much faster than the others. On these two data sets, the run time of SQLP exceeds $5000$ seconds, so we don’t plot them.

Conclusion

Many machine learning problems can be reduced to SDP formulation, which would be expensive for large-scale learning. In this paper, we take the scalable learning of SDP problem into consideration. More specifically, we reformulate unconstrained SDP as a biconvex problem with the addition of a Courant penalty, which can be easily optimized using alternating accelerated gradient descent algorithm. Furthermore, we show that when the penalty parameter is larger than our theoretical bound, any local minimizer of the biconvex surrogate will provide a global minimizer of the original SDP problem if the condition of rank deficiency is satisfied. Experiments on two SDP-related problems: nonparametric kernel learning, colored maximum variance unfolding and spectral decomposition, demonstrate that the proposed algorithm is both accurate and scalable.

For the future, we would generalize the biconvex surrogate to more complex circumstances.

Footnotes

Downloaded from http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html.
Downloaded from https://www.cc.gatech.edu/~lsong/code.html..

References

Cem Aksoylar, Lorenzo Orecchia, and Venkatesh Saligrama. Connected subgraph detection with mirror descent on SDPs. In Proceedings of the 34th International Conference on Machine Learning, pages 51–59, 2017.
A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
S. Burer and C. Choi. Computational enhancements in low-rank semidefinite programming. Optimization Methods and Software, 21(3):493–512, 2006.
S. Burer and R.D.C. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming (Series B), 95:329–357, 2003.
A. D’aspremont, El L. Ghaoui, M. I. Jordan, and G. R. G. Lanckriet. A direct formulation for sparse PCA using semidefinite programming. SIAM Review, 49(3):434–48, 2007.
Murat A. Erdogdu, Yash Deshpande, and Andrea Montanari. Inference in graphical models via semidefinite programming hierarchies. In Advances in Neural Information Processing Systems 30, pages 416–424. 2017.
Rong Ge, Chi Jin, and Yi Zheng. No spurious local minima in nonconvex low rank problems: A unified geometric analysis. In Proceedings of the 34th International Conference on Machine Learning, pages 1233–1242, 2017.
Saeed Ghadimi and Guanghui Lan. Accelerated gradient methods for nonconvex nonlinear and stochastic programming. Mathematical Programming, 156(1-2):59–99, 2016.
J. Gorski, F. Pfeuffer, and K. Klamroth. Biconvex sets and optimization with biconvex functions: A survey and extensions. Mathematical Methods of Operations Research, 66(3):373–407, 2007.
En-Liang Hu, Bo Wang, and Song-Can Chen. BCDNPKL: Scalable non-parametric kernel learning using block coordinate descent. In Proceedings of the 28th International Conference on Machine Learning, pages 209–216, Bellevue, WA, USA, June 2011.
M. Jaggi. Convex optimization without projection steps. Preprint arXiv:1108.1170, 2011.
Michel Journeé, F. Bach, Pierre-Antoine Absil, and Rodolphe Sepulchre. Low-rank optimization on the cone of positive semidefinite matrices. SIAM Journal on Optimization, pages 2327–2351, 2010.
Soeren Laue. A hybrid algorithm for convex semidefinite optimization. In John Langford and Joelle Pineau, editors, Proceedings of the 29th International Conference on Machine Learning (ICML-12), ICML ’12, pages 177–184, New York, NY, USA, July 2012. Omnipress.
Huan Li and Zhouchen Lin. Accelerated proximal gradient methods for nonconvex programming. In In Advances in Neural Information Processing Systems (NIPS), pages 379–387, 2015.
Qiuwei Li and Gongguo Tang. The nonconvex geometry of low-rank matrix optimizations with general objective functions. CoRR, abs/1611.03060, 2016.
Xin Liu, Zaiwen Wen, and Yin Zhang. An efficient gauss-newton algorithm for symmetric low-rank product matrix approximations. SIAM Journal on Optimization, 25(3):1571–1608, 2015.
Y. Nesterov and A. Nemirovski. Interior-Point Polynomial Algorithms in Convex Programming. SIAM, 1994.
J. Nocedal and S.J. Wright. Numerical Optimization. Springer, 2006.
G. Obozinski, B. Taskar, and M. I. Jordan. Joint covariate selection and joint subspace selection for multiple classification problems. Statistics and Computing, 20(2):231–252, 2009.
Michael L. Overton and Robert S. Womersley. On the sum of the largest eigenvalues of a symmetric matrix. SIAM Journal on Matrix Analysis and Applications, 13(1):41–45, 1992.
W. M. Rand. Objective criteria for the evaluation of clustering methods. Journal of the American Statistical Association, 66:846–850, 1971.
Martin Royer. Adaptive clustering through semidefinite programming. In Advances in Neural Information Processing Systems 30, pages 1793–1801. 2017.
L. Song, A.J. Smola, K.M. Borgwardt, and A. Gretton. Colored maximum variance unfolding. In Advances in Neural Information Processing Systems 20, 2008.
N. Srebro, J.D.M. Rennie, and T.S. Jaakola. Maximum-margin matrix factorization. In Advances in Neural Information Processing Systems 17, pages 1329–1336, 2005.
Joel A Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher. Fixed-rank approximation of a positive-semidefinite matrix from streaming data. In Advances in Neural Information Processing Systems 30, pages 1225–1234. 2017.
Kilian Q. Weinberger, Fei Sha, and Lawrence K. Saul. Learning a kernel matrix for nonlinear dimensionality reduction. In Proceedings of the Twenty-first International Conference on Machine Learning, pages 839–846, 2004.
Zaiwen Wen, Chao Yang, Xin Liu, and Yin Zhang. Trace-penalty minimization for large-scale eigenspace computation. J. Sci. Comput., 66(3):1175–1203, 2016.
X.-M. Wu, A. M.-C. So, Z. Li, and S.-Y.R. Li. Fast graph Laplacian regularized kernel learning via semidefinite-quadratic-linear programming. In Advances in Neural Information Processing Systems 22, pages 1964–1972, 2009.
E.P. Xing, A.Y. Ng, M.I. Jordan, and S.J. Russell. Distance metric learning with application to clustering with side-information. In Advances in Neural Information Processing Systems, pages 505–512, 2002.
Y. Xu and W. Yin. A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal on Imaging Sciences, 6(3):1758–1789, 2013.
J. Zhuang, I.W. Tsang, and S.C.H. Hoi. A family of simple non-parametric kernel learning algorithms. Journal of Machine Learning Research, 12:1313–1347, 2011.

Comments 0

Request Comment

You are adding the first comment!

How to quickly get a good reply:

Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.

The feedback must be of minimum 40 characters and the title a minimum of 5 characters

Add comment

Cancel

Loading ...

Back to home

403150

Biconvex Landscape In SDP-Related Learning

Abstract

Introduction

Related Work

Notation

Problem Formulation

Motivation

Problem Formulation

Contributions

Proposed Algorithm

Biconvex Penalty

Definition 1.

Proposition 1.

Assumption 1.

The Theoretical Bound

Theorem 1.

Corollary 1.

Remark 1.

Our Algorithm

Applications

NPKL: nonparametric kernel learning

Biconvex reformulation

Estimating the bound of γ

Computing the optimal stepsizes

Proposition 2.

CMVU:colored maximum variance unfolding

Complexity Analysis

Experiments

Results on NPKL

Results on CMVU

Conclusion

Footnotes

References

Request Comment

Estimating the bound of $γ$