# An Efficient Approach to Solve the Large-Scale Semidefinite Programming Problems

Authors: Yongbin Zheng; Yuzhuang Yan; Sheng Liu; Xinsheng Huang; Wanying Xu
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Solving the large-scale problems with semidefinite programming (SDP) constraints is of great importance in modeling and model reduction of complex system, dynamical system, optimal control, computer vision, and machine learning. However, existing SDP solvers are of large complexities and thus unavailable to deal with large-scale problems. In this paper, we solve SDP using matrix generation, which is an extension of the classical column generation. The exponentiated gradient algorithm is also used
to solve the special structure subproblem of matrix generation. The numerical experiments show that our approach is efficient and scales very well with the problem dimension. Furthermore, the proposed algorithm is applied for a clustering problem. The experimental results on real datasets imply that the proposed approach outperforms the traditional interior-point SDP solvers in terms of efficiency and scalability.


---

## Body

## 1. Introduction

Semidefinite programming (SDP) is a technique widely used in modeling of complex systems and some important issues in computer vision and machine learning. Examples of application include model reduction [1], modeling of nonlinear systems [2, 3], optimal control [4], clustering [5–7], robust Euclidean embedding [8], kernel matrix learning [9], and metric learning [10]. It can be seen in [5, 7] that SDP relaxation can produce more accurate estimates than spectral methods. However, the SDP optimization suffers extremely high time complexity. Current SDP solvers utilizing interior-point methods scale as O(n4.5) or worse [1, 6, 10], where n is the number of rows (or columns) of the semidefinite matrix. As a result, SDP can only run on small datasets. It is meaningful to propose a scalable SDP solver.Inspired by column generation [11], which is a classical technology in optimization literature for solving large scale problem, Shen et al. proposed matrix generation in [10] for solving a semidefinite metric learning problem. By using matrix generation, the particular semidefinite program in metric learning was converted into a sequence of linear programs, and it is possible to use the well-developed linear programming technology.In this paper, we propose the matrix generation-based iteration approach to solve general SDP optimization problems. At each iteration, the exponentiated gradient (EG) algorithm [12] is applied to solve the special structure subproblem. Experiments are also carried out on some real datasets to evaluate the proposed method's efficacy and efficiency.The method proposed here can be seen as an extension of column generation to solve SDP problems. The proposed matrix generation method also has the drawback of column generation. It converges slowly of one wants to achieve high accuracy. This is the so-called “tailing effect”. In practice, for many applications, we do not need a very accurate solution. Typically, a moderately accurate solution suffices. On the other hand, the proposed method can also be considered as a generalization of EG to the matrix case. EG is used to solve problems with a vector optimization variable, and the vector must be on a simplex. Here, we generalize EG in the sense that the proposed method solves optimization with a semidefinite matrix whose eigenvalues are on a simplex.We present our main results in the next section.


## 2. The Algorithm

### 2.1. Notation

The following notations will be used throughout the paper:(i)
a bold lower-case letter(x): a column vector,(ii)
an upper-case letter(X): a matrix,(iii)
A≽0: a positive semidefinite (p.s.d.) matrix,(iv)
a≽b: the component-wise inequality between two vectors,(v)
ℝm×n: the vector space of real matrices of size m×n,(vi)
𝕊: the space of real matrices,(vii)
𝕊n: the space of symmetric matrices of size n×n,(viii)
𝕊+n: the space of symmetric positive semidefinite matrices of size n×n,(ix)
〈A,B〉=Tr(A⊤B): the inner product defined on the above spaces,(x)
Tr(·): the trace of a matrix,(xi)
Δn:=x∈ℝn∣x≽0,1⊤x=1: the n-dimensional simplex set for vectors,(xii)
𝒫n:={X∈𝕊+n∣X≽0,Tr(X)=1}: the density matrix set,(xiii)
𝒬n:={X∈𝕊+n∣X≽0,Tr(X)=1,rank(X)=1}: the dyad set of matrices.


### 2.2. Solving Optimization Problem with SDP Constraints Using Matrix Generation

We are interested in the following general convex problem:(2.1)min⁡Xf(X),s.t.X∈Pn.
Here, f(·) is a convex function defined in 𝕊+.We need to extend the Fenchel conjugate [13] to matrices. For a function f(·):ℝn→(-∞,∞], the Fenchel conjugate f*(·):ℝn→(-∞,∞] is the lower semicontinuous and convex function(2.2)f*(u)=sup⁡x∈Rn{u⊤x-f(x)}.
Now, we impose the following two conditions on f(·):  (i)
f(·) is differentiable everywhere,(ii)
the gradientf′(x) is continuous and monotonically increasing at each direction.In other words,f(·) is differentiable and strictly convex. In this case, the Fenchel conjugate is also called Legendre transform and has the following important result:(2.3)u=f′(x).
By analogy, we can define a Fenchel conjugate for a matrix(2.4)f*(U)=sup⁡X∈S{〈U,X〉-f(X)},
as in the vector space ℝn. f*(·) defined in (2.4) must be lower semicontinuous and convex, because f*(·) is a supremum of linear continuous functions of U.We are ready to reformulate our problem (2.1). It is proved in [10] that a p.s.d. matrix can always be decomposed as a linear convex combination of a set of rank-one matrices, we decompose X∈𝕊+ into a convex combination of a set of dyads(2.5)X=∑j=1JθjZj,
with θ=[θ1,…,θJ]⊤∈ΔJ, Zj∈𝒬n,  for all j. So, we can write (2.1) into
(2.6)min⁡θf(X),s.t.X=∑j=1JθjZj,θ≽0,1⊤θ=1,Zj∈Qn,∀j.
Note that here, θ and X are the optimization variables. Although X is redundant in (2.6), we need to keep it to arrive at the important Lagrange dual problem later.This above problems is still very hard to solve, since it has nonconvex rank constraints, and the variableJ is indefinite, because there are an indefinite number of rank-one matrices. However, if we somehow know matrices Zj (j=1,…) a priori, then all the constraints imposed on Zj (j=1,…) can be dropped, and the problem becomes a linear program.We apply matrix generation to solve the problem. Matrix generation is an extension of column generation to nonpolyhedral semidefinite constraints for solving difficult large-scale optimization problems [10]. It is a method to avoid considering all variables of a problem explicitly, and only a small subset of the entire variable set is considered. Once the problem with the small subset variables is solved, the question of if there are any other variables that can be included to improve the solution should be answered. For convex programs, the primal variables correspond to the dual constraints. The dual can be used to solve the above subproblem.The Lagrangian is(2.7)L(θ,X,U,t,p)=f(X)-〈U,X-∑j=1JθjZj〉-p⊤θ-t(1⊤θ-1),
with p≽0, and the dual is(2.8)inf⁡θ,XL=inf⁡X(f(X)-〈U,X〉)+(φ-p-t1︸must be0)⊤θ+t=t-f*(U),
where we have defined φ with φj=〈U,Zj〉, j=1⋯J. The dual problem of (2.6) is(2.9)max⁡t,Ut-f*(U),s.t.φ≽t1.
Essentially the dual is(2.10)max⁡t,Ut-f*(U),s.t.〈U,Zj〉≥t,∀j=1⋯J.We now only consider a small subset of the variables in the primal; that is, only a subset ofZ (denoted by Z¯) is used. The LP solved using Z¯ is usually termed restricted master problem (RMP). Because the primal variables correspond to the dual constraints, solving RMP is equivalent to solve a relaxed version of the dual problem [10]. If all the constraints that we have not added to the dual problem are kept, solving the restricted problem is equivalent to solve the original problem. Otherwise, if there exists at least one constraint that is violated, the violated constraints correspond to variables in primal that are not in RMP. Adding these variables to RMP leads to a new RMP problem that needs to be reoptimized. Here, we have a iterative algorithm that either finds a new Z′ such that(2.11)〈U,Z′〉<t,
where t is the solution of the current restricted problem, or it guarantees that such a Z′ does not exist. To make convergence fast, we find the one by solving the following optimization problem:(2.12)Z′=arg minZ⁡{〈U,Z′〉},s.t.Z∈Qn.
The optimal value Z′ is exactly vv⊤, where v is the eigenvector of U that corresponds to the smallest eigenvalue. For the time being, Z′ is added to the dyads, and the only variable to estimate is θ. The primal can be writhen as(2.13)min⁡θf(θ),s.t.θ∈ΔJ.
We generalize the EG algorithm to the matrix case and use it to solve the above problem.At optimality, because of (2.3), we have the connection between the primal and dual variables(2.14)U⋆=f′(X⋆).
We also know that at least one of the dual constraints takes the equality at optimality, which enable us to obtain t⋆(2.15)t⋆=min⁡j=1⋯J{〈U⋆,Zj〉},t⋆ can be used as the iteration stopping criteria.


### 2.3. The Algorithm Implementation

#### 2.3.1. The Matrix Generation

We propose the implementation of the algorithm for the general convex problem (2.1). To be more general, we extend the SDP constraint X∈𝒫n={X∈𝕊+n∣X≽0,Tr(X)=1} to X∈{X∈𝕊+n∣X≽0,Tr(X)=k}, as shown in(2.16)(2.16)min⁡Xf(X),s.t.X≽0,Tr(X)=k,
where k is a positive integer.The detailed implementation of the algorithm is proposed in Algorithm1. We can observe, at each iteration, that a new matrix is generated, hence the algorithm is named matrix generation. Besides, at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed. This makes the algorithm efficient.Algorithm 1:The matrix generation.
Initialize:(i) The maximum number of integrationsJmax.(ii) The pre-set tolerance valueε (e.g., 10-5).(iii)t0*=0.Iteration 1:(1) randomly selectZ1∈𝒬n; θ=k.(2) then compute:X=θZ1, U*=f′(X), t1*=〈U,Z1〉.While  iterationJ=2,3,…,Jmax  do(1) Ifabs(tJ*-tJ-1*)<ε then break (problem solved);(2) Find a newZJ by finding the eigenvector v corresponding to the smallest eigenvalue ofU* (v is normalized by v⟵(v/∥v∥2) and ZJ=v*v⊤).(3) Find a newθ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}, by solving the followingspectral structure problem by EG algorithm as described in Section2.3.3:minθf(θ),s.t.θ∈RJ,1⊤θ=k,θ≻̲0wheref(θ) is obtained from f(X) by substituting X with ∑j=1JθjZj(4) ComputeX=∑j=1JθjZj, U*=f′(X), tJ*=min{〈U*,ZJ〉}Output: The p.s.d. matrix X=∑j=1JθjZj and minimal value f(X).


#### 2.3.2. The EG Algorithm

In [12], it is shown that EG style updates can converge very quickly compared to other methods. It has been successfully applied to structured prediction problems such as structured support vector machines and conditional random field. We generalize the EG algorithm to matrix to solve the special structure convex optimization problem in MG algorithm. We first give a brief introduction to the EG algorithm [12, 14, 15]. The EG algorithm efficiently solves the convex optimization problem [16](2.17)min⁡θf(θ),s.t.θ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ⪰0},
under the assumption that the object function f(·) is a convex Lipschitz continuous function with Lipschitz constant Lf w.r.t. a fixed given norm ∥·∥. The mathematical definition of Lf is that |f(x)-f(z)|≤Lf∥x-z∥ holds for any x, z in the domain of f(·). The detail of the EG algorithm is shown in Algorithm 2.Algorithm 2:The EG algorithm: solving the special structure problem shown in (2.17).
Initialize:(i)θ0∈ the interior of ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}.(ii) The maximum number of integrationsimaxWhile  iterationi=1,2,…,imax  do(1) Generate the sequence{θi} with:θji=θji-1exp[-τifj′(θi-1)]∑j=1Jθji-1exp[-τifj′(θi-1)]whereτi is step-size, which can be determined by τi=(2logJ/Lf)(1/i) following [15],f′(θ)=[f1′(θ),…,fJ′(θ)]T is the gradient of f(·).(2) Stop if some stopping criteria are met.Output: θi and f(θ).


#### 2.3.3. Evaluation of the Algorithm

In this section, we evaluate the convergence, running speed, and memory consumption of the algorithm by solving(2.18)min⁡X‖K-X‖F2,s.t.X≽0,X1⊤=1,X=X⊤,Tr(X)=k,
where ∥X∥F=∑ijXij2 is the Frobenius norm, K is a n×n symmetric p.s.d. matrix, and k is a positive integer. This problem is important in the issue of the affinity matrix normalization of spectral clustering [17].The above problem is approximate to(2.19)min⁡X(‖K-X‖F2+δ‖X1⊤-1‖22),s.t.X≽0,Tr(X)=k,
where δ is a scalar, for example, δ=103.Before reporting the results, we compare the proposed algorithm with the convex optimization solver, namely, SDPT3 [18] which is used as internal solvers in the disciplined convex programming software CVX [19, 20]. The CVX toolbox can solve standard problems such as linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and semidefinite programs (SDPs).(1) Convergence
We randomly generate a50×50 p.s.d. matrix K in (2.18); let k be 5. In Figure 1, we plot the optimal value obtained by the proposed algorithm at each iteration (the blue curve). The red dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox. We can see that the algorithm converges to the near-optimal solutions quickly.Figure 1
The convergence of the proposed algorithm. The blue curve shows the optimal value obtained by the proposed algorithm at each iteration, and the dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox.(2) Running Time
We solve the problem in (2.16) using the proposed algorithm and the CVX toolbox, respectively, and the running time is shown in Figure 2. The blue curve is the running time of the proposed algorithm versus the matrix dimension, and the red curve is the running time of CVX versus the matrix dimension. We can see that the running time of CVX increases quickly with the matrix dimension and the running time of the proposed algorithm is a fraction of the CVX and scales very well with the matrix dimension. The proposed algorithm has a clear advantage on computational efficiency. One reason is that at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed in the proposed algorithm.Figure 2
Running time of the proposed algorithm. The blue curve shows the running time of the proposed algorithm versus the matrix dimension, and the red curve shows the running time of CVX versus the matrix dimension.(3) Memory Consumption
Figure3 shows the memory consumption of the proposed algorithm and the CXV toolbox. The blue curve is the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve is the approximate memory consumption of CVX versus the matrix dimension. It can be seen that the memory consumption of CVX increases quickly with the matrix dimension n and the memory consumption of the proposed algorithm is a fraction of the CVX and scales very well with the dimension.Figure 3
Memory consumption of the proposed algorithm. The blue curve shows the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve shows the approximate memory consumption of CVX versus the matrix dimension.In summary, the proposed algorithm can converge to the near-optimal solutions quickly and scale very well with the dimension, which are very important for solving large-scale SDP problems. As a result, the proposed method is potential to be applied to computer vision to deal with large matrix of 3D visual data [21] and clustering in data analysis. In next section, we will show how to apply the proposed method to clustering problems.


## 2.1. Notation

The following notations will be used throughout the paper:(i)
a bold lower-case letter(x): a column vector,(ii)
an upper-case letter(X): a matrix,(iii)
A≽0: a positive semidefinite (p.s.d.) matrix,(iv)
a≽b: the component-wise inequality between two vectors,(v)
ℝm×n: the vector space of real matrices of size m×n,(vi)
𝕊: the space of real matrices,(vii)
𝕊n: the space of symmetric matrices of size n×n,(viii)
𝕊+n: the space of symmetric positive semidefinite matrices of size n×n,(ix)
〈A,B〉=Tr(A⊤B): the inner product defined on the above spaces,(x)
Tr(·): the trace of a matrix,(xi)
Δn:=x∈ℝn∣x≽0,1⊤x=1: the n-dimensional simplex set for vectors,(xii)
𝒫n:={X∈𝕊+n∣X≽0,Tr(X)=1}: the density matrix set,(xiii)
𝒬n:={X∈𝕊+n∣X≽0,Tr(X)=1,rank(X)=1}: the dyad set of matrices.


## 2.2. Solving Optimization Problem with SDP Constraints Using Matrix Generation

We are interested in the following general convex problem:(2.1)min⁡Xf(X),s.t.X∈Pn.
Here, f(·) is a convex function defined in 𝕊+.We need to extend the Fenchel conjugate [13] to matrices. For a function f(·):ℝn→(-∞,∞], the Fenchel conjugate f*(·):ℝn→(-∞,∞] is the lower semicontinuous and convex function(2.2)f*(u)=sup⁡x∈Rn{u⊤x-f(x)}.
Now, we impose the following two conditions on f(·):  (i)
f(·) is differentiable everywhere,(ii)
the gradientf′(x) is continuous and monotonically increasing at each direction.In other words,f(·) is differentiable and strictly convex. In this case, the Fenchel conjugate is also called Legendre transform and has the following important result:(2.3)u=f′(x).
By analogy, we can define a Fenchel conjugate for a matrix(2.4)f*(U)=sup⁡X∈S{〈U,X〉-f(X)},
as in the vector space ℝn. f*(·) defined in (2.4) must be lower semicontinuous and convex, because f*(·) is a supremum of linear continuous functions of U.We are ready to reformulate our problem (2.1). It is proved in [10] that a p.s.d. matrix can always be decomposed as a linear convex combination of a set of rank-one matrices, we decompose X∈𝕊+ into a convex combination of a set of dyads(2.5)X=∑j=1JθjZj,
with θ=[θ1,…,θJ]⊤∈ΔJ, Zj∈𝒬n,  for all j. So, we can write (2.1) into
(2.6)min⁡θf(X),s.t.X=∑j=1JθjZj,θ≽0,1⊤θ=1,Zj∈Qn,∀j.
Note that here, θ and X are the optimization variables. Although X is redundant in (2.6), we need to keep it to arrive at the important Lagrange dual problem later.This above problems is still very hard to solve, since it has nonconvex rank constraints, and the variableJ is indefinite, because there are an indefinite number of rank-one matrices. However, if we somehow know matrices Zj (j=1,…) a priori, then all the constraints imposed on Zj (j=1,…) can be dropped, and the problem becomes a linear program.We apply matrix generation to solve the problem. Matrix generation is an extension of column generation to nonpolyhedral semidefinite constraints for solving difficult large-scale optimization problems [10]. It is a method to avoid considering all variables of a problem explicitly, and only a small subset of the entire variable set is considered. Once the problem with the small subset variables is solved, the question of if there are any other variables that can be included to improve the solution should be answered. For convex programs, the primal variables correspond to the dual constraints. The dual can be used to solve the above subproblem.The Lagrangian is(2.7)L(θ,X,U,t,p)=f(X)-〈U,X-∑j=1JθjZj〉-p⊤θ-t(1⊤θ-1),
with p≽0, and the dual is(2.8)inf⁡θ,XL=inf⁡X(f(X)-〈U,X〉)+(φ-p-t1︸must be0)⊤θ+t=t-f*(U),
where we have defined φ with φj=〈U,Zj〉, j=1⋯J. The dual problem of (2.6) is(2.9)max⁡t,Ut-f*(U),s.t.φ≽t1.
Essentially the dual is(2.10)max⁡t,Ut-f*(U),s.t.〈U,Zj〉≥t,∀j=1⋯J.We now only consider a small subset of the variables in the primal; that is, only a subset ofZ (denoted by Z¯) is used. The LP solved using Z¯ is usually termed restricted master problem (RMP). Because the primal variables correspond to the dual constraints, solving RMP is equivalent to solve a relaxed version of the dual problem [10]. If all the constraints that we have not added to the dual problem are kept, solving the restricted problem is equivalent to solve the original problem. Otherwise, if there exists at least one constraint that is violated, the violated constraints correspond to variables in primal that are not in RMP. Adding these variables to RMP leads to a new RMP problem that needs to be reoptimized. Here, we have a iterative algorithm that either finds a new Z′ such that(2.11)〈U,Z′〉<t,
where t is the solution of the current restricted problem, or it guarantees that such a Z′ does not exist. To make convergence fast, we find the one by solving the following optimization problem:(2.12)Z′=arg minZ⁡{〈U,Z′〉},s.t.Z∈Qn.
The optimal value Z′ is exactly vv⊤, where v is the eigenvector of U that corresponds to the smallest eigenvalue. For the time being, Z′ is added to the dyads, and the only variable to estimate is θ. The primal can be writhen as(2.13)min⁡θf(θ),s.t.θ∈ΔJ.
We generalize the EG algorithm to the matrix case and use it to solve the above problem.At optimality, because of (2.3), we have the connection between the primal and dual variables(2.14)U⋆=f′(X⋆).
We also know that at least one of the dual constraints takes the equality at optimality, which enable us to obtain t⋆(2.15)t⋆=min⁡j=1⋯J{〈U⋆,Zj〉},t⋆ can be used as the iteration stopping criteria.


## 2.3. The Algorithm Implementation

### 2.3.1. The Matrix Generation

We propose the implementation of the algorithm for the general convex problem (2.1). To be more general, we extend the SDP constraint X∈𝒫n={X∈𝕊+n∣X≽0,Tr(X)=1} to X∈{X∈𝕊+n∣X≽0,Tr(X)=k}, as shown in(2.16)(2.16)min⁡Xf(X),s.t.X≽0,Tr(X)=k,
where k is a positive integer.The detailed implementation of the algorithm is proposed in Algorithm1. We can observe, at each iteration, that a new matrix is generated, hence the algorithm is named matrix generation. Besides, at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed. This makes the algorithm efficient.Algorithm 1:The matrix generation.
Initialize:(i) The maximum number of integrationsJmax.(ii) The pre-set tolerance valueε (e.g., 10-5).(iii)t0*=0.Iteration 1:(1) randomly selectZ1∈𝒬n; θ=k.(2) then compute:X=θZ1, U*=f′(X), t1*=〈U,Z1〉.While  iterationJ=2,3,…,Jmax  do(1) Ifabs(tJ*-tJ-1*)<ε then break (problem solved);(2) Find a newZJ by finding the eigenvector v corresponding to the smallest eigenvalue ofU* (v is normalized by v⟵(v/∥v∥2) and ZJ=v*v⊤).(3) Find a newθ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}, by solving the followingspectral structure problem by EG algorithm as described in Section2.3.3:minθf(θ),s.t.θ∈RJ,1⊤θ=k,θ≻̲0wheref(θ) is obtained from f(X) by substituting X with ∑j=1JθjZj(4) ComputeX=∑j=1JθjZj, U*=f′(X), tJ*=min{〈U*,ZJ〉}Output: The p.s.d. matrix X=∑j=1JθjZj and minimal value f(X).


### 2.3.2. The EG Algorithm

In [12], it is shown that EG style updates can converge very quickly compared to other methods. It has been successfully applied to structured prediction problems such as structured support vector machines and conditional random field. We generalize the EG algorithm to matrix to solve the special structure convex optimization problem in MG algorithm. We first give a brief introduction to the EG algorithm [12, 14, 15]. The EG algorithm efficiently solves the convex optimization problem [16](2.17)min⁡θf(θ),s.t.θ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ⪰0},
under the assumption that the object function f(·) is a convex Lipschitz continuous function with Lipschitz constant Lf w.r.t. a fixed given norm ∥·∥. The mathematical definition of Lf is that |f(x)-f(z)|≤Lf∥x-z∥ holds for any x, z in the domain of f(·). The detail of the EG algorithm is shown in Algorithm 2.Algorithm 2:The EG algorithm: solving the special structure problem shown in (2.17).
Initialize:(i)θ0∈ the interior of ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}.(ii) The maximum number of integrationsimaxWhile  iterationi=1,2,…,imax  do(1) Generate the sequence{θi} with:θji=θji-1exp[-τifj′(θi-1)]∑j=1Jθji-1exp[-τifj′(θi-1)]whereτi is step-size, which can be determined by τi=(2logJ/Lf)(1/i) following [15],f′(θ)=[f1′(θ),…,fJ′(θ)]T is the gradient of f(·).(2) Stop if some stopping criteria are met.Output: θi and f(θ).


### 2.3.3. Evaluation of the Algorithm

In this section, we evaluate the convergence, running speed, and memory consumption of the algorithm by solving(2.18)min⁡X‖K-X‖F2,s.t.X≽0,X1⊤=1,X=X⊤,Tr(X)=k,
where ∥X∥F=∑ijXij2 is the Frobenius norm, K is a n×n symmetric p.s.d. matrix, and k is a positive integer. This problem is important in the issue of the affinity matrix normalization of spectral clustering [17].The above problem is approximate to(2.19)min⁡X(‖K-X‖F2+δ‖X1⊤-1‖22),s.t.X≽0,Tr(X)=k,
where δ is a scalar, for example, δ=103.Before reporting the results, we compare the proposed algorithm with the convex optimization solver, namely, SDPT3 [18] which is used as internal solvers in the disciplined convex programming software CVX [19, 20]. The CVX toolbox can solve standard problems such as linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and semidefinite programs (SDPs).(1) Convergence
We randomly generate a50×50 p.s.d. matrix K in (2.18); let k be 5. In Figure 1, we plot the optimal value obtained by the proposed algorithm at each iteration (the blue curve). The red dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox. We can see that the algorithm converges to the near-optimal solutions quickly.Figure 1
The convergence of the proposed algorithm. The blue curve shows the optimal value obtained by the proposed algorithm at each iteration, and the dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox.(2) Running Time
We solve the problem in (2.16) using the proposed algorithm and the CVX toolbox, respectively, and the running time is shown in Figure 2. The blue curve is the running time of the proposed algorithm versus the matrix dimension, and the red curve is the running time of CVX versus the matrix dimension. We can see that the running time of CVX increases quickly with the matrix dimension and the running time of the proposed algorithm is a fraction of the CVX and scales very well with the matrix dimension. The proposed algorithm has a clear advantage on computational efficiency. One reason is that at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed in the proposed algorithm.Figure 2
Running time of the proposed algorithm. The blue curve shows the running time of the proposed algorithm versus the matrix dimension, and the red curve shows the running time of CVX versus the matrix dimension.(3) Memory Consumption
Figure3 shows the memory consumption of the proposed algorithm and the CXV toolbox. The blue curve is the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve is the approximate memory consumption of CVX versus the matrix dimension. It can be seen that the memory consumption of CVX increases quickly with the matrix dimension n and the memory consumption of the proposed algorithm is a fraction of the CVX and scales very well with the dimension.Figure 3
Memory consumption of the proposed algorithm. The blue curve shows the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve shows the approximate memory consumption of CVX versus the matrix dimension.In summary, the proposed algorithm can converge to the near-optimal solutions quickly and scale very well with the dimension, which are very important for solving large-scale SDP problems. As a result, the proposed method is potential to be applied to computer vision to deal with large matrix of 3D visual data [21] and clustering in data analysis. In next section, we will show how to apply the proposed method to clustering problems.


## 2.3.1. The Matrix Generation

We propose the implementation of the algorithm for the general convex problem (2.1). To be more general, we extend the SDP constraint X∈𝒫n={X∈𝕊+n∣X≽0,Tr(X)=1} to X∈{X∈𝕊+n∣X≽0,Tr(X)=k}, as shown in(2.16)(2.16)min⁡Xf(X),s.t.X≽0,Tr(X)=k,
where k is a positive integer.The detailed implementation of the algorithm is proposed in Algorithm1. We can observe, at each iteration, that a new matrix is generated, hence the algorithm is named matrix generation. Besides, at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed. This makes the algorithm efficient.Algorithm 1:The matrix generation.
Initialize:(i) The maximum number of integrationsJmax.(ii) The pre-set tolerance valueε (e.g., 10-5).(iii)t0*=0.Iteration 1:(1) randomly selectZ1∈𝒬n; θ=k.(2) then compute:X=θZ1, U*=f′(X), t1*=〈U,Z1〉.While  iterationJ=2,3,…,Jmax  do(1) Ifabs(tJ*-tJ-1*)<ε then break (problem solved);(2) Find a newZJ by finding the eigenvector v corresponding to the smallest eigenvalue ofU* (v is normalized by v⟵(v/∥v∥2) and ZJ=v*v⊤).(3) Find a newθ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}, by solving the followingspectral structure problem by EG algorithm as described in Section2.3.3:minθf(θ),s.t.θ∈RJ,1⊤θ=k,θ≻̲0wheref(θ) is obtained from f(X) by substituting X with ∑j=1JθjZj(4) ComputeX=∑j=1JθjZj, U*=f′(X), tJ*=min{〈U*,ZJ〉}Output: The p.s.d. matrix X=∑j=1JθjZj and minimal value f(X).


## 2.3.2. The EG Algorithm

In [12], it is shown that EG style updates can converge very quickly compared to other methods. It has been successfully applied to structured prediction problems such as structured support vector machines and conditional random field. We generalize the EG algorithm to matrix to solve the special structure convex optimization problem in MG algorithm. We first give a brief introduction to the EG algorithm [12, 14, 15]. The EG algorithm efficiently solves the convex optimization problem [16](2.17)min⁡θf(θ),s.t.θ=[θ1,…,θJ]⊤∈ΔJ={θ∈RJ:1⊤θ=k,θ⪰0},
under the assumption that the object function f(·) is a convex Lipschitz continuous function with Lipschitz constant Lf w.r.t. a fixed given norm ∥·∥. The mathematical definition of Lf is that |f(x)-f(z)|≤Lf∥x-z∥ holds for any x, z in the domain of f(·). The detail of the EG algorithm is shown in Algorithm 2.Algorithm 2:The EG algorithm: solving the special structure problem shown in (2.17).
Initialize:(i)θ0∈ the interior of ΔJ={θ∈RJ:1⊤θ=k,θ≻̲0}.(ii) The maximum number of integrationsimaxWhile  iterationi=1,2,…,imax  do(1) Generate the sequence{θi} with:θji=θji-1exp[-τifj′(θi-1)]∑j=1Jθji-1exp[-τifj′(θi-1)]whereτi is step-size, which can be determined by τi=(2logJ/Lf)(1/i) following [15],f′(θ)=[f1′(θ),…,fJ′(θ)]T is the gradient of f(·).(2) Stop if some stopping criteria are met.Output: θi and f(θ).


## 2.3.3. Evaluation of the Algorithm

In this section, we evaluate the convergence, running speed, and memory consumption of the algorithm by solving(2.18)min⁡X‖K-X‖F2,s.t.X≽0,X1⊤=1,X=X⊤,Tr(X)=k,
where ∥X∥F=∑ijXij2 is the Frobenius norm, K is a n×n symmetric p.s.d. matrix, and k is a positive integer. This problem is important in the issue of the affinity matrix normalization of spectral clustering [17].The above problem is approximate to(2.19)min⁡X(‖K-X‖F2+δ‖X1⊤-1‖22),s.t.X≽0,Tr(X)=k,
where δ is a scalar, for example, δ=103.Before reporting the results, we compare the proposed algorithm with the convex optimization solver, namely, SDPT3 [18] which is used as internal solvers in the disciplined convex programming software CVX [19, 20]. The CVX toolbox can solve standard problems such as linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and semidefinite programs (SDPs).(1) Convergence
We randomly generate a50×50 p.s.d. matrix K in (2.18); let k be 5. In Figure 1, we plot the optimal value obtained by the proposed algorithm at each iteration (the blue curve). The red dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox. We can see that the algorithm converges to the near-optimal solutions quickly.Figure 1
The convergence of the proposed algorithm. The blue curve shows the optimal value obtained by the proposed algorithm at each iteration, and the dash line shows the ground truth obtained by directly solving the original problem in (2.16) using the CVX toolbox.(2) Running Time
We solve the problem in (2.16) using the proposed algorithm and the CVX toolbox, respectively, and the running time is shown in Figure 2. The blue curve is the running time of the proposed algorithm versus the matrix dimension, and the red curve is the running time of CVX versus the matrix dimension. We can see that the running time of CVX increases quickly with the matrix dimension and the running time of the proposed algorithm is a fraction of the CVX and scales very well with the matrix dimension. The proposed algorithm has a clear advantage on computational efficiency. One reason is that at each iteration, only the most significant eigenvalue and its corresponding eigenvector are needed in the proposed algorithm.Figure 2
Running time of the proposed algorithm. The blue curve shows the running time of the proposed algorithm versus the matrix dimension, and the red curve shows the running time of CVX versus the matrix dimension.(3) Memory Consumption
Figure3 shows the memory consumption of the proposed algorithm and the CXV toolbox. The blue curve is the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve is the approximate memory consumption of CVX versus the matrix dimension. It can be seen that the memory consumption of CVX increases quickly with the matrix dimension n and the memory consumption of the proposed algorithm is a fraction of the CVX and scales very well with the dimension.Figure 3
Memory consumption of the proposed algorithm. The blue curve shows the approximate memory consumption of the proposed algorithm versus the matrix dimension, and the red curve shows the approximate memory consumption of CVX versus the matrix dimension.In summary, the proposed algorithm can converge to the near-optimal solutions quickly and scale very well with the dimension, which are very important for solving large-scale SDP problems. As a result, the proposed method is potential to be applied to computer vision to deal with large matrix of 3D visual data [21] and clustering in data analysis. In next section, we will show how to apply the proposed method to clustering problems.


## 3. The Application

The problem of partitioning data points into a number of distinct sets, known as the clustering problem, is one of the most fundamental and significant topics in data analysis and machine learning. From an optimization viewpoint, clustering aims at the minimal sum-of-squares (MSSC) based on some definition of similarly or distance. But the variable is a 0-1 discrete indicator matrixY∈ℝm×n (n is the number of samples, and m is the number of clusters), thus the optimization becomes a mixed integer programming with nonlinear objective, which is NP-hard. The computation time of NP-hard problems is too heavy to use in practical engineering applications and should be accelerated by some approximation methods. K-means and spectral clustering are the two most well-known algorithms to approximately solve the clustering problem. However, the former is merely able to achieve local optimum through an indirect way, and the later, which can obtain a global optimum through eigenvalue decomposition, is unfortunately a weak relaxation of the N-cut problem [7]. Recent research of semidefinite programming (SDP) relaxation on clustering proposed that the minimal sum-of-squares (MSSC) [22] is possible to be relaxed as SDP problems, which is more direct than K-means and tighter than spectral clustering [7].In this section, we apply the proposed algorithm to solve the SDP relaxed clustering problem (0-1 SDP for clustering) [22](3.1)min⁡XTr(W(I-X)),s.t.X1=1,Tr(X)=m,X≥0,X⪰0,
where W is the so-called affinity matrix whose entries represent the similarities or closeness among the entities in the dataset, and m is the class number of the dataset. According to [22], it is reasonable to relax the nonnegative constraint in  (3.1)(3.2)min⁡XTr(W(I-X)),s.t.X1=1,Tr(X)=m,X⪰0.We solve the problem in (3.2) using the algorithm described in Algorithm 1. To test the algorithm, we carry out some experiments on four real datesets collected from UCI machine learning repository  (http://www.ics.uci.edu/~mlearn/MLRepository.html.) The datasets are listed in Table 1 together with some their characteristics. k, Size, and Dim are the class number, the number of instances, and the number of attributes of each dataset, respectively. All the experiments are done by using Matlab on a personal computer with a 2.93 GHz processor and a 2 G memory. The accuracy of clustering each dataset of our algorithm is compared with K-means, spectral clustering, and CVX (here CVX means the SDP problem in (3.2) is solved using the CVX toolbox). The results are shown in Table 1. We can see that(i)
the proposed algorithm performs better than K-means,(ii)
the proposed algorithm achieves higher accuracy than the spectral clustering that is because the proposed algorithm is the SDP relaxation of the clustring problem and is tighter than the spectral relaxation,(iii)
the proposed algorithm is as good as CVX. This means it can solve the SDP problem as accurately as the CVX toolbox.Table 1
UCI datasets used together with some characteristics and the clustering results achieved using the different methods.
DatasetkSizeDims.AccuracyK-meansSpectral clusteringSDP relaxationOursCVXSoybean447360.7230.7450.7660.766Iris315050.8930.9330.9470.953Wine3178130.7020.7250.7300.730SPECTF heart2267440.6250.8020.8160.816We also compare the speed of the proposed algorithm with CVX, and the results are shown in Table2. It shows that the proposed algorithm is much faster than CVX.Table 2
Time consumptions (seconds) of the proposed algorithm and CVX on clustering the selected UCI datasets.
SoybeanIrisWineSPECTF heartOurs1.342.8794.2221.29CVX1.214.6296.9528.00As observed from the experiments, it can be found that the proposed algorithm performs very well on these datasets and can efficiently solve the 0-1 SDP problem for clustering. Both speed and accuracy are improved as compared with traditional interior-point SDP solvers. We plan to further extend the approach to more general problems and look for other applications such as modeling [23, 24] and stabilizing some kinds of dynamical systems [25] model reduction.


## 4. Conclusion

This paper proposed an approach to solve the optimization problem with SDP constraints using matrix generation and exponentiated gradient. The process is faster and more scalable than the traditional SDP approach. The results of numerical experiments show that the method scales very well with the problem dimensions. We also applied the approach to solve the SDP relaxation clustering problem. Experimental results on real datasets show that the algorithm outperforms much in both speed and accuracy as compared with interior-point SDP solvers.



---

*Source: 764760-2011-08-03.xml*