{'fraction_non_alphanumeric': 0.05569956244945358, 'fraction_numerical': 0.016133380337183503, 'mean_word_length': 4.283663854497644, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large-scale Binary Quadratic Optimization Using Semidefinite Relaxation and Applications

In computer vision, many problems such as image segmentation, pixel labelling, and scene parsing can be formulated as binary quadratic programs (BQPs). For submodular problems, cuts based methods can be employed to efficiently solve large-scale problems. However, general nonsubmodular problems are significantly more challenging to solve. Finding a solution when the problem is of large size to be of practical interest, however, typically requires relaxation. Two standard relaxation methods are widely used for solving general BQPs--spectral methods and semidefinite programming (SDP), each with their own advantages and disadvantages. Spectral relaxation is simple and easy to implement, but its bound is loose. Semidefinite relaxation has a tighter bound, but its computational complexity is high, especially for large scale problems. In this work, we present a new SDP formulation for BQPs, with two desirable properties. First, it has a similar relaxation bound to conventional SDP formulations. Second, compared with conventional SDP methods, the new SDP formulation leads to a significantly more efficient and scalable dual optimization approach, which has the same degree of complexity as spectral methods. We then propose two solvers, namely, quasi-Newton and smoothing Newton methods, for the dual problem. Both of them are significantly more efficiently than standard interior-point methods. In practice, the smoothing Newton solver is faster than the quasi-Newton solver for dense or medium-sized problems, while the quasi-Newton solver is preferable for large sparse/structured problems. Our experiments on a few computer vision applications including clustering, image segmentation, co-segmentation and registration show the potential of our SDP formulation for solving large-scale BQPs.


INTRODUCTION
Binary quadratic programs (BQPs) are a class of combinatorial optimization problems with binary variables, quadratic objective function and linear/quadratic constraints. They appear in a wide variety of applications in computer vision, such as image segmentation/pixel labelling, image registration/matching, image denoising/restoration. Moreover, Maximum a Posteriori (MAP) inference problems for Markov Random Fields (MRFs) can be formulated as BQPs too. There are a long list of references to applications formulated as BQPs or specifically MRF-MAP problems. Readers may refer to [1], [2], [3], [4], [5], [6] and the references therein for detailed studies. Unconstrained BQPs with submodular pairwise terms can be solved exactly and efficiently using graph cuts [7], [8], [9]. However solving general BQP problems is known to be NPhard (see [10] for exceptions). In other words, it is unlikely to find polynomial time algorithms to exactly solve these problems. Alternatively, relaxation approaches can be used to produce a feasible solution close to the global optimum in polynomial time. In order to accept such a relaxation we require a guarantee that the divergence between the solutions to the original problem and the relaxed problem is bounded. The quality of the relaxation thus depends upon the tightness of the bounds. Developing an efficient relaxation algorithm with a tight relaxation bound that can achieve a good solution (particularly for large problems) is thus of great practical importance. There are a number of relaxation methods for BQPs (in particular MRF-MAP inference problems) in the literature, including linear programming (LP) relaxation [11], [12], [13], [14], quadratic programming relaxation [15], second order cone relaxation [16], [17], [18], spectral relaxation [19], [20], [21], [22] and SDP relaxation [17], [23].
Spectral methods are effective for many computer vision applications, such as image segmentation [19], [20] and motion segmentation [24]. The optimization of spectral methods eventually lead to the computation of top eigenvectors. Nevertheless, spectral methods may produce loose relaxation bounds in many cases [25], [26], [27]. Moreover, the inherent quadratic programming formulation of spectral methods is difficult to incorporate certain types of additional constraints [21]. SDP relaxation has been shown that it leads to tighter approximation than other relaxation methods for many combinatorial optimization problems [28], [29], [30], [31]. In particular for the max-cut problem, Goemans and Williamson [32] achieve the stateof-the-art 0.879 approximation ratio using SDP relaxation. SDP relaxation has also been used in a range of vision problems, such as image segmentation [33], restoration [34], [35], graph matching [36], [37] and co-segmentation [38]. In a standard SDP problem, a linear function of a symmetric matrix X is optimized, subject to linear (in)equality constraints and the constraint of X being positive semidefinite (p.s.d.). The standard SDP problem and its Lagrangian dual problem are written as: (SDP-P) min where m = |I eq | + |I in |, and I eq (I in ) denotes the indexes of linear (in)equality constraints. The p.s.d. constraint X ∈ S n + is convex, so SDP problems are convex optimization problems and the above two formulations are equivalent if a feasible solution exists. The SDP problem (1) can be considered as a semi-infinite LP problem, as the p.s.d. constraint can be converted to an infinite number of linear constraints: X, aa ≥ 0, ∀a ∈ R n . Through SDP, these infinite number of linear constraints can be handled in finite time.
It is widely accepted that interior-point methods [39], [40] are very robust and accurate for general SDP problems up to a moderate size (see SeDuMi [41], SDPT3 [42] and MOSEK [43] for implementations). However, its high computational complexity and memory requirement hampers the application of SDP methods to large-scale problems. Approximate nonlinear programming methods [44], [45], [46] are proposed for SDP problems based on low-rank factorization, which may converge to a local optimum. Augmented Lagrangian methods [47], [48] and the variants [49], [50] have also been developed. As gradient-descend based methods [51], they may converge slowly. The spectral bundle method [52] and the log-barrier algorithm [53] can be used for large-scale problems as well. A drawback is that they can fail to solve some SDP problems to satisfactory accuracies [48].
In this work, we propose a regularized SDP relaxation approach to BQPs. Preliminary results of this paper appeared in [54]. Our main contributions are as follows.
1) Instead of directly solving the standard SDP relaxation to BQPs, we propose a quadratically regularized version of the original SDP formulation, which can be solved efficiently and achieve a solution quality comparable to the standard SDP relaxation. 2) We proffer two algorithms to solve the dual problem, based on quasi-Newton (referred to as SDCut-QN) and smoothing Newton (referred to as SDCut-SN) methods respectively. The sparse or low-rank structure of specific problems are also exploited to speed up the computation. The proposed solvers require much lower computational cost and storage memory than standard interior-point methods. In particular, SDCut-QN has a lower computational cost in each iteration while needs more iterations to converge. On the other hand, SDCut-SN converges quadratically with higher computational complexity per iteration. In our experiments, SDCut-SN is faster for dense or medium-sized problems, and SDCut-QN is more efficient for large-scale sparse/structured problems. 3) We demonstrate the efficiency and flexibility of our proposed algorithms by applying them to a variety of computer vision tasks. We show that due to the capability of accommodating various constraints, our methods can encode problem-dependent information. More specifically, the formulation of SDCut allows multiple additional linear and quadratic constraints, which enables a broader set of applications than what spectral methods and graph-cut methods can be applied to.
Notation A matrix (column vector) is denoted by a bold capital (lower-case) letter. R n denotes the space of real-valued n × 1 vectors. R n + and R n − represent the non-negative and non-positive orthants of R n respectively. S n denotes the space of n × n symmetric matrices, and S n + represents the corresponding cone of positive semidefinite (p.s.d.) matrices. For two vectors, x ≤ y indicates the element-wise inequality; trace(X), rank(X) and diag(X) denote the trace, rank and the main diagonal elements of X respectively. Diag(x) denotes a diagonal matrix with the elements of vector x on the main diagonal. X 2 F denotes the Frobenius norm of X. The inner product of two matrices is defined as X, Y . I n indicates the n × n identity matrix. 0 and 1 denote all-zero and all-one column vectors respectively.∇f(·) and ∇ 2 f(·) stand for the first-order and second-order derivatives of function f(·) respectively.

BQPS AND THEIR SDP RELAXATION
Let us consider a binary quadratic program of the following form: Note that BQP problems can be considered as special cases of quadratically constrained quadratic program (QCQP), as the constraint x ∈ {1, −1} n is equivalent to x 2 i = 1, ∀i = 1, · · · , n. Problems over x ∈ {0, 1} n can be also expressed as {1, −1}problems (3) by replacing x with y = 2x − 1.
Solving (3) is in general NP-hard, so relaxation methods are considered in this paper. Relaxation to (3) can be done by extending the feasible set to a larger set, such that the optimal value of the relaxation is a lower bound on the optimal value of (3). The SDP relaxation to (3) can be expressed as: Note that constraint (4e) is equivalent to X − xx ∈ S n , which is the convex relaxation to the nonconvex constraint X − xx = 0. In other words, (4) is equivalent to (3), by replacing constraint (4e) with X = xx or by adding the constraint rank( 1 x x X ) = 1. The objective function and constraints (apart from the p.s.d. constraint) of the problem (4) are all linear with respect to X = 1 x x X , so (4) can be expressed in the homogenized form shown in (1) with respect to X. For simplicity, we consider the homogeneous problem (1), instead of (4), in the sequel.
Note that the SDP solution does not offer a feasible solution to the BQP (3) directly, unless it is of rank 1. A rounding procedure is required to extract a feasible BQP solution from the SDP solution, which will be discussed in Section 4.3.

SDCUT FORMULATION
A regularized SDP formulation is considered in this work: where γ > 0 is a prescribed parameter (its practical value is discussed in Section 5.1). Compared to (1), the formulation (5) adds into the objective function a Frobenius-norm term with respect to X. The reasons for choosing this particular formulation are two-fold: i) The solution quality of (5) can be as close to that of (4) as desired by making γ sufficiently large. ii) A simple dual formulation can be derived from (5), which can be optimized using quasi-Newton or inexact generalized Newton approaches.
In the following, a few desirable properties of (5) are demonstrated, where X denotes the optimal solution to (1) and X γ denotes the optimal solution to (5) with respect to γ. The proofs can be found in Section 7. Proposition 1. The following results hold: The above results show that the solution quality of (5) can be monotonically improved towards that of (4), by making γ sufficiently large.
Proposition 2. The dual problem of (5) can be simplified to where λ i , p i , i = 1, · · · , n are eigenvalues and the corresponding eigenvectors of C(u). Supposing problem (5) is feasible and denoting u as the dual optimal solution, we have: The simplified dual (6) is convex and contains only simple box constraints. Furthermore, its objective function d γ (·) has the following important properties. Proposition 3. d γ (·) is continuously differentiable but not necessarily twice differentiable, and its gradient is given by where Φ : Based on the above result, the dual problem can be solved by quasi-Newton methods directly. Furthermore, we also show in Section 4.2 that, the second-order derivatives of d γ (·) can be smoothed such that inexact generalized Newton methods can be applied.
2γ yields a lower-bound on the optimum of the BQP (3).
The above result is important as the lower-bound can be used to examine how close between an approximate binary solution and the global optimum.

Related Work
Considering the original SDP dual problem (2), we can find that its p.s.d. constraint, that is Other forms of penalty terms may be employed in the dual. The spectral bundle method of [52] penalizes λ max (X) and the log-barrier function is used in [53]. It is shown in [48] that these two first-order methods may converge slowly for some SDP problems. Note that the objective function of the spectral bundle methods is not necessarily differentiable (λ max (·) is differentiable if and only if it has multiplicity one). The objective function of our formulation is differentiable and its twice derivatives can be smoothed, such that classical methods can be easily used for solving our problems, using quasi-Newton and inexact generalized Newton methods.
Consider a proximal algorithm for solving SDP with only equality constraints (see [47], [48], [49], [50], [55]): where Φ[X] := [ B 1 , X , · · · , B m , X ] . Our algorithm is equivalent to solving the inner problem, that is, evaluating G γ (Y), with a fixed γ and Y = 0. In other words, our methods attempt to solve the original SDP relaxation approximately, with a faster speed. After rounding, typically, the resulting solutions of our algorithms are already close to those of the original SDP relaxation.
Our method is mainly motivated by the work of Shen et al. [56], which presented a fast dual SDP approach to Mahalanobis metric learning. They, however, focused on learning a real-valued metric for nearest neighbour classification. Here, in contrast, we are interested in discrete combinatorial optimization problems arising in computer vision. Krislock et al. [57] have independently formulated a similar SDP problem for the max-cut problem, which is simpler than the problems that we solve here. Moreover, they focus on globally solving the max-cut problem using branch-andbound.

SOLVING THE DUAL PROBLEM
Based on Proposition 3, first-order methods (for example gradient descent, quasi-Newton), which only require the calculation of the objective function and its gradients, can be directly applied to solving (6). It is difficult in employing standard Newton methods, however, as they require the calculation of second-order derivatives. In the following two sections, we present two algorithms for solving the dual (6), which are based on quasi-Newton and inexact generalized Newton methods respectively.

Quasi-Newton Methods
One main advantage of quasi-Newton methods over Newton methods is that the inversion of the Hessian matrix is approximated by analyzing successive gradient vectors, and thus that there is no need to explicitly compute the Hessian matrix and its inverse, which can be very expensive. Therefore the per-iteration computation cost of quasi-Newton methods is less than that of standard Newton methods. Algorithm 1 SDCut-QN: Solving (6) using quasi-Newton methods. Input: A, Φ, b, γ, u 0 , Kmax, τ > 0.
The quasi-Newton algorithm for (6) (referred to as SDCut-QN) is summarized in Algorithm 1. In Step 1, the dual problem (6) is solved using L-BFGS-B [58], which only requires the calculation of the dual objective function (6) and its gradient (8). At each iteration, a descent direction for ∆u is computed based on the gradient ∇d γ (u) and the approximated inverse of the Hessian matrix: H ≈ (∇ 2 d γ (u)) −1 . A step size ρ is found using line search. The algorithm is stopped when the difference between successive dual objective values is smaller than a pre-set tolerance.
After solving the dual using L-BFGS-B, the primal optimal variable X is calculated from the dual optimal u based on Equation (7) in Step 2.
Finally in Step 3, the primal optimal variable X is discretized and factorized to produce the feasible binary solution x , which will be described in Section 4.3.
Now we have an upper-bound and a lower-bound (see Propsition 4) on the optimum of the original BQP (3) (referred to as p ): 2γ . These two values are used to measure the solution quality in the experiments.

Smoothing Newton Methods
As d γ (u) is a concave function, the dual problem (6) is equivalent to finding u ∈ D such that u − u , −∇d γ (u ) ≥ 0, ∀u ∈ D, which is known as variational inequality [59]. D := R |Ieq| ×R |Iin| + is used to denote the feasible set of the dual problem. Thus (6) is also equivalent to finding a root of the following equation: where is continuous but not continuously differentiable, as both Π D and Π S n + have the same smoothness property. Therefore, standard Newton methods cannot be applied directly to solving (10). In this work, we use the inexact smoothing Newton method in [60] to solve the smoothed Newton equation: whereF( , u) is a smoothing function of F(u), which is constructed as follows.
Firstly, the smoothing functions for Π D and Π S n + are respectively written as: Algorithm 2 SDCut-SN: Solving (6) using smoothing Newton methods.
Step 1.2: Solve the following linear system up to certain accuracy Step 1.3: Line Search Step 2: u = u k+1 , X = γΠ S n + (C(u )).
Algorithm 3 Randomized Rounding Procedure: x = Round(X ) Input: The SDP solution X , which is decomposed to a set of vectors v 1 . . . vn ∈ R r where r = rank(X ). for k = 0, 1, 2, . . . , K do Step 1: Random sampling: obtain a real 1- Step 2: Discretization: z is discretized to a feasible BQP solution (see Table 2 for problem-specific methods). Output: x is assigned to the best feasible solution.
where λ i and p i are the ith eigenvalue and the corresponding eigenvector of X. φ( , v) is the Huber smoothing function that we adopt here to replace max(0, v): Note Lipschitz continuous on R, R × R m , R × S n respectively, and they are continuously differentiable when = 0. ThenF( , u) is defined as: which has the same smoothness property asΠ D andΠ S n + . The presented inexact smoothing Newton method (referred to as SDCut-SN) is shown in Algorithm 2. In Step 1.2, the Newton linear system (16) is solved approximately using conjugate gradient (CG) methods when |I in | = 0 and using biconjugate gradient stabilized (BiCGStab) methods [61] otherwise. In Step 1.3, we carry out a search in the direction [∆ k ; ∆u k ] for an appropriate step size ρ l such that the norm of E( , u) is decreased.

Randomized Rounding Procedure
In this section, we describe a randomized rounding procedure (see Algorithm 3) for obtaining a feasible binary solution from the relaxed SDP solution X .
Suppose that X is decomposed into a set of r-dimensional vectors v 1 . . . v n , such that X ij = v i v j . This decomposition can be easily obtained through the eigen-decomposition of X : We can see that these vectors reside on the r-dimensional unit sphere S r := {v ∈ R r , v v = 1}, and the angle between two vectors v i and v j defines how likely the corresponding two variables x i and x j will be separated (assigned with different labels). To transform these vectors into binary solutions, they are firstly projected onto a random 1dimensional line y ∼ N(0, I r ) in Step 1 of Algorithm 3, that is, z = [v 1 . . . v n ] y. Note that Step 1 is equivalent to sampling z from the Gaussian distribution N(0, X ), which has a probabilistic interpretation [62], [63]: X is the optimal solution to the problem where Σ denotes a covariance matrix. The proof is simple: since (17) is equivalent to (1). In other words, z solves the BQP in expectation. As the eigen-decomposition of X is already known when computing Π S n + (C(u )) at the last descent step, there is no extra computation for obtaining v 1 . . . v n . Due to the low-rank structure of SDP solutions (see Section 4.4), the computational complexity of sampling z is linear in the number of variables n.
Note that the above random sampling procedure does not guarantee that a feasible solution can always be found. In particular, this procedure will certainly fail when equality constraints are imposed on the problems [62]. But for all the problems considered in this work, each random sample z can be discretized to a "nearby" feasible solution (Step 2 of Algorithm 3). The discretization step is problem dependant, which is discussed in Table 2.

Speeding Up the Computation
In this section, we discuss several techniques for the eigendecompostion of C(u), which is one of the computational bottleneck for our algorithms. Low-rank Solution In our experiments, we observe that the final p.s.d. solution typically has a low-rank structure and r = rank(Π S n + (C(u))) usually decreases sharply such that r n for most of descent iterations in both our algorithms. Actually, it is known (see [64] and [65]) that any SDP problem with m linear constraints has an optimal solution X ∈ S n + , such that rank(X )(rank(X ) + 1)/2 ≤ m. It means that the rank of X is roughly bounded by √ 2m. Then Lanczos methods can be used to efficiently calculate the r positive eigenvalues of C(u) and the corresponding eigenvectors. Lanczos methods rely only on the product of the matrix C(u) and a column vector. This simple interface allows us to exploit specific structures of the coefficient matrices A and B i , i = 1, · · · , m. Specific Problem Structure In many cases, A and B i are sparse or structured. Such that the computational complexity and memory requirement of the matrix-vector product with respect to C(u) can be considered as linear in n, which are assumed as O(nt 1 ) and O(nt 2 ) respectively. The iterative Lanczos methods are faster than standard eigensolvers when r n and C(u) is sparse/structured, which require O(nr 2 + nt 1 r) flops and O(nr + nt 2 ) bytes at each iteration of Lanczos factorization, given that the number of Lanczos basis vectors is set to a small multiple (1 ∼ 3) of r. ARPACK [66], an implementation of Lanczos algorithms, is employed in this work for the eigen-decomposition of sparse or structured matrices. The DSYEVR function in LAPACK [67] is used for dense matrices.
Warm Start A good initial point is crucial for the convergence speed of iterative Lanczos methods. In quasi-Newton and smoothing Newton methods, the step size ∆u = u k+1 − u k tends to decrease with descent iterations. It means that C(u k+1 ) and C(u k ) may have similar eigenstructures, which inspires us to use a random linear combination of eigenvectors of C(u k ) as the starting point of the Lanczos process for C(u k+1 ). Parallelization Due to the importance of eigen-decomposition, its parallelization has been well studied and there are several off-theshelf parallel eigensolvers (such as SLEPc [68], PLASMA [69] and MAGMA [70]). Therefore, our algorithms can also be easily parallelized by using these off-the-shelf parallel eigensolvers.

Convergence Speed, Computational Complexity and Memory Requirement
SDCut-QN In general, quasi-Newton methods converge superlinearly given that the objective function is at least twice differentiable (see [71], [72], [73]). However, the dual objective function in our case (6) is not necessarily twice differentiable. So the theoretical convergence speed of SDCut-QN is unknown.
At each iteration of L-BFGS-B, both of the computational complexity and memory requirement of L-BFGS-B itself are O(m). The only computational bottleneck of SDCut-QN is on the computation of the projection Π S n + (C(u)), which is discussed in Section 4.4. SDCut-SN The inexact smoothing Newton method SDCut-SN is quadratically convergent under the assumption that the constraint nondegenerate condition holds at the optimal solution (see [60]). There are two computationally intensive aspects of SDCut-SN: i). the CG algorithms for solving the linear system (16). In the appendix, we show that the Jacobian-vector product requires O(m + n 2 r) flops at each CG iteration, where r = rank(Π S n + (C(u))). ii). All eigenpairs of C(u) are needed to obtain Jacobian matrices implicitly, which takes O(n 3 ) flops using DSYEVR function in LAPACK.
From Table 1, we can see that the computational costs and memory requirements for both SDCut-QN and SDCut-SN are linear in m, which means that our methods are much more scalable to m than interior-point methods. In terms of n, our methods is also more scalable than interior-point methods and comparable to spectral methods. Especially for sparse/structured matrices, the computational complexity of SDCut-QN is linear in n. As SDCut-SN cannot significantly benefit from sparse/structured matrices, it needs more time than SDCut-QN in each descent iteration for such matrices. However, SDCut-SN has a fast convergence rate than SDCut-QN. In the experiment section, we compare the speeds of SDCut-SN and SDCut-QN in different cases.

APPLICATIONS
We now show how we can attain good solutions on various vision tasks with the proposed methods. The two proposed methods, SDCut-QN and SDCut-SN, are evaluated on several computer vision applications. The BQP formulation of different applications and the corresponding rounding heuristics are demonstrated in Table 2. The corresponding SDP relaxation can be obtained based on (4). In the experiments, we also compare our methods to spectral methods [19], [20], [21], [22], graph cuts based methods [7], [8], [9] and interior-point based SDP methods [41], [42], [43]. The upper-bounds (that is, the objective value of BQP solutions) and the lower-bounds (on the optimal objective value of BQPs) achieved by different methods are demonstrated, and the runtimes are also compared. The code is written in Matlab, with some key subroutines implemented in C/MEX. We have used the L-BFGS-B [58] for the optimization in SDCut-QN. All of the experiments are evaluated on a core of Intel Xeon E5-2680 2.7GHz CPU (20MB cache). The maximum number of descent iterations of SDCut-QN and SDCut-SN are set to 50 and 500 respectively. As shown in Algorithm 1 and Algorithm 2, the same stopping criterion is used for SDCut-QN and SDCut-SN, and the tolerance τ is set to 10 7 eps where eps is the machine precision. The initial values of the dual variables u i , i ∈ I eq are set to 0, and u i , i ∈ I in are set to a small positive number. The selection of parameter γ will be discussed in the next section.

Graph Bisection
Graph bisection is a problem of separating the nodes of a weighted graph G = (V, E) into two disjoint sets with equal cardinality, while minimizing the total weights of the edges being cut. V denotes the set of nodes and E denotes the set of non-zero edges. The BQP formulation of graph bisection can be found in (18) of Table 2. To enforce the feasibility (two partitions with equal size), the randomized score vector z in Algorithm 3 is dicretized by thresholding the median value (see Table 2).
To show that the proposed SDP methods have better solution quality than spectral methods we compare the graph-bisection results of RatioCut [76], Normalized Cut (NCut) [19] and SDCut-QN on two artificial 2-dimensional datasets.
As shown in Fig. 1, the first data set (the first row) contains two sets of points with different densities, and the second set contains an outlier. RatioCut and NCut fail to offer satisfactory results on both of the data sets, possibly due to the poor approximation of spectral relaxation. In contrast, our SDCut-QN achieves desired results on these data sets.
Secondly, to demonstrate the impact of the parameter γ, we test SDCut-QN and SDCut-SN on a random graph with γ ranging from 10 2 to 10 4 (A and B i in (1) are scaled such that A 2 F = B 2 F = 1). The graph is generated with 1000 vertices and all possible edges are assigned a non-zero weight uniformly sampled from (0, 1]. As the resulting affinity matrices are dense, the DSYEVR routine in LAPACK package is used for eigen-decomposition. In Fig. 2 otherwise, where dij denotes the Euclidean distance between i and j. Discretization: x = sign(z − median(z)).
Image segmentation with partial grouping constraints (Sec. 5.2) where fi denotes the local feature of pixel i. The weighted partial grouping pixels are defined as s f = Pt f /(1 Pt f ) and s b = Pt b /(1 Pt b ) for foreground and background respectively, where t f , t b ∈ {0, 1} n are two indicator vectors for manually labelled pixels and P = Diag(W1) −1 W is the normalized affinity matrix used as smoothing terms [20]. The overlapped non-zero elements between s f and s b are removed. κ ∈ (0, 1] denotes the degree of belief. Discretization: see (25).
W is the affinity matrix as defined above. q is the target K-bin color histogram; ti ∈ {0, 1} n is the indicator vector for every color bin; δ is the prescribed upper-bound on the Euclidean distance between the obtained histogram and q. Note that (20b) is equivalent to a quadratic constraint on x and can be expressed as xBx + a x ≤ b. Constraint (20b) is penalized in the objective function with a weight (multiplier) α > 0 in this work: The definition of A can be found in [38]. s is the number of images, ni is the number of pixels for i-th image, and n = s i=1 ni. ti ∈ {0, 1} n is the indicator vector for the i-th image. κ ∈ (0, 1]. Discretization: see (27).
The unary terms (h ∈ R n ) and pairwise terms (H ∈ R n×n ) are learned using decision tree fields [75].
Discretization: x = sign(z). lower-bounds, number of iterations and time achieved by SDCut-QN and SDCut-SN, with respect to different values of γ. There are several observations: i) With the increase of γ, upper-bounds become smaller and lower-bounds become larger, which implies a tighter relaxation. ii) Both SDCut-QN and SDCut-SN take more iterations to converge when γ is larger. iii) SDCut-SN uses fewer iterations than SDCut-QN. The above observations coincide with the analysis in Section 4.5. Using a larger parameter γ yields better solution quality, but at the cost of slower convergence speed. The choice of a good γ is data dependant. To reduce the difficulty of the choice of γ, the matrices A and B i of Equation (1) are scaled such that the Frobenius norm is 1 in the following experiments.
Thirdly, experiments are performed to evaluate another two factors affecting the speed of our methods: the sparsity of the affinity matrix W and the matrix size n. The numerical results corresponding to dense and sparse affinity matrices are shown in Table 3 and Table 4 respectively. The sparse affinity matrices are generated from random graphs with 8-neighbour connection. In these experiments, the size of matrix W is varied from 200 to 5000. ARPACK is used by SDCut-QN for partial eigendecomposition of sparse problems, and DSYEVR is used for other cases. For both SDCut-QN and SDCut-SN, the number of iterations does not grow significantly with the increase of n. However, the running time is still correlated with n, since an eigen-decompostion of an n × n matrix needs to be computed at each iteration for both of our methods. We also find that the second-order method SDCut-SN uses significantly fewer iterations than the first-order method SDCut-QN. For dense affinity matrices, SDCut-SN runs consistently faster than SDCut-QN. In contrast for sparse affinity matrices, SDCut-SN is only faster than SDCut-QN on problems of size up to n ≥ 2000. That is because the Lanczos method used by SDCut-QN (for partial eigen-decompostion) scales much better for large sparse matrices than the standard factorization method (DSYEVR) used by SDCut-SN (for full eigen-    Finally, we evaluate SDCut-QN on a large dense graph with 10000 nodes. The speed performance is compared on a single CPU core (using DSYEVR function of LAPACK as eigensolver) and a hybrid CPU+GPU workstation (using the DSYEVDX 2STAGE function of MAGMA as eigensolver). The results are shown in Table 5 and we can see that the parallelization brings a 10-fold speedup over running on a single CPU core. The lower-/upper-bounds are almost identical as there is no difference apart from the implementation of eigen-decompostion.

Constrained Image Segmentation
We consider image segmentation with two types of quadratic constraints (with respect to x): partial grouping constraints [20] and histogram constraints [77]. The affinity matrix W is sparse, so ARPACK is used by SDCut-QN for eigen-decomposition.
Besides interior-point SDP methods, we also compare our methods with graph-cuts [7], [8], [9] and two constrained spectral clustering method proposed by Maji et al. [78] (referred to as BNCut) and Wang and Davidson [79] (referred to as SMQC). BNCut and SMQC can encode only one quadratic constraint, but it is difficult (if not impossible) to generalize them to multiple quadratic constraints. Partial Grouping Constraints The corresponding BQP formulation is Equation (19) Fig. 4. SDCut-SN uses fewer iterations than SDCut-QN and is faster than all other SDP based methods. Graph cuts and SMQC exhibit worse F-measure scores than SDP based methods.
obtained from any random sample z as follows: where Fig. 3 illustrates the result for image segmentation with partial grouping constraints on the Berkeley dataset [80]. All the test images are over-segmented into about 760 superpixels. We find that BNCut did not accurately segment foreground, as it only incorporates a single set of grouping pixels (foreground). In contrast, our methods are able to accommodate multiple sets of grouping pixels and segment the foreground more accurately. In Table 6, we compare the CPU time and the upper-bounds of SDCut-QN, SeDuMi and SDPT3. SDCut-QN achieves objective values similar to that of SeDuMi and SDPT3, yet is over 10 times faster. Histogram Constraints Given a random sample z, a feasible solution x to the corresponding BQP formulation (20) can be obtained through andz is obtained by sorting z in descending order. For graph cuts methods, the histogram constraint is encoded as unary terms: ϕ i = − ln Pr(f i |fore)/Pr(f i |back) , i = 1, 2, . . . , n. Pr(f i |fore) and Pr(f i |back) are probabilities for the color of the ith pixel belonging to foreground and background respectively. Fig. 4 and Table 7 demonstrate the results for image segmentation with histogram constraints. We can see that unary terms (the second row in Fig. 4) are not ideal especially when the color distribution of foreground and background are overlapped. For example in the first image, the white collar of the person in the foreground have similar unary terms with the white wall in the background. The fourth row of Fig. 4 shows that the unsatisfactory unary terms degrade the segmentation results of graph cuts methods significantly.
The average F-measure of all evaluated methods are reported in Table 7. Our methods outperforms graph cuts and SMQC in terms of F-measure. As for the running time, SDCut-SN is faster    than all other SDP-based methods (that is, SDCut-QN, SeDuMi, SDPT3 and MOSEK). As expected, SDCut-SN uses much less (1/6) iterations than SDCut-QN. SDCut-QN and SDCut-SN have comparable upper-bounds and lower-bounds than interior-point methods.
From Table 7, we can find that SMQC is faster than our methods. However, SMQC does not scale well to large problems since it needs to compute full eigen-decomposition. We also test SDCut-QN and SMQC on problems with a larger number of superpixels (9801). Both of the algorithms achieve similar segmentation results, but SDCut-QN is much faster than SMQC (23m21s vs. 4h9m).

Image Co-segmentation
The task of image co-segmentation [38] aims to partition a common object from multiple images simultaneously. In this work, the Weizman horses 1 and MSRC 2 datasets are tested. There are 6 ∼ 10 images in each of four classes, namely "car-front", "car-back", "face" and "horse". Each image is oversegmented to 400 ∼ 700 superpixels. The number of binary variables n is then increased to 4000 ∼ 7000. The BQP formulation for image co-segmentation can be found in Table 2 (see [38] for details). The matrix A can be decomposed into a sparse matrix and a structural matrix, such that ARPACK can be used by SDCut-QN. Each vector z = [z (1) , · · · , z (K) ] (where z i corresponds to the i-th image) randomly sampled from N(0, X ) is discretized to a feasible BQP solution as follows: where θ (i) can be obtained as (26). We compare our methods with the low-rank factorization method [38] (referred to as LowRank) and interior-point methods. As we can see in Table 8, SDCut-QN takes 10 times more iterations than SDCut-SN, but still runs faster than SDCut-SN especially when the size of problem is large (see "face" data). The reason is that SDCut-QN can exploit the specific structure of matrix A in eigen-decomposition. SDCut-QN runs also 5 times faster than LowRank. All methods provide similar upper-bounds (primal objective values), and the score vectors shown in Fig. 5 also show that the evaluated methods achieve similar visual results.

Graph Matching
In the graph matching problems considered in this work, each of the K source points must be matched to one of the L target points, where L ≥ K. The optimal matching should maximize both of the local feature similarities between matched-pairs and the structure similarity between the source and target graphs.
The BQP formulation of graph matching can be found in Table 2, which can be relaxed to: L j=1 x (i−1)L+j = 1, ∀i ∈ K, (28c) X (i−1)L+j,(i−1)L+k = 0, ∀i ∈ K, j = k ∈ L, (28d) X (j−1)L+i,(k−1)L+i = 0, ∀i ∈ L, j = k ∈ K, (28e) where K = {1, · · · , K} and L = {1, · · · , L}. A feasible binary solution x is obtained by solving the following linear program (see [37] for details): Two-dimensional points are randomly generated for evaluation. Table 9 shows the results for different problem sizes: n ranges from 201 to 3201 and m ranges from 3011 to 192041. SDCut-SN and SDCut-QN achieves exactly the same upper-bounds as interior-point methods and comparable lower-bounds. Regarding the running time, SDCut-SN takes much less number of iterations to converge and is relatively faster (within 2 times) than SDCut-QN. Our methods run significantly faster than interior-point methods. Taking the case K × L = 25 × 50 as an example, SDCut-SN and SDCut-QN converge at around 4 minutes and interiorpoint methods do not converge within 24 hours. Furthermore, interior-point methods runs out of 100G memory limit when the number of primal constraints m is over 10 5 . SMAC [22], a spectral     method incorporating affine constrains, is also evaluated in this experiment, which provides worse upper-bounds and error ratios.

Image Deconvolution
Image deconvolution with a known blurring kernel is typically equivalent to solving a regularized linear inverse problem (see (23) in Table 2). In this experiment, we test our algorithms on two binary 30 × 75 images blurred by an 11 × 11 Gaussian kernel. LP based methods such as TRWS [12] and MPLP [13] are also evaluated. Note that the resulting models are difficult for graph cuts or LP relaxation based methods in that it is densely connected and contains a large portion of non-submodular pairwise potentials. We can see from Fig. 6 and Table 10 that QPBO [7], [8], [9] leaves most of pixels unlabelled and LP methods (TRWS and MPLP) achieves worse segmentation accuracy. SDCut-SN achieves a 10fold speedup over interior-point methods while keep comparable upper-/lower-bounds. Using much less runtime, SDCut-SN still yields significantly better upper-/lower-bounds than LP methods.

Chinese Character Inpainting
The MRF models for Chinese character inpainting are obtained from the OpenGM benchmark [82], in which the unary terms and pairwise terms are learned using decision tree fields [75]. As there are non-submodular terms in these models, they cannot be solved exactly using graph cuts. In this experiments, all models are firstly reduced using QPBO and different algorithms are compared on the 100 reduced models. Our approach is compared to LPbased methods, including TRWS, MPLP. From the results shown in Table 11, we can see that SDCut-SN runs much faster than interior-point methods (SeDuMi and MOSEK) and has similar upper-bounds and lower-bounds. SDCut-SN is also better than TRWS and MPLP in terms of upper-bound and lower-bound. An extension of MPLP (refer to as MPLP-C) [81], [14], which adds violated cycle constraints iteratively, is also evaluated in this experiment. In MPLP-C, 1000 LP iterations are performed initially and then 20 cycle constraints are added at every 20    LP iteration. MPLP-C performs worse than SDCut-SN under the runtime limit of 5 minutes, and outperforms SDCut-SN with a much longer runtime limit (1 hour). We also find that SDCut-SN achieves better lower-bounds than MPLP-C on the instances with more edges (pairwise potential terms). Note that the time complexity of MPLP-C (per LP iteration) is proportional to the number of edges, while the time complexity of SDCut-SN (see Table 1) is less affected by the edge number. It should be also noticed that SDCut-SN uses much less runtime than MPLP-C, and its bound quality can be improved by adding linear constraints (including cycle constraints) as well [83].

CONCLUSION
In this paper, we have presented a regularized SDP algorithm (SD-Cut) for BQPs. SDCut produces bounds comparable to the conventional SDP relaxation, and can be solved much more efficiently. Two algorithms are proposed based on quasi-Newton methods (SDCut-QN) and smoothing Newton methods (SDCut-SN) respectively. Both SDCut-QN and SDCut-SN are more efficient than classic interior-point algorithms. To be specific, SDCut-SN is faster than SDCut-QN for small to medium sized problems. If the matrix to be eigen-decomposed, C(u), has a special structure (for example, sparse or low-rank) such that matrix-vector products can be computed efficiently, SDCut-QN is much more scalable to large problems. The proposed algorithms have been applied to several computer vision tasks, which demonstrate their flexibility in accommodating different types of constraints. Experiments also show the computational efficiency and good solution quality of SDCut. We have made the code available online 3 .
As a pointwise minimum of affine functions of t, θ(t) is concave and continuous. It is also easy to find that θ(0) = 0 and θ(t) is monotonically increasing on R + . So for any > 0, there is a t > 0 (and equivalently γ = 1 2t > 0) such that |p(X )−p(X γ )| < .  (ii) By the definition of X γ1 and X γ2 , it is clear that

Proof of Proposition 2
Proof. The Lagrangian of the primal problem (5) is: where u ∈ R |Ieq| ×R |Iin| + and Z ∈ S n + are Lagrangian multipliers. Supposing (5) and (31) are feasible, strong duality holds and ∇ X L(X , u , Z ) = 0, where X , u and Z are optimal solutions. Then we have that By substituting X to (31), we have the dual: The variable Z can be further eliminated as follows. Given a fixed u, (33) can be simplified to: min Z∈S n + Z + C(u) 2 F , which is proved to has the solution Z = Π S n + (−C(u)) (see [84] or Section 8.1.1 of [85]). Note that C(u) = Π S n + (C(u)) − Π S n + (−C(u)), so Z + C(u) = Π S n + (C(u)). By substituting Z into (33) and (32), we have the simplified dual (6) and Equation (7).

The Spherical Constraint
Before proving Proposition 4, we first give the following theorem.
Theorem 5. (The spherical constraint). For any X ∈ S n + , we have the inequality X F ≤ trace(X), in which the equality holds if and only if rank(X) = 1.
Proof. The proof given here is an extension of the one in [89]. We have X 2 F = trace(XX ) = n i=1 λ 2 i ≤ (trace(X)) 2 , where λ i ≥ 0 denotes the i-th eigenvalue of X. Note that , if and only if there is only one non-zero eigenvalue of X, that is, rank(X) = 1.

Proof of Proposition 4
Proof. Firstly, we have the following inequalities: where the second inequality is based on Theorem 5. For the BQP (3) that we consider, it is easy to see trace(X ) = n. Furthermore, p γ (X ) ≥ p γ (X γ ) holds by definition. Then we have that It is known that the optimum of the original SDP problem (4) is a lower-bound on the optimum of the BQP (3) (denoted by p ): p(X ) ≤ p . Then according to (35), we have p γ (X γ ) − n 2 2γ ≤ p(X ) ≤ p . Finally based on the strong duality, the primal objective value is not smaller than the dual objective value in the feasible set (see for example [85]

APPENDIX
In this section, we present some computational details.

A.1.1 Euclidean Projection onto the P.S.D. Cone
Theorem A.1. The Euclidean projection of a symmetric matrix X ∈ S n onto the positive semidefinite cone S n + , is given by where λ i , p i , i = 1, · · · , n are the eigenvalues and the corresponding eigenvectors of X.
Proof. This result is well-known and its proof can be found in [84] or Section 8.1.1 of [85].

A.1.2 Derivatives of Separable Spectral Functions
A spectral function F(X) : S n → R is a function which depends only on the eigenvalues of a symmetric matrix X, and can be written as f(λ) for some symmetric function f : R n → R, where λ = [λ i , · · · , λ n ] denotes the vector of eigenvalues of X. A function f(·) is symmetric means that f(x) = f(Ux) for any permutation matrix U and any x in the domain of f(·). Such symmetric functions and the corresponding spectral functions are called separable, when f(x) = n i=1 g(x i ) for some function g : R → R. It is known (see [86], [87], [88], for example) that a spectral function has the following properties: Theorem A.2. A separable spectral function F(·) is k-times (continuously) differentiable at X ∈ S n , if and only if its corresponding function g(·) is k-times (continuously) differentiable at λ i , i = 1, . . . , n, and the first-and second-order derivatives of F(·) are given by ∇F(X) = P diag ∇g(λ 1 ), ∇g(λ 2 ), . . . , ∇g(λ n ) P , if λ i = λ j , i, j = 1, . . . , n. λ = [λ 1 , · · · , λ n ] and P = [p 1 , · · · , p n ] are the collection of eigenvalues and the corresponding eigenvectors of X.

A.2.1 Smoothing Function
In this section, we show how to constrcut a smoothing function of F(u) (see (10)). First, the smoothing functions for Π D and Π S n + are written as follows respectively:Π where λ i and p i are the ith eigenvalue and the corresponding eigenvector of X. φ( , v) is the Huber smoothing function that we adopt here to replace max(0, v): if v < −0.5 . which has the same smooth property asΠ D andΠ S n + . (16) The linear system (16) can be decomposed to two parts:

A.2.2 Solving the Linear System
where ∇ F and ∇ uF denote the partial derivatives ofF with respect to and u respectively. One can firstly obtain the value of ∆ k by (A.8b) and then solve the linear system (A.8c) using CG-like algorithms.
Since the Jacobian matrix ∇ uF ( k , u k ) ∈ R m×m is nonsymmetric when inequality constraints exist, biconjugate gradient stabilized (BiCGStab) methods [61] are used for (A.8c) with respect to |I in | = 0, and classic conjugate gradient methods are used when |I in | = 0.
The computational bottleneck of CG-like algorithms is on the Jacobian-vector products at each iteration. We discuss in the following the computational complexity of it in our specific cases. Firstly, we give the partial derivatives of smoothing functions φ( , v) : R×R → R,Π D ( , v) : R × R m → R m andΠ S n + ( , X) : R × S n → S n : In general cases, computing (A.13a) and (A.13b) needs O(mn 2 + n 3 ) flops. However, based on the observation that most of B i , i = 1, . . . , m contain only O(1) elements and r = rank(Π S n + (C(u))) n, the computation cost can be dramatically reduced. Firstly, super sparse B i s lead to the computation cost of Φ and Ψ reduced from O(mn 2 ) to O(m+n). Secondly, note that [Ω( , λ)] ij = 0, ∀λ i , λ j < 0. Given r n and is small enough, the matrix Ω only contains non-zero elements in the first r columns and rows. Thus the matrix multiplication in (A.13a), (A.13b) and (A.7) can be computed in O(n 2 r) flops rather than the usual O(n 3 ) flops.
In summary, the computation cost of the right hand side of Equ. (A.8c) and the Jacobian-vector product (A.13b) can be reduced from O(mn 2 + n 3 ) to O(m + n 2 r) in our cases.