# Fast Large-Scale Discrete Optimization Based on Principal Coordinate Descent

Huan Xiong &Mengyang Yu &Li Liu &Fan Zhu &Fumin Shen &Ling Shao

###### Abstract

Binary optimization, a representative subclass of discrete optimization, plays
an important role in mathematical optimization and has various applications in
computer vision and machine learning. Usually, binary optimization problems
are NP-hard and difficult to solve due to the binary constraints, especially
when the number of variables is very large. Existing methods often suffer from
high computational costs or large accumulated quantization errors, or are only
designed for specific tasks. In this paper, we propose a fast algorithm to
find effective approximate solutions for general binary optimization problems.
The proposed algorithm iteratively solves minimization problems related to the
linear surrogates of loss functions, which leads to the updating of some
binary variables most impacting the value of loss functions in each step. Our
method supports a wide class of empirical objective functions with/without
restrictions on the numbers of $1$s and $-1$s in the binary variables.
Furthermore, the theoretical convergence of our algorithm is proven, and the
explicit convergence rates are derived, for objective functions with Lipschitz
continuous gradients, which are commonly adopted in practice. Extensive
experiments on several binary optimization tasks and large-scale datasets
demonstrate the superiority of the proposed algorithm over several state-of-
the-art methods in terms of both effectiveness and efficiency.

##  1 Introduction

Binary optimization problems are generally formulated as follows:

$\displaystyle\min_{\mathbf{x}}~{}f(\mathbf{x}),\,\,\,\,\,\,\,\,\text{s.t.}~{}~{}\mathbf{x}\in\\{\pm
1\\}^{n}.$ (1)

Problem (1) appears naturally in several fields of computer vision and machine
learning, including clustering [36], graph bisection [37, 43], image denoising
[5], dense subgraph discovery [2, 4, 43, 46], multi-target tracking [34], and
community discovery [12]. In many application scenarios, such as binary
hashing [11, 21, 32, 33, 35, 39, 42], Problem (1) needs to be solved for
millions of binary variables, which makes the size $2^{n}$ of the feasible set
very large (far larger than the number of atoms in the universe). Usually, it
is difficult to find the optimal solution. Therefore, providing a fast
algorithm to approximately solve Problem (1) is very important in practice.

Furthermore, additional constraints on the numbers of $1$s and $-1$s in the
binary variables $\mathbf{x}$ are adopted in many cases. For example, in
binary hashing [32, 33] and graph bisection [37, 43], the balance condition is
often required, which means that the numbers of $1$s and $-1$s are equal to
each other. On the other hand, dense subgraph discovery [2, 4, 43, 46] and
information theoretic clustering [36] require that the numbers of $1$s and
$-1$s in $\mathbf{x}$ are some fixed integers.

To handle these previously mentioned constraints, in this paper, we focus on
the following binary optimization problem:

$\displaystyle\min_{\mathbf{x}}~{}f(\mathbf{x}),\,\,\,\,\,\,\,\,\text{s.t.}~{}~{}\mathbf{x}\in\sl{\Omega_{r}},$
(2)

where $\mathbf{x}$ is a binary vector of length $n$, $f(\cdot)$ is a
differentiable objective function (which may be nonconvex), $r\geq-1$ is a
given integer, and the restriction $\Omega_{r}$ on $\mathbf{x}$ is defined as

$\displaystyle\Omega_{r}=\begin{cases}\\{\pm 1\\}^{n},\ \ \
\,\text{if}\,\,\,r=-1;\\\ \\{\mathbf{x}\in\\{\pm
1\\}^{n}:\mathbf{1}^{\top}\mathbf{x}=2r-n\\},\ \ \
\,\text{if}\,\,\,r\in\mathbb{N}_{\geq 0},\end{cases}$ (3)

where $\mathbb{N}_{\geq 0}$ denote the set of nonnegative integers. When
$r=-1$, Problem (2) is a binary optimization problem without further
constraints. When $r\in\mathbb{N}_{\geq 0}$, Problem (2) becomes an
optimization problem with the restriction that there are exactly $r$ $1$s in
the binary vector $\mathbf{x}$. For instance, when $r=n/2$, the constraint
$\mathbf{1}^{\top}\mathbf{x}=2r-n=0$ implies that the number of $1$s is equal
to the number of $-1$s in $\mathbf{x}$.

In general, Problem (2) is NP-hard due to the binary constraints [15]. Many
algorithms, such as continuous relaxation, equivalent optimization, signed
gradient optimization and direct discrete optimization, have been proposed to
solve it approximately (for details, please refer to Section 2). However, they
usually suffer from high computational costs or large accumulated quantization
errors, or are only designed for specific tasks.

To overcome these difficulties, in this paper, we propose a novel and fast
optimization algorithm, termed Discrete Principal Coordinate Descent (DPCD),
to approximately solve Problem (2). The time complexity for the binary
optimization problem is relatively high when directly applying signed gradient
methods (updating all variables at each time based on the gradients). On the
contrary, the proposed DPCD focuses on the principal coordinates most
impacting the value of the loss function. At each iteration, DPCD can
adaptively decide the number of principal coordinates that need to be
optimized, which can be regarded as analogous to the adaptive learning rate in
the normal gradient descent methods. Different from other binary optimization
algorithms in the literature, our DPCD method supports a large family of
empirical objective functions with/without restrictions on the numbers of $1$s
and $-1$s in the binary variables. Furthermore, we prove theoretical
convergence of DPCD for loss functions with Lipschitz continuous gradients,
which cover almost every loss function in practice. Explicit convergence rates
are also derived. Extensive experiments on two binary optimization tasks:
dense subgraph discovery and binary hashing, show the superiority of our
method over state-of-the-art methods in terms of both efficiency and
effectiveness.

##  2 Related work

A very rich literature and a wide range of promising methods exist in binary
optimization. We briefly review three classes of representative and related
methods.

#### Continuous relaxation methods.

An intuitive method to approximately solve Problem (1) is to relax the binary
constraints to continuous variables, then threshold the continuous solutions
to binary vectors. For instance, in the Linear Programming (LP) relaxation
[13, 16], the binary constraint is substituted with the box constraint, i.e.,
$\mathbf{x}\in[-1,1]^{n}$, which can be approximately solved by continuous
optimization methods such as the interior-point method [24]. On the other
hand, the Semi-Definite Programming (SDP) relaxation [37] replaces the binary
constraint with some positive semi-definite matrix constraint. In Spectral
relaxation [19, 27], the binary constraint is relaxed to some $\ell_{2}$-ball,
which is non-convex. One of the advantages of such continuous relaxation
methods is that the relaxed problems can be approximately solved efficiently
by existing continuous optimization solvers. However, the relaxation is
usually too loose, and the thresholding often yields large quantization
errors.

#### Equivalent optimization methods.

Unlike relaxation methods, equivalent optimization methods replace the binary
constraint with some equivalent forms, which are much easier to handle. For
example, motivated by linear and spectral relaxations, Wu and Ghanem [41]
replaced the binary constraint with the intersection of the box $[-1,1]^{n}$
and the sphere $\\{\mathbf{x}:\|\mathbf{x}\|_{2}^{2}=n\\}$, and then applied
the Alternating Direction Method of Multipliers (ADMM) [7, 17, 38] to solve
the optimization problem iteratively. Other methods in this direction include
the MPEC-ADM and MPEC-EPM methods ([43, 45]), the $\ell_{0}$ norm
reformulation [22, 44], the $\ell_{2}$ box non-separable reformulation [25],
and the piecewise separable reformulation [48]. Usually, these equivalent
optimization methods guarantee the convergence to some stationary and feasible
points, but the convergence speed is often too slow, resulting in high
computational costs for large-scale optimization problems.

#### Signed gradient methods.

In the Signed Gradient Method (SGM) [21], a linear surrogate of the objective
function $f(\mathbf{x})$ is given at each iteration. Then, the minimization
(actually a maximization problem was studied in the original paper [21], we
state an equivalent form here) of this surrogate function gives the updating
rule for Problem (1) as: $\mathbf{x}^{k+1}=-{\text{sgn}}(\nabla
f(\mathbf{x}^{k}))$. The sequence obtained by this updating rule is guaranteed
to converge if the objective function is concave. However, even for a very
simple non-concave function, SGM may generate a divergent sequence and never
converge (please refer to Lemma 1). Furthermore, SGM cannot handle binary
problems with restriction on the number of $1$s since the number of $1$s may
change during each iteration. A stochastic version of this method was given in
Adaptive Discrete Minimization (ADM) [20], in which an adaptive ratio $\psi$
was selected at each iteration, then some random $\psi n$ entries of
$\mathbf{x}$ were updated by
$\mathbf{x}^{k+1}_{i}=-{\text{sgn}}(\nabla_{i}f(\mathbf{x}^{k}))$. Although
ADM works well for certain loss functions, it fails when the value of the loss
function depends largely on only a few variables, since the random selecting
procedure may skip such important variables.

#### Discrete optimization methods.

In the field of image hashing, many direct discrete optimization methods, such
as DCC [32], SADH [33], ARE [14], ITQ [10] and FastHash [18], have been
proposed, which aim to optimize binary variables directly (SGM can also be
seen as a discrete optimization method). For example, the Coordinate Descent
(CD) method [40] is widely used for solving optimization problems with smooth
and convex constraints. Motivated by this method, several discrete cyclic
coordinate descent (DCC) methods (e.g., RDCM [23], FSDH [11], and SDH [32])
have been proposed to handle the binary constraint directly. The main idea is
that, at each iteration, a subproblem with most entries of the binary
variables fixed is considered, and the loss function is minimized with respect
to the remaining entries. Although such methods can work well for specific
loss functions, most of them are usually difficult to be extended to handle
general binary optimization problems. Furthermore, they often suffer from
expensive computational costs.

##  3 Proposed algorithm

In this section, we present in detail the DPCD algorithm for solving Problem
(2), which is a general binary optimization problem with/without restrictions
on the numbers of $1$s and $-1$s. We also provide a theoretical convergence
analysis of the proposed method.

###  3.1 Notation and preliminaries

We first introduce some notation and preliminaries. A vector is represented by
some lowercase bold character, while a matrix is represented by some uppercase
bold character. Let $\mathbf{x}_{i}$ and $\mathbf{A}_{ij}$ denote the $i$-th
and $(i,j)$-th entries of a vector $\mathbf{x}$ and a matrix $\mathbf{A}$,
respectively. The transpose of a matrix $\mathbf{A}$ is represented by
$\mathbf{A}^{\top}$. We use $\langle\cdot,\cdot\rangle$ to denote the
Euclidean inner product. Let
$\|\mathbf{A}\|=\sqrt{\sum_{ij}\mathbf{A}_{ij}^{2}}$ and
$\|\mathbf{A}\|_{1}=\sum_{ij}|\mathbf{A}_{ij}|$ be the Frobenius norm and
$1$-norm of a matrix $\mathbf{A}$, respectively. The gradient of a
differentiable function $f(\mathbf{x})$ is denoted by $\nabla
f(\mathbf{x})=(\nabla_{1}f(\mathbf{x}),\nabla_{2}f(\mathbf{x}),\cdots,\nabla_{n}f(\mathbf{x}))$.
Let
${\text{sgn}}(\mathbf{x})=({\text{sgn}}(\mathbf{x}_{1}),{\text{sgn}}(\mathbf{x}_{2}),\cdots,{\text{sgn}}(\mathbf{x}_{n}))$
denote the element-wise sign function where ${\text{sgn}}(\mathbf{x}_{i})=1$
for $\mathbf{x}_{i}\geq 0$ and $-1$ otherwise. The Hamming distance between
two binary vectors $\mathbf{y}$ and $\mathbf{z}$ of equal length is defined by
$d_{H}(\mathbf{y},\mathbf{z})$, which is the number of positions at which the
corresponding entries are different. For a set $S$, let $\\#S$ denote the
number of elements in $S$.

###  3.2 Main algorithm

The proposed DPCD algorithm runs iteratively between principal coordinate
update and neighborhood search.

#### Principal coordinate update.

The basic idea is that, in the $k$-th iteration, we change the sign of some
adaptively chosen entries of the binary vector $\mathbf{x}^{k}$, such that the
value of the loss function should decrease steeply after each change. To
achieve this goal, we focus on $L$-principal coordinates (see the following
definition), which have major influences on the value change of the loss
function.

######  Definition 1.

Let $f(\mathbf{x})$ be a differentiable function and $L>0$ be a positive
constant. A coordinate index $i$ is called an $L$-principal coordinate of
$\mathbf{x}\in\\{-1,1\\}^{n}$ if the product
$\mathbf{x}_{i}\cdot\nabla_{i}f(\mathbf{x})\geq L$.

One motivation of our method is SGM [21]. In the $k$-th iteration of SGM, the
linear surrogate of the objective function $f(\mathbf{x})$ is given as:

$\displaystyle\hat{f}^{k}(\mathbf{x})=f(\mathbf{x}^{k})+\langle\nabla
f(\mathbf{x}^{k}),\mathbf{x}-\mathbf{x}^{k}\rangle.$ (4)

Then $\mathbf{x}^{k+1}$ is obtained by minimizing this surrogate function:

$\displaystyle\mathbf{x}^{k+1}=\arg\min_{\mathbf{x}\in\\{\pm
1\\}^{n}}\hat{f}^{k}(\mathbf{x})=-{\text{sgn}}(\nabla f(\mathbf{x}^{k})).$ (5)

The sequence obtained by Eq. (5) is guaranteed to converge if $f(\mathbf{x})$
is concave. However, from Lemma 1 in Subsection 3.3 we know, for non-concave
functions, SGM [21] may generate divergent sequences and never converge, since
changing too many entries of $\mathbf{x}$ at a time may increase the value of
the loss function. To overcome this difficulty, the proposed DPCD method only
changes the signs of entries whose absolute values of directional derivatives
are large (entries with principal coordinates). This yields convergence for a
wide class of loss functions (please refer to Lemma 1 and Theorem 1). Also, in
the proposed algorithm, the constraint $\Omega_{r}$ is always satisfied after
each iteration.

To be more specific, given $\mathbf{x}^{k}$ at the $k$-th iteration, we first
calculate the gradient $\nabla
f(\mathbf{x}^{k})=(\nabla_{1}f(\mathbf{x}^{k}),\nabla_{2}f(\mathbf{x}^{k}),\cdots,\nabla_{n}f(\mathbf{x}^{k}))$
for the differentiable loss function $f(\mathbf{x})$. Next, we derive some
proper thresholds $L_{1},L_{2}$ based on $\nabla f(\mathbf{x}^{k})$. When
$\nabla f$ is $L_{0}$-Lipschitz continuous on $[-1,1]^{n}$, where $L_{0}$ is
easy to calculate, we simply set

$\displaystyle L_{1}=L_{2}=L_{0}+\epsilon$ (6)

for some sufficiently small positive constant $\epsilon>0$, i.e., we consider
$(L_{0}+\epsilon)$-principal coordinates. For example, in Lemma 1, it is easy
to see that $L_{0}=1$, then we take $L_{1}=L_{2}=1+\epsilon$ for some small
$\epsilon>0$ in the algorithm. When $L_{0}$ does not exist or is difficult to
compute, we let $L_{1}$ and $L_{2}$ be the averages of the absolute values of
the positive and negative entries in the gradient $\nabla f(\mathbf{x}^{k})$,
respectively, i.e.,

$\displaystyle
L_{1}=\frac{1}{n_{1}}\sum_{\nabla_{i}f(\mathbf{x}^{k})>0}\nabla_{i}f(\mathbf{x}^{k})\,\,\text{and}\,\,L_{2}=-\frac{1}{n_{2}}\sum_{\nabla_{i}f(\mathbf{x}^{k})<0}\nabla_{i}f(\mathbf{x}^{k}),$
(7)

where $n_{1}$ and $n_{2}$ are the numbers of positive and negative entries in
$\nabla f(\mathbf{x}^{k})$, respectively. With the given thresholds $L_{1}$
and $L_{2}$, we set $S^{k+}=\\{1\leq i\leq
n:\nabla_{i}f(\mathbf{x}^{k})>\alpha_{1}L_{1},\,\mathbf{x}^{k}_{i}=1\\}$ and
$S^{k-}=\\{1\leq i\leq
n:\nabla_{i}f(\mathbf{x}^{k})<-\alpha_{2}L_{2},\,\mathbf{x}^{k}_{i}=-1\\}$ to
be the sets of $L_{1}$-principal coordinates with positive partial derivatives
and $L_{2}$-principal coordinate with negative partial derivatives,
respectively, where $\alpha_{1}$ and $\alpha_{2}$ are some parameters in
$[0.1,10]$ which will be learned depending on the tasks and datasets. If the
restriction condition is $\mathbf{x}\in\Omega_{-1}$, we update
$\mathbf{x}^{k+1}_{i}$ by solving $\min_{\mathbf{x}\in\\{\pm
1\\}^{n}}\hat{f}^{k}(\mathbf{x})$ in Eq. (5) with respect to $i\in S^{k+}\cup
S^{k-}$ (other entries of $\mathbf{x}$ are fixed) and derive:

$\displaystyle\mathbf{x}^{k+1}_{i}=-\text{sgn}(\nabla_{i}f(\mathbf{x}^{k}))=-\mathbf{x}^{k}_{i}.$
(8)

In other words, we change the sign of $\mathbf{x}^{k}_{i}$ for
$\mathbf{x}^{k}_{i}=1$ with $\alpha_{1}L_{1}$-principal coordinates, and for
$\mathbf{x}^{k}_{i}=-1$ with $\alpha_{2}L_{2}$-principal coordinates. If the
number of $1$s in $\mathbf{x}$ is required to be fixed (i.e., the restriction
condition is $\mathbf{x}\in\Omega_{r}$ for some $r\in\mathbb{N}$), we update
Eq. (8) with respect to the $m$ largest absolute values in
$\\{|\nabla_{i}f(\mathbf{x}^{k})|:i\in S^{k+}\\}$ and
$\\{|\nabla_{j}f(\mathbf{x}^{k})|:j\in S^{k-}\\}$, respectively, where
$m=\min\\{\\#S^{k+},\,\\#S^{k-}\\}$ (such procedure guarantees that
$\mathbf{x}^{k}\in\Omega_{r}$ implies $\mathbf{x}^{k+1}\in\Omega_{r}$). When
the complexity of gradient calculations is low, the updating is very fast. A
complexity analysis of the proposed DPCD is given in Subsection 4.2, which
shows that the algorithm complexity of DPCD for supervised discrete hashing is
linear and thus the algorithm runs very fast for large-scale datasets.

#### Neighborhood search.

We add an optional heuristic neighborhood search after the principal
coordinate update to avoid saddle points. In practice, we run one neighborhood
search after $T$ principal coordinate updates, where $T$ is between 10 to 20.
First, we define the concept of _$m$ -neighbors_ for a point
$\mathbf{x}\in\Omega_{r}$ where $m\in\mathbb{N}$ is some small positive
integers. When $r=-1$, the set of $m$-neighbors for $\mathbf{x}\in\Omega_{-1}$
is denoted by $N_{-1}(\mathbf{x}):=\\{\mathbf{y}\in\\{\pm
1\\}^{n}:0<d_{H}(\mathbf{y},\mathbf{x})\leq m\\}$, which is the set of points
with a Hamming distance at most $m$ from $\mathbf{x}$. When $r\geq 0$, the set
of $m$-neighbors for $\mathbf{x}\in\Omega_{r}$ is denoted by
$N(\mathbf{x}):=\\{\mathbf{y}\in\\{\pm
1\\}^{n}:0<d_{H}(\mathbf{y},\mathbf{x})\leq
2m,\,\,\sum_{i=1}^{n}\mathbf{y}_{i}=\sum_{i=1}^{n}\mathbf{x}_{i}\\}$, which is
the set of points obtained by interchanging at most $m$ pairs of entries $1$
and $-1$ in $\mathbf{x}$. For instance, when $m=1$, we have
$N_{-1}((1,-1,1,1)):=\\{(-1,-1,1,1),(1,1,1,1),(1,-1,-1,1),(1,-1,1,-1)\\}$ and
$N((1,-1,1,1)):=\\{(-1,1,1,1),(1,1,-1,1),(1,1,1,-1)\\}$.

In a neighbor search for some $\mathbf{x}^{*}\in\Omega_{-1}$ (or $\Omega_{r}$
with $r\in\mathbb{N}_{\geq 0}$), the aim is to find the point $\mathbf{y}\in
N_{-1}(\mathbf{x}^{*})\cup\\{\mathbf{x}^{*}\\}$ (or
$N(\mathbf{x}^{*})\cup\\{\mathbf{x}^{*}\\}$) with the minimal function value
$f(\mathbf{y})-f(\mathbf{x}^{*})$ (or equivalently, $f(\mathbf{y})$). In
practice, we can sample from $N_{-1}(\mathbf{x}^{*})$ (or $N(\mathbf{x}^{*})$)
instead of iterating over all points, if $n$ is large or calculating
$f(\mathbf{y})-f(\mathbf{x}^{*})$ is slow. This neighborhood search step is
helpful for finding a local minimum point. The proposed algorithm is
summarized in Algorithm 1.

Input: Loss function $f(\mathbf{x})$, code length $n$, the restriction
$\Omega_{r}$ where $r=-1\,\text{or}\,r\in\mathbb{N}_{\geq 0}$, parameters
$\alpha_{1}$, $\alpha_{2}$.

Output: Binary codes $\mathbf{x}^{*}$.

Initialize $\mathbf{x}^{*}$ by the sign of some random vector according to
$\Omega_{r}$; $\mathbf{x}^{1}=\mathbf{x}^{*}$ and $k=1$;

while  _not converge or not reach maximum iterations_ do

Calculate the gradient $\nabla
f(\mathbf{x}^{k})=(\nabla_{1}f(\mathbf{x}^{k}),\nabla_{2}f(\mathbf{x}^{k}),\cdots,\nabla_{n}f(\mathbf{x}^{k}))$;

Derive proper thresholds $L_{1},L_{2}$ by Eq. (6) or Eq. (7);

Build sets
$S^{k+}=\\{i:\nabla_{i}f(\mathbf{x}^{k})>\alpha_{1}L_{1},\,\mathbf{x}^{k}_{i}=1\\}$
and
$S^{k-}=\\{i:\nabla_{i}f(\mathbf{x}^{k})<-\alpha_{2}L_{2},\,\mathbf{x}^{k}_{i}=-1\\}$;

if  _the restriction condition is $\mathbf{x}\in\Omega_{-1}$ (i.e., $r=-1$)_
then

Update
$\mathbf{x}^{k+1}_{i}=-\text{sgn}(\nabla_{i}f(\mathbf{x}^{k}))=-\mathbf{x}^{k}_{i}$
for $i\in S^{k+}\cup S^{k-}$;

else

Sort $\\{|\nabla_{i}f(\mathbf{x}^{k})|:i\in S^{k+}\\}$ and
$\\{|\nabla_{j}f(\mathbf{x}^{k})|:j\in S^{k-}\\}$ in descending order as
$|\nabla_{i_{1}}f(\mathbf{x}^{k})|\geq|\nabla_{i_{2}}f(\mathbf{x}^{k})|\geq|\nabla_{i_{3}}f(\mathbf{x}^{k})|\geq\cdots$,
and
$|\nabla_{j_{1}}f(\mathbf{x}^{k})|\geq|\nabla_{j_{2}}f(\mathbf{x}^{k})|\geq|\nabla_{j_{3}}f(\mathbf{x}^{k})|\geq\cdots$,
respectively;

Update $\mathbf{x}^{k+1}_{i_{l}}=-\mathbf{x}^{k}_{i_{l}}$ and
$\mathbf{x}^{k+1}_{j_{l}}=-\mathbf{x}^{k}_{j_{l}}$ for $1\leq
l\leq\min\\{\\#S^{k+},\,\\#S^{k-}\\}$;

(Optional) Neighborhood search for $\mathbf{x}^{k+1}$;

$k=k+1$;

Return $\mathbf{x}^{*}=\mathbf{x}^{k+1}$.

Algorithm 1 Discrete Principal Coordinate Descent (DPCD)

###  3.3 Convergence comparison: DPCD vs. SGM

One of the differences between our DPCD and SGM [21] is the choice of
$\mathbf{x}_{i}^{k}$ that should be updated by Eq. (8) at the $k$-th
iteration. SGM updates Eq. (8) for each $i$, while the proposed DPCD only
changes the sign of $\mathbf{x}_{i}^{k}$ when the coordinate indexes is
$L$-principal for some $L$. This difference is crucial to the convergences of
the algorithms. More specifically, SGM can only guarantee convergence for the
minimization of concave functions, while DPCD converges in finite steps for
any functions with Lipschitz continuous gradients. We show the superiority of
the proposed DPCD by the following example. The case $n=2$ is presented in
Figure 1.

Figure 1: The optimization routes of DPCD and SGM on a 2-d example. The black
circles are the contour lines of the objective function. Our method DPCD can
converge to the optimum from any initialization points in only one step, while
SGM always oscillates between two points. In fact, each $\beta_{i}$ can be in
the open interval (-1,1), and the optimum is decided by the signs of all
$\beta_{i}$.

######  Lemma 1.

Consider the problem

$\min_{\mathbf{x}\in\\{\pm
1\\}^{n}}f(\mathbf{x}):=\min_{(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n})\in\\{\pm
1\\}^{n}}~{}\frac{1}{2}\sum_{i=1}^{n}(\mathbf{x}_{i}+{\beta_{i}})^{2},$

where $0<\beta_{i}<1$. Let $\beta=\min_{i}\beta_{i}$. It holds that: (1) SGM
generates a divergent sequence for any initial point; (2) With parameters
$\alpha_{1}=\alpha_{2}=1$ and $0<\epsilon<\beta$, the proposed DPCD method
always converges to the optimal solution.

###### Proof.

First, we apply the method SGM to this problem. The gradient of
$f(\mathbf{x})$ is straightforward: $\nabla
f(\mathbf{x})=(\mathbf{x}_{1}+{\beta_{1}},\mathbf{x}_{2}+{\beta_{2}},\cdots,\mathbf{x}_{n}+{\beta_{n}})$.
Since $0<\beta_{i}<1$, we have $\mathbf{x}=\text{sgn}(\nabla f(\mathbf{x}))$
for any $\mathbf{x}\in\\{\pm 1\\}^{n}$. Then, by the updating rule fo SGM,
$\mathbf{x}^{k+1}=-\text{sgn}(\nabla f(\mathbf{x}^{k}))=-\mathbf{x}^{k}$.
Therefore, starting from any point $\mathbf{x}\in\\{\pm 1\\}^{n}$, SGM always
generates a divergent sequence
$\mathbf{x},-\mathbf{x},\mathbf{x},-\mathbf{x},\mathbf{x},-\mathbf{x},\cdots$.

On the other hand, since $\|\nabla f(\mathbf{y})-\nabla
f(\mathbf{z})\|=\|\mathbf{y}-\mathbf{z}\|$ for any
$\mathbf{y},\mathbf{z}\in\\{\pm 1\\}^{n}$, $f(\mathbf{x})$ is a function with
$1$-Lipschitz continuous gradient. Therefore, in DPCD, we set
$1<L_{1}=L_{2}=1+\epsilon<1+\beta$. Then by definition,

$\displaystyle S^{k+}$ $\displaystyle=\\{1\leq i\leq
n:\nabla_{i}f(\mathbf{x}^{k})>1+\epsilon,\,\mathbf{x}^{k}_{i}=1\\}$
$\displaystyle=\\{1\leq i\leq
n:\mathbf{x}^{k}_{i}+\beta_{i}>1+\epsilon,\,\mathbf{x}^{k}_{i}=1\\}$
$\displaystyle=\\{1\leq i\leq n:\,\mathbf{x}^{k}_{i}=1\\},$

and

$\displaystyle S^{k-}$ $\displaystyle=\\{1\leq i\leq
n:\nabla_{i}f(\mathbf{x}^{k})<-1-\epsilon,\,\mathbf{x}^{k}_{i}=-1\\}$
$\displaystyle=\\{1\leq i\leq
n:\mathbf{x}^{k}_{i}+\beta_{i}<-1-\epsilon,\,\mathbf{x}^{k}_{i}=-1\\}$
$\displaystyle=\emptyset.$

This implies that

$\displaystyle\mathbf{x}^{k+1}_{i}=\begin{cases}-{\text{sgn}}(\nabla_{i}f(\mathbf{x}^{k}))=-\mathbf{x}^{k}_{i}=-1,\,\text{if}\,\,i\in
S^{k+}=S^{k+}\cup S^{k-}\\\ \mathbf{x}^{k}_{i}=-1,\,\text{if}\,\,i\notin
S^{k+}.\end{cases}$ (9)

Then, we have $\mathbf{x}^{k+1}_{i}=-1$ for any $k\geq 2$ and $i$, and thus
$\mathbf{x}^{2}=\mathbf{x}^{3}=\cdots=\mathbf{x}^{k}=(-1,-1,\cdots,-1)$ for
$k\geq 2$. It is easy to check that $(-1,-1,\cdots,-1)$ is indeed the optimal
point for the problem in Lemma $1$. Therefore, the proposed DPCD converges to
the optimal solution in only one updating step for any initial point. ∎

###  3.4 Theoretical convergence results

For simplicity, we ignore the neighborhood search part in the convergence
analysis. Actually, the neighborhood search does not have any influence on the
convergence since it always generates some binary vectors with non-increasing
values of the loss function. Thus, we can only focus on the principal
coordinate update part of the proposed DPCD method.

We derive the following convergence results. When we say the algorithm
converges in $T$ steps, we mean that the binary vector $\mathbf{x}^{T+1}$
obtained in $(T+1)$-th iteration equals $\mathbf{x}^{T}$ in $T$-th iteration
(then the algorithm can stop here). Moreover, it is easy to see that our
algorithm can converge to some local optimum with the help of neighborhood
search (without neighborhood search, it only converges to some fixed binary
vectors).

######  Theorem 1.

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a differentiable function such
that $\nabla f$ is $L_{0}$-Lipschitz continuous on $[-1,1]^{n}$. Setting the
thresholds $L_{1}=L_{2}=L_{0}+\epsilon$ where $\epsilon>0$, and ignoring the
neighborhood search, Algorithm 1 always converges in
$\frac{f_{max}-f_{min}}{2\epsilon}$ steps at most, where
$f_{max}:=\max_{\mathbf{x}\in\\{\pm 1\\}^{n}}f(\mathbf{x})$ and
$f_{min}:=\min_{\mathbf{x}\in\\{\pm 1\\}^{n}}f(\mathbf{x})$.

###### Proof.

Let $L=L_{1}=L_{2}=L_{0}+\epsilon$. Since $\nabla f$ is $L_{0}$-Lipschitz, we
have for any $\mathbf{x},\mathbf{y}\in[-1,1]^{n}$,

$\|\nabla f(\mathbf{x})-\nabla f(\mathbf{y})\|\leq
L_{0}\|\mathbf{x}-\mathbf{y}\|.$

Then, by Cauchy-Schwarz inequality,

$(\nabla f(\mathbf{x})-\nabla
f(\mathbf{y})^{\top}(\mathbf{x}-\mathbf{y}))\leq\|\nabla f(\mathbf{x})-\nabla
f(\mathbf{y})\|\cdot\|\mathbf{x}-\mathbf{y}\|\leq
L_{0}\rVert\mathbf{x}-\mathbf{y}\rVert^{2}.$

By the gradient monotonicity equivalence of convexity [6, Page 40], this
yields that $g(x)=\frac{L_{0}}{2}\|\mathbf{x}\|^{2}-f(\mathbf{x})$ is a convex
function on $[-1,1]^{n}$. Therefore, due to the first-order equivalence of
convexity [8, Page 69],

$\frac{L_{0}}{2}\|\mathbf{y}\|^{2}-f(\mathbf{y})\geq\frac{L_{0}}{2}\|\mathbf{x}\|^{2}-f(\mathbf{x})+(L_{0}\mathbf{x}-\nabla
f(\mathbf{x}))^{\top}(\mathbf{y}-\mathbf{x}),$

which implies that

$f(\mathbf{y})\leq f(\mathbf{x})+\nabla
f(\mathbf{x})^{\top}(\mathbf{y}-\mathbf{x})+\frac{L_{0}}{2}\lVert\mathbf{y}-\mathbf{x}\rVert^{2},~{}\forall\mathbf{x},\mathbf{y}\in[-1,1]^{n}.$

Let $\mathbf{x}=\mathbf{x}^{k}$ and $\mathbf{y}=\mathbf{x}^{k+1}$ in the above
inequality, where $\mathbf{x}^{k+1}$ is determined by $\mathbf{x}^{k}$ and (8)
in the paper. We have

$\displaystyle f(\mathbf{x}^{k+1})\leq f(\mathbf{x}^{k})+\nabla
f(\mathbf{x}^{k})^{\top}(\mathbf{x}^{k+1}-\mathbf{x}^{k})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}-\mathbf{x}^{k}\rVert^{2}.$
(10)

If $\mathbf{x}^{k+1}=\mathbf{x}^{k}$, we obtain
$\mathbf{x}^{j+1}=\mathbf{x}^{j}$ for each $j\geq k$ due to the updating rule.

If $\mathbf{x}^{k+1}_{i}\neq\mathbf{x}^{k}_{i}$ for some $1\leq i\leq n$, by
the updating rule we have,
$\mathbf{x}^{k+1}_{i}=-{\text{sgn}}(\nabla_{i}f(\mathbf{x}^{k}))=-\mathbf{x}^{k}_{i}$
and $\|\nabla_{i}f(\mathbf{x}^{k})\|>L$. Then,

$\nabla_{i}f(\mathbf{x}^{k})(\mathbf{x}^{k+1}_{i}-\mathbf{x}^{k}_{i})=-2\,\nabla_{i}f(\mathbf{x}^{k})\,{\text{sgn}}(\nabla_{i}f(\mathbf{x}^{k}))=-2\,\|\nabla_{i}f(\mathbf{x}^{k})\|<-2L$

and
$\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{i}-\mathbf{x}^{k}_{i}\rVert^{2}=\frac{L_{0}}{2}\cdot
4=2L_{0}.$ Therefore,

$\displaystyle\nabla_{i}f(\mathbf{x}^{k})^{\top}(\mathbf{x}^{k+1}_{i}-\mathbf{x}^{k}_{i})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{i}-\mathbf{x}^{k}_{i}\rVert^{2}<-2(L-L_{0}).$
(11)

If $\mathbf{x}^{k+1}\neq\mathbf{x}^{k}$, we know
$\\{j:\,\mathbf{x}^{k+1}_{j}\neq\mathbf{x}^{k}_{j}\\}$ is not empty. Then, by
Eq. (10) and Eq. (11) we have,

$\displaystyle f(\mathbf{x}^{k+1})-f(\mathbf{x}^{k})\leq$ $\displaystyle\nabla
f(\mathbf{x}^{k})^{\top}(\mathbf{x}^{k+1}-\mathbf{x}^{k})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}-\mathbf{x}^{k}\rVert^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}\left(\nabla_{j}f(\mathbf{x}^{k})(\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j}\rVert^{2}\right)$
$\displaystyle=$
$\displaystyle\sum_{j:\,\mathbf{x}^{k+1}_{j}=\mathbf{x}^{k}_{j}}\left(\nabla_{j}f(\mathbf{x}^{k})(\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j}\rVert^{2}\right)$
$\displaystyle+\sum_{j:\,\mathbf{x}^{k+1}_{j}\neq\mathbf{x}^{k}_{j}}\left(\nabla_{j}f(\mathbf{x}^{k})(\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j}\rVert^{2}\right)$
$\displaystyle=$ $\displaystyle
0+\sum_{j:\,\mathbf{x}^{k+1}_{j}\neq\mathbf{x}^{k}_{j}}\left(\nabla_{j}f(\mathbf{x}^{k})^{\top}(\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j})+\frac{L_{0}}{2}\lVert\mathbf{x}^{k+1}_{j}-\mathbf{x}^{k}_{j}\rVert^{2}\right)$
$\displaystyle<$ $\displaystyle-2(L-L_{0})$ $\displaystyle=$
$\displaystyle-2\epsilon.$ (12)

This means that, in each updating step, the function value decreases at least
$2\epsilon$, which completes the proof since the feasible set of Problem (2)
in the main paper is finite. ∎

The above theorem is very versatile since most loss functions in practice are
$L_{0}$-Lipschitz continuous on $[-1,1]^{n}$ for some
$L_{0}\in\mathbb{R}^{+}$. Furthermore, since $n$ is finite, $\\{\pm 1\\}^{n}$
is a finite set, thus the above $f_{max}$ and $f_{min}$ always exist and are
finite. When $f(\mathbf{x})$ is quadratic, we obtain the following direct
corollary.

######  Corollary 1.

Let
$f(\mathbf{x})=\mathbf{x}^{\top}\mathbf{A}\mathbf{x}+\mathbf{c}^{\top}\mathbf{x}+d$
be some quadratic function where $\mathbf{A}\in\mathbb{R}^{n\times n}$,
$\mathbf{x},\mathbf{c}\in\mathbb{R}^{n}$, and $d\in\mathbb{R}$, in Theorem 1.
Then Algorithm 1 always converges in
$\frac{\|\mathbf{A}\|_{1}+\|\mathbf{c}\|_{1}}{\epsilon}$ steps at most.

###### Proof.

Since $\mathbf{x}\in\\{\pm 1\\}^{n}$, we have
$-\|\mathbf{A}\|_{1}-\|\mathbf{c}\|_{1}+d\leq
f(\mathbf{x})\leq\|\mathbf{A}\|_{1}+\|\mathbf{c}\|_{1}+d$. Then
$f_{max}-f_{min}\leq 2(\|\mathbf{A}\|_{1}+\|\mathbf{c}\|_{1})$. By Theorem $1$
we complete the proof. ∎

##  4 Experiments

In this section, we compare the proposed DPCD algorithm with several state-of-
the-art methods on two binary optimization tasks: dense subgraph discovery and
binary hashing. All codes are implemented in MATLAB using a workstation with
an Intel 8-core 2.6GHz CPU and 32GB RAM.

###  4.1 Dense subgraph discovery

Table 1: Statistics for the graphs used in the dense subgraph discovery
experiments.

Graph | # Nodes | # Arcs | # Arcs/# Nodes
---|---|---|---
uk-2007-05 | 100000 | 3050615 | 30.506
dblp-2010 | 326186 | 1615400 | 4.952
eswiki-2013 | 972933 | 23041488 | 23.683
hollywood-2009 | 1139905 | 113891327 | 99.913

#### Optimization problem.

Dense subgraph discovery [29, 30, 43, 46] has many applications in graph
mining, such as real-time story identification [3], finding correlated genes
[47] and graph visualization [1]. Let $\sl{G}$ be a given undirected weighted
graph with $n$ nodes, and $k$ be a given positive integer such that $1\leq
k\leq n$. The aim is to find the maximum density subgraph (the subgraph with
the maximal sum of edge weights) with cardinality $k$. Let
$\mathbf{W}\in\mathbb{R}^{n\times n}$ be the symmetric adjacency matrix of the
graph $G$, where $\mathbf{W}_{ij}$ denotes the weight of the edge between
vertices $i$ and $j$. Then the problem can be formulated as the following
optimization problem

$\max_{\mathbf{x}\in\\{0,1\\}^{n}}\,\,\mathbf{x}^{\top}\mathbf{W}\mathbf{x},\quad\,\,\,\text{s.t.}~{}~{}\mathbf{x}^{\top}\mathbf{1}=k.$
(13)

Note that, in this case, the variables are $0$ or $1$ instead of $-1$ or $1$.
In order to translate the problem to the form in Problem (2), the substitution
$\mathbf{x}=\frac{1}{2}(\mathbf{y}+\mathbf{1})$ is adopted. Therefore, Problem
(13) is equivalent to:

$\min_{\mathbf{y}\in\\{-1,1\\}^{n}}\,\,-\mathbf{y}^{\top}\mathbf{W}\mathbf{y}-2\mathbf{y}^{\top}\mathbf{W}\mathbf{1}-\mathbf{1}^{\top}\mathbf{W}\mathbf{1},~{}~{}~{}\text{s.t.}~{}~{}\mathbf{y}^{\top}\mathbf{1}=2k-n.$
(14)

In this way, the above problem can be approximately solved using our Algorithm
1.

(a) uk-2007-05

(b) dblp-2010

(c) eswiki-2013

(d) hollywood-2009

Figure 2: Experimental results for dense subgraph discovery with
$k\in\\{200,400,800,1600,3200,6400\\}$. The metric is
$\mathbf{x}^{\top}\mathbf{W}\mathbf{x}/k$.

#### Graph datasets.

The experiments for dense subgraph discovery are conducted on four large-scale
graphs uk-2007-05, dblp-2010, eswiki-2013, and hollywood-2009 from The
Laboratory for Web Algorithmics111http://law.di.unimi.it/datasets.php. Table 1
gives a brief description of each graph. For example, hollywood-2009 contains
roughly 1.13 million nodes and 113 million edges.

Table 2: The CPU time comparison (seconds) of dense subgraph discovery on the
four graphs with $k=1600$.

Method | uk-2007-05 | dblp-2010 | eswiki-2013 | hollywood-2009
---|---|---|---|---
DPCD | 1.77 | 6.19 | 14.17 | 21.80
LP | 5.02 | 7.27 | 58.96 | 115.38
RAVI | 2.28 | 6.40 | 104.94 | 77.95
L2box-ADMM | 35.84 | 74.54 | 335.83 | 610.72
MPEC-EPM | 45.21 | 158.11 | 1454.07 | 2114.20
MPEC-ADM | 36.76 | 81.64 | 421.35 | 720.66

#### DPCD vs. state-of-the-art methods.

The proposed DPCD is compared with the methods LP [13], RAVI [29], L2box-ADMM
[41], MPEC-EPM and MPEC-ADM [43], using the same objective function (14). The
cardinality $k$ is in the set $\\{200,400,800,1600,3200,6400\\}$. For the DPCD
method, we run one $5$-neighborhood search ($m=5$ in the neighborhood search
setting) after $10$ principal coordinate updates, and set the maximum
iteration number for the principal coordinate update part to $100$. $L_{1}$
and $L_{2}$ are updated by Eq. (7). We tune the parameters $\alpha_{1}$ and
$\alpha_{2}$ from $\\{0.1,0.2,0.3,\ldots,0.9,1,2,3,\ldots,10\\}$ by cross-
validation according to datasets. For other methods, we adopt the
implementations and parameters suggested by the authors. The experimental
results in Figure 2 are reported in terms of
$\mathbf{x}^{\top}\mathbf{W}\mathbf{x}/k$, which is the density of the
subgraph with $k$ vertices. We can see that DPCD finds a denser subgraph than
all compared methods in each case. Among the other methods, MPEC-EPM
consistently outperforms LP in all the experiments. RAVI generally leads to
solutions with low density. Furthermore, we provide CPU time comparisons for
the six methods on the four graphs. From Table 2 we can see that, the proposed
DPCD achieves the fastest runtime on all graphs. This is due to the fast
updates in Algorithm 1 (the main complexity is to calculate the gradient of
the loss function, which can be done quickly). Also, the efficiency of our
method becomes more obvious as the graph size increases. LP and RAVI are
faster than the other three methods, since MPEC-EPM needs to run the LP
procedure multiple times, and L2-box ADMM and MPEC-ADM usually need more
iterations to converge.

Table 3: Comparison between DPCD with and without neighborhood search.

Method (Subgraph) | uk-2007-05 | dblp-2010 | eswiki-2013 | hollywood-2009
---|---|---|---|---
DPCD (loss function) | 96.573 | 45.238 | 158.340 | 1944.921
DPCD-0 (loss function) | 93.5585 | 44.2512 | 147.815 | 1866.240
DPCD (run time) | 1.7689 | 6.1092 | 14.1745 | 21.8038
DPCD-0 (run time) | 1.1522 | 5.5186 | 9.9596 | 12.8730

#### With and without the neighborhood search.

We conduct the comparison of the proposed method with its variant without
neighborhood search technique, on the task of dense subgraph discovery with
cardinality $1600$ and in terms of loss functions and run time. DPCD-0 refers
to our algorithm without neighborhood search. From Table 3 we can see that
DPCD-0 is slightly faster, while DPCD usually achieves better results:

###  4.2 Binary hashing

Table 4: Evaluation of DPCD and five general binary optimization methods with
the same supervised loss function. The CIFAR-10 dataset is adopted. Results
are reported in terms of MAP, Precision@500 and the training time.

Method | MAP | Precision@500 | Training time (seconds)
---|---|---|---
$32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits
DPCD | 0.7019 | 0.7088 | 0.7126 | 0.6337 | 0.6353 | 0.6370 | 3.76 | 6.01 | 9.39
DCC | 0.5941 | 0.6193 | 0.6314 | 0.5486 | 0.5766 | 0.5894 | 11.18 | 36.03 | 158.87
SGM | 0.6856 | 0.6986 | 0.7013 | 0.6177 | 0.6308 | 0.6360 | 7.87 | 10.83 | 16.30
LP | 0.5237 | 0.5468 | 0.5459 | 0.4704 | 0.4972 | 0.4866 | 4.52 | 7.96 | 13.64
L2box-ADMM | 0.6399 | 0.6724 | 0.6830 | 0.5929 | 0.6095 | 0.6162 | 43.07 | 90.10 | 200.94
MPEC-EPM | 0.5823 | 0.6253 | 0.6276 | 0.5385 | 0.5738 | 0.5790 | 36.36 | 124.54 | 260.22

Binary hashing aims to encode high-dimensional data points, such as images and
videos, into compact binary hash codes such that the similarities between the
original data points and hash codes are preserved. This can be used to provide
a constant or sub-linear search time and reduce the storage cost dramatically
for such data points. The efficiency and effectiveness of binary hashing make
it a popular technique in machine learning, information retrieval and computer
vision [11, 14, 21, 33, 35]. In a typical binary hashing task,
$\mathbf{X}=(\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n})^{\top}\in\mathbb{R}^{n\times
d}$ denotes a matrix of the original data points, where
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ is the $i$-th sample point, $n$ is the
number of samples, and $d$ is the dimension of each sample. In the supervised
setting, we let $\mathbf{Y}\in\mathbb{R}^{n\times c}$ be the label matrix,
i.e., $\mathbf{Y}_{i,j}=1$ if $\mathbf{x}_{i}$ belongs to the $j$-th class and
$0$ otherwise. The aim is to map $\mathbf{X}$ to some
$\mathbf{B}=(\mathbf{b}_{1},\mathbf{b}_{2},\ldots,\mathbf{b}_{n})^{\top}\in\\{\pm
1\\}^{n\times r}$, i.e., map each $x_{i}$ to a binary code
$\mathbf{b}_{i}\in\\{\pm 1\\}^{r}$ for some small integer $r$ and preserve
some similarities between the original data points in $\mathbf{X}$ and hash
codes in $\mathbf{B}$.

#### Image datasets.

Three large-scale image datasets,
CIFAR-10222http://www.cs.toronto.edu/kriz/cifar.html.,
ImageNet333http://www.image-net.org/., and NUS-
WIDE444http://lms.comp.nus.edu.sg/research/NUS-WIDE.htm., are used in the
binary hashing experiments. CIFAR-10 has 60k images, which are divided into 10
classes with 6k images each. We use a 384-dimensional GIST feature vector [26]
to represent each image. 59k images are selected as the training set and the
test set contains the remaining 1k images. A subset of the ImageNet, ILSVRC
2012, contains about 1.2 million images with 1k categories. As in [14, 32], we
use 4096-dimensional deep feature vectors for each image, take 127K training
images from the 100 largest classes, and 50K images from the validation set as
the test set. NUS-WIDE contains about 270K images with 81 labels. The images
may have multiple labels. The 500-dimensional Bag-of-Words features are used
here [9]. We adopt the 21 most frequent labels with the corresponding 193K
images. For each label, 100 images are randomly selected as the test set and
the remaining as the training set.

#### DPCD vs. general binary optimization methods.

To illustrate the efficiency and effectiveness of our algorithm, we compare
DPCD with several general state-of-the-art binary optimization methods DCC
[32], SGM [21], LP [13], L2box-ADMM [41], and MPEC-EPM [43], on the dataset
CIFAR-10. Various loss functions are designed for binary hashing (for
examples, see Tables 5 and 6). For fair comparison, we adopt the widely used
supervised objective function [11, 32]

$\displaystyle
f(\mathbf{B},\mathbf{W})=\frac{1}{2}\|\mathbf{Y}-\mathbf{B}\mathbf{W}\|_{2}^{2}+\frac{\delta}{2}\|\mathbf{W}\|_{2}^{2}$
(15)

for each method. Thus the optimization problem becomes:

$\displaystyle\min_{\mathbf{B},\mathbf{W}}~{}\frac{1}{2}\|\mathbf{Y}-\mathbf{B}\mathbf{W}\|_{2}^{2}+\frac{\delta}{2}\|\mathbf{W}\|_{2}^{2}~{}~{}~{}~{}~{}~{}\text{s.t.}~{}~{}~{}\mathbf{B}\in\\{\pm
1\\}^{n\times r},\mathbf{W}\in\mathbb{R}^{r\times c},$ (16)

where $\beta$ is a regularization parameter, and
$\mathbf{W}\in\mathbb{R}^{r\times c}$ is the projection matrix (see [32])
which will be learned jointly with $\mathbf{B}$. The whole optimization runs
iteratively over $\mathbf{B}$ and $\mathbf{W}$. When $\mathbf{W}$ is fixed, we
apply DPCD algorithm to $\mathbf{B}$. The key step is to calculate the
gradient of $f(\mathbf{B},\mathbf{W})$ as:

$\displaystyle\nabla_{\mathbf{B}}f(\mathbf{B},\mathbf{W})=(\mathbf{B}\mathbf{W}-\mathbf{Y})\mathbf{W}^{\top}.$
(17)

Then $L_{1},L_{2}$ can be obtained by Eq. (7). After deriving $S^{k+}$ and
$S^{k-}$, we update $\mathbf{B}$ by Eq. (8). When $\mathbf{B}$ is fixed,
$\mathbf{W}$ can be updated by

$\displaystyle\mathbf{W}={\arg}\min_{\mathbf{W}^{*}\in\mathbb{R}^{r\times
c}}f(\mathbf{B},\mathbf{W}^{*})=(\mathbf{B}^{\top}\mathbf{B}+\beta\mathbf{I}_{r})^{-1}\mathbf{B}^{\top}\mathbf{Y}.$
(18)

Finally, we adopt the linear hash function
$h(\mathbf{X})=\text{sgn}(\mathbf{X}\mathbf{P})$ to encode $\mathbf{X}$ onto
binary codes, where $\mathbf{P}\in\mathbb{R}^{d\times r}$ can be derived by:

$\displaystyle\mathbf{P}={\arg}\min_{\mathbf{P}^{*}\in\mathbb{R}^{d\times
r}}\|\mathbf{X}\mathbf{P}^{*}-\mathbf{B}\|^{2}=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{B}.\hskip
12.0pt$ (19)

During the test phase, for a query item, first we use the above linear hash
function to derive the hash code, then adopt nearest neighbor search under
Hamming distance to find its similar items. For the DPCD method, we run one
$5$-neighborhood search after $10$ principal coordinate updates, and tune the
parameters $\alpha_{1}$ and $\alpha_{2}$ from
$\\{0.1,0.2,0.3,\ldots,0.9,1,2,3,\ldots,10\\}$ by cross-validation according
to datasets and binary code lengths, and set the maximum iteration number for
$\mathbf{B}$ to $20$ each time when $\mathbf{W}$ is fixed. We run at most five
iterations for updating $\mathbf{B}$ and $\mathbf{W}$ iteratively. For other
methods, we adopt the implementations and parameters suggested by the authors.
Ground truths are defined by the label information from the datasets. The
experimental results are reported in terms of mean average precision (MAP),
Precision@500 (Precision@500 refers to the ratio of the number of retrieved
true positive items among 500 nearest neighbors to 500.) and training time
efficiency (we ignore the comparison of test time here since the test parts
are similar for each algorithm). The experiments are conducted on the CIFAR-10
dataset. From Table 4 we can see that the proposed DPCD outperforms all other
methods in MAP and Precision@500. For instance, on CIFAR-10 with 96 bits, DPCD
outperforms DCC by 8.1$\%$, and MPEC-EPM by 8.5$\%$ in terms of MAP. It is
clear that increasing the number of bits yields better performance for all
methods. Also, the training time for the proposed DPCD method is always faster
than other compared methods. For example, DPCD runs 20 times faster than MPEC-
EPM on $64$ bits, which verifies that one major advantage of our method is the
fast optimization process.

#### Algorithm complexity analysis.

Now we discuss the complexity of the above supervised DPCD algorithm. The
calculation of the gradient $\nabla_{\mathbf{B}}f(\mathbf{B}^{k},\mathbf{W})$
is obtained by Eq. (17), thus the complexity is $O(nrc)$. Then $L_{1}$ and
$L_{2}$ are derived by Eq. (7), which has complexity $O(nr)$ since there are
$nr$ additions in Eq. (7). Furthermore, $S^{k+}$ and $S^{k-}$ can be
determined by running through all $\nabla_{ij}f(\mathbf{B}^{k})$ where $1\leq
i\leq n,\,1\leq j\leq r$, which also has complexity $O(nr)$. Similarly, the
update for $\mathbf{B}$ has complexity $O(nr)$. Let $T$ be the maximum
iteration number during the $\mathbf{B}$ updating step (the number of
principal coordinate update parts of DPCD). Then, the total time complexity of
updating $\mathbf{B}$ is $O(Tnrc)$. When $\mathbf{B}$ is given, the complexity
for updating $\mathbf{W}$ in Eq. (18) is $O(ncr+r^{3}+cr^{2})$. The complexity
for calculating the matrix $\mathbf{P}$ in Eq. (19) is $O(d^{3}+nd^{2}+ndr)$.
Suppose that there are at most $t$ iterations for updating $\mathbf{B}$ and
$\mathbf{W}$ iteratively. Since $r,c\ll n$, the total complexity for
supervised DPCD is $O\left(tTncr+d^{3}+nd^{2}+ndr\right)$, which is a linear
function of $n$.

#### DPCD vs. specific binary hashing methods.

The proposed method is not only fast, but also very versatile, and can thus be
used to handle many different loss functions. To demonstrate this, we apply
DPCD to several widely used loss functions, and compare the results with the
corresponding hashing methods that were designed for each specific loss
function. The optimization process is similar to the supervised case. Table 5
illustrates a comparison of DPCD with three unsupervised methods, SADH-L [33],
ARE [14], and ITQ [10], on the ImageNet dataset, while Table 6 shows a
comparison with supervised hashing methods, SDH [32], FSDH [11], and FastHash
[18], on NUS-WIDE, using their specific loss functions. The proposed DPCD
shows increased performance over the original methods and achieves higher MAP
and Precision@500 in most cases, especially for the supervised loss functions.
In terms of the training time, DPCD outperforms all methods except FSDH. The
proposed method can significantly decrease the training time for unsupervised
loss functions due to the fast updating of the binary codes $\mathbf{B}$.
Finally, we conclude that DPCD is a fast and effective optimization method for
large-scale image retrieval tasks.

Table 5: Evaluation of the proposed DPCD method and three unsupervised
methods. The ImageNet dataset is adopted. Results are reported in terms of
MAP, Precision@500 and training time.

Method | Loss Function | MAP | Precision@500 | Training time (s)
---|---|---|---|---
$32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits
SADH-L | $\min_{\mathbf{B}}\mathbf{Tr}(\mathbf{B}^{\top}\mathbf{L}\mathbf{B})$ | 0.2448 | 0.3104 | 0.3294 | 0.3692 | 0.4588 | 0.4835 | 121.90 | 384.01 | 738.33
DPCD | 0.2612 | 0.3085 | 0.3441 | 0.3945 | 0.4508 | 0.4992 | 19.03 | 33.50 | 48.29
ARE | $\min_{\mathbf{B}}\|\mathbf{B}\mathbf{B}^{\top}-r\mathbf{X}\mathbf{X}^{\top}\|^{2}$ | 0.2509 | 0.2997 | 0.3276 | 0.3626 | 0.4478 | 0.4724 | 244.46 | 287.19 | 332.95
DPCD | 0.2808 | 0.3316 | 0.3607 | 0.3970 | 0.4856 | 0.5112 | 6.59 | 8.37 | 11.24
ITQ | $\min_{\mathbf{B},\mathbf{R}}\|\mathbf{B}-\mathbf{X}\mathbf{W}\mathbf{R}\|^{2}$ | 0.3209 | 0.4075 | 0.4388 | 0.4269 | 0.5208 | 0.5581 | 26.01 | 32.97 | 33.39
DPCD | 0.3093 | 0.4155 | 0.4467 | 0.4148 | 0.5283 | 0.5654 | 3.35 | 7.23 | 8.02

Table 6: Evaluation of the proposed DPCD method and three supervised methods.
The NUS-WIDE dataset is adopted. Results are reported in terms of MAP,
Precision@500 and training time.

Method | Loss Function | MAP | Precision@500 | Training time (s)
---|---|---|---|---
$32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits |  $32$ bits |  $64$ bits |  $96$ bits
SDH | $\min_{\mathbf{B},\mathbf{W}}\|\mathbf{Y}-\mathbf{B}\mathbf{W}\|^{2}+{\delta}\|\mathbf{W}\|^{2}$ | 0.5716 | 0.5827 | 0.5920 | 0.6041 | 0.6056 | 0.6189 | 24.53 | 107.50 | 486.95
DPCD | 0.6124 | 0.6230 | 0.6392 | 0.6287 | 0.6357 | 0.6502 | 6.91 | 16.39 | 22.53
FSDH | $\min_{\mathbf{B},\mathbf{W}}\|\mathbf{B}-\mathbf{Y}\mathbf{W}\|_{2}^{2}+{\delta}\|\mathbf{W}\|_{2}^{2}$ | 0.5690 | 0.5639 | 0.5676 | 0.5918 | 0.5860 | 0.5988 | 2.76 | 3.07 | 5.95
DPCD | 0.6159 | 0.6170 | 0.6253 | 0.6260 | 0.6297 | 0.6335 | 5.42 | 8.96 | 10.99
FastHash | $\min_{\mathbf{B}}\|\mathbf{B}\mathbf{B}^{\top}-r\mathbf{Y}\|^{2}$ | 0.5174 | 0.5398 | 0.5403 | 0.5867 | 0.6014 | 0.6180 | 1381.75 | 4668.13 | 10605.82
DPCD | 0.5621 | 0.5579 | 0.5767 | 0.5991 | 0.6145 | 0.6226 | 6.42 | 7.88 | 11.79

##  5 Conclusion and future work

This paper presents a novel fast optimization method, called Discrete
Principal Coordinate Descent (DPCD), to approximately solve binary
optimization problems with/without restrictions on the numbers of $1$s and
$-1$s in the variables. We derive several theoretical results on the
convergence of the proposed algorithm. Experiments on dense subgraph discovery
and binary hashing demonstrate that our method generally outperforms state-of-
the-art methods in terms of both solution quality and optimization efficiency.

In the future, we plan to extend our algorithm to a more general framework.
Our methods can be seen as a discrete version of the normal gradient descent
methods. Since the gradient descent has several useful variants such as
momentum, Adam and Adagrad methods [28, 31], it would be possible to combine
our methods with these gradient-based methods and propose some discrete
versions of them. We would like to explore more on this direction in the
future.

## References

  * [1] J Ignacio Alvarez-Hamelin, Luca Dall’Asta, Alain Barrat, and Alessandro Vespignani.  Large scale networks fingerprinting and visualization using the k-core decomposition.  In Advances in Neural Information Processing Systems (NIPS), pages 41–50, 2006.
  * [2] Brendan PW Ames.  Guaranteed recovery of planted cliques and dense subgraphs by convex relaxation.  Journal of Optimization Theory and Applications, 167(2):653–675, 2015.
  * [3] Albert Angel, Nikos Sarkas, Nick Koudas, and Divesh Srivastava.  Dense subgraph maintenance under streaming edge weight updates for real-time story identification.  Proceedings of the VLDB Endowment, 5(6):574–585, 2012.
  * [4] Oana Denisa Balalau, Francesco Bonchi, TH Chan, Francesco Gullo, and Mauro Sozio.  Finding subgraphs with maximum total density and limited overlap.  In Proceedings of the Eighth ACM International Conference on Web Search and Data Mining, pages 379–388. ACM, 2015.
  * [5] Shujun Bi, Xiaolan Liu, and Shaohua Pan.  Exact penalty decomposition method for zero-norm minimization based on mpec formulation.  SIAM Journal on Scientific Computing, 36(4):A1451–A1477, 2014.
  * [6] Jonathan Borwein and Adrian S Lewis.  Convex analysis and nonlinear optimization: theory and examples.  Springer Science & Business Media, 2010.
  * [7] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine learning, 3(1):1–122, 2011.
  * [8] Stephen Boyd and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [9] Tat-Seng Chua, Jinhui Tang, Richang Hong, Haojie Li, Zhiping Luo, and Yan-Tao Zheng.  Nus-wide: A real-world web image database from national university of singapore.  In Proc. of ACM Conf. on Image and Video Retrieval (CIVR’09), Santorini, Greece., July 8-10, 2009.
  * [10] Yunchao Gong, Svetlana Lazebnik, Albert Gordo, and Florent Perronnin.  Iterative quantization: A procrustean approach to learning binary codes for large-scale image retrieval.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(12):2916–2929, 2013.
  * [11] Jie Gui, Tongliang Liu, Zhenan Sun, Dacheng Tao, and Tieniu Tan.  Fast supervised discrete hashing.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(2):490–496, 2018.
  * [12] Lifang He, Chun-Ta Lu, Jiaqi Ma, Jianping Cao, Linlin Shen, and Philip S Yu.  Joint community and structural hole spanner detection via harmonic modularity.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 875–884. ACM, 2016.
  * [13] Cho-Jui Hsieh, Nagarajan Natarajan, and Inderjit S Dhillon.  Pu learning for matrix completion.  In ICML, pages 2445–2453, 2015.
  * [14] Mengqiu Hu, Yang Yang, Fumin Shen, Ning Xie, and Heng Tao Shen.  Hashing with angular reconstructive embeddings.  IEEE Transactions on Image Processing, 27(2):545–555, 2018.
  * [15] David S Johnson and Michael R Garey.  Computers and intractability: A guide to the theory of NP-completeness, volume 1.  WH Freeman San Francisco, 1979.
  * [16] Nikos Komodakis and Georgios Tziritas.  Approximate labeling via graph cuts based on linear programming.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(8):1436–1453, 2007.
  * [17] Guoyin Li and Ting Kei Pong.  Global convergence of splitting methods for nonconvex composite optimization.  SIAM Journal on Optimization, 25(4):2434–2460, 2015.
  * [18] Guosheng Lin, Chunhua Shen, Qinfeng Shi, Anton Van den Hengel, and David Suter.  Fast supervised hashing with decision trees for high-dimensional data.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1963–1970, 2014.
  * [19] Guosheng Lin, Chunhua Shen, David Suter, and Anton Van Den Hengel.  A general two-step approach to learning-based hashing.  In Proceedings of the IEEE international conference on computer vision, pages 2552–2559, 2013.
  * [20] Li Liu, Ling Shao, Fumin Shen, and Mengyang Yu.  Discretely coding semantic rank orders for supervised image hashing.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1425–1434, 2017.
  * [21] Wei Liu, Cun Mu, Sanjiv Kumar, and Shih-Fu Chang.  Discrete graph hashing.  In Advances in Neural Information Processing Systems (NIPS), pages 3419–3427, 2014.
  * [22] Zhaosong Lu and Yong Zhang.  Sparse approximation via penalty decomposition methods.  SIAM Journal on Optimization, 23(4):2448–2478, 2013.
  * [23] Yadan Luo, Yang Yang, Fumin Shen, Zi Huang, Pan Zhou, and Heng Tao Shen.  Robust discrete code modeling for supervised hashing.  Pattern Recognition, 75:128–135, 2018.
  * [24] Sanjay Mehrotra.  On the implementation of a primal-dual interior point method.  SIAM Journal on Optimization, 2(4):575–601, 1992.
  * [25] Walter Murray and Kien-Ming Ng.  An algorithm for nonlinear optimization problems with binary variables.  Computational Optimization and Applications, 47(2):257–288, 2010\.
  * [26] Aude Oliva and Antonio Torralba.  Modeling the shape of the scene: A holistic representation of the spatial envelope.  International Journal of Computer Vision, 42(3):145–175, 2001.
  * [27] Carl Olsson, Anders P Eriksson, and Fredrik Kahl.  Solving large scale binary quadratic problems: Spectral methods vs. semidefinite programming.  In 2007 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8. IEEE, 2007.
  * [28] Ning Qian.  On the momentum term in gradient descent learning algorithms.  Neural networks, 12(1):145–151, 1999.
  * [29] Sekharipuram S Ravi, Daniel J Rosenkrantz, and Giri Kumar Tayi.  Heuristic and special case algorithms for dispersion problems.  Operations Research, 42(2):299–310, 1994.
  * [30] Polina Rozenshtein, Francesco Bonchi, Aristides Gionis, Mauro Sozio, and Nikolaj Tatti.  Finding events in temporal networks: Segmentation meets densest-subgraph discovery.  In 2018 IEEE International Conference on Data Mining (ICDM), pages 397–406. IEEE, 2018.
  * [31] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  arXiv preprint arXiv:1609.04747, 2016.
  * [32] Fumin Shen, Chunhua Shen, Wei Liu, and Heng Tao Shen.  Supervised discrete hashing.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 37–45, 2015.
  * [33] Fumin Shen, Yan Xu, Li Liu, Yang Yang, Zi Huang, and Heng Tao Shen.  Unsupervised deep hashing with similarity-adaptive and discrete optimization.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2018\.
  * [34] Xinchu Shi, Haibin Ling, Junling Xing, and Weiming Hu.  Multi-target tracking by rank-1 tensor approximation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2387–2394, 2013.
  * [35] Jingdong Wang, Ting Zhang, Nicu Sebe, Heng Tao Shen, et al.  A survey on learning to hash.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(4):769–790, 2018.
  * [36] Meihong Wang and Fei Sha.  Information theoretical clustering via semidefinite programming.  In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 761–769, 2011.
  * [37] Peng Wang, Chunhua Shen, Anton van den Hengel, and Philip HS Torr.  Large-scale binary quadratic optimization using semidefinite relaxation and applications.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(3):470–485, 2017.
  * [38] Yu Wang, Wotao Yin, and Jinshan Zeng.  Global convergence of admm in nonconvex nonsmooth optimization.  Journal of Scientific Computing, pages 1–35, 2015.
  * [39] Yair Weiss, Antonio Torralba, and Rob Fergus.  Spectral hashing.  In Advances in Neural Information Processing Systems (NIPS), pages 1753–1760, 2009.
  * [40] Stephen J Wright.  Coordinate descent algorithms.  Mathematical Programming, 151(1):3–34, 2015.
  * [41] Baoyuan Wu and Bernard Ghanem.  lp-box admm: A versatile framework for integer programming.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2018\.
  * [42] Huan Xiong, Mengyang Yu, Li Liu, Fan Zhu, Jie Qin, Fumin Shen, and Ling Shao.  A generalized method for binary optimization: Convergence analysis and applications.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [43] Ganzhao Yuan and Bernard Ghanem.  Binary optimization via mathematical programming with equilibrium constraints.  arXiv preprint arXiv:1608.04425, 2016.
  * [44] Ganzhao Yuan and Bernard Ghanem.  Sparsity constrained minimization via mathematical programming with equilibrium constraints.  arXiv preprint arXiv:1608.04430, 2016.
  * [45] Ganzhao Yuan and Bernard Ghanem.  An exact penalty method for binary optimization based on mpec formulation.  In AAAI, pages 2867–2875, 2017.
  * [46] Xiao-Tong Yuan and Tong Zhang.  Truncated power method for sparse eigenvalue problems.  Journal of Machine Learning Research, 14(Apr):899–925, 2013.
  * [47] Bin Zhang and Steve Horvath.  A general framework for weighted gene co-expression network analysis.  Statistical applications in genetics and molecular biology, 4(1), 2005.
  * [48] Zhongyuan Zhang, Tao Li, Chris Ding, and Xiangsun Zhang.  Binary matrix factorization with applications.  In Data Mining, 2007. ICDM 2007. Seventh IEEE International Conference on, pages 391–400. IEEE, 2007.