# An ODE approach to multiple choice polynomial programming

Sihong Shao111CAPT, LMAM and School of Mathematical Sciences, Peking
University, Beijing 100871, China. Email<EMAIL_ADDRESS>Yishan
Wu222CAPT, LMAM and School of Mathematical Sciences, Peking University,
Beijing 100871, China. Email<EMAIL_ADDRESS>
###### Abstract

We propose an ODE approach to solving multiple choice polynomial programming
(MCPP) after assuming that the optimum point can be approximated by the
expected value of so-called thermal equilibrium as usually did in simulated
annealing. The explicit form of the feasible region and the affine property of
the objective function are both fully exploited in transforming the MCPP
problem into the ODE system. We also show theoretically that a local optimum
of the former can be obtained from an equilibrium point of the latter.
Numerical experiments on two typical combinatorial problems, MAX-$k$-CUT and
the calculation of star discrepancy, demonstrate the validity of our ODE
approach, and the resulting approximate solutions are of comparable quality to
those obtained by the state-of-the-art heuristic algorithms but with much less
cost. This paper also serves as the first attempt to use a continuous
algorithm for approximating the star discrepancy.

AMS subject classifications: 90C59 35Q90, 90C09, 90C23,

Keywords: Pseudo-Boolean optimization; Multiple choice constraint; Continuous
approach; MAX-CUT; Star dicrepancy

##  1 Introduction

We consider the following pseudo-Boolean optimization problem

$\displaystyle\min_{x\in\\{0,1\\}^{n}}$ $\displaystyle f(x),$ (1.1) s.t.
$\displaystyle\sum_{i\in I_{j}}x_{i}=1,\ j=1,2,\dots m,$

where $f$ is a polynomial function, $x$ is an $n$-dimensional Boolean vector,
the indices $[n]:=\\{1,2,\dots,n\\}$ are divided into $m$ disjoint subsets
$I_{1},I_{2},\dots,I_{m}$, and the cardinality of each $I_{j}$, denoted by
$d_{j}:=|I_{j}|$, must be greater than $1$. Then $x$ is accordingly divided
into $m$ vectors $x^{(1)},x^{(2)},\dots,x^{(m)}$ where each $x^{(j)}$ picks
out all the entries in $I_{j}$ of $x$. The constraint means that, for each
$j$, there exists exact one element that equals to $1$ in the subvector
$x^{(j)}$, implying that only one is determined from $d_{j}$ choices and thus
$m$ items chosen from $n$ choices in total. A group of entries of $x$ in a
single $I_{j}$ represents a decision out of finite choices. More precisely, we
use multiple choice polynomial programming (MCPP) to call the problem (1.1)
[8], which is capable of dealing with various problems in diverse disciplines
[1], for example, the MAX-$k$-CUT [16] and star discrepancy [6] problems. As a
typical $0$-$1$ programming problem, it can also be treated with standard
integer programming solvers, such as IBM-CPLEX, SCIP, GLPK, etc. However, in
doing so we have to transform MCPP instances into linear forms by defining new
variable to represent each monomial, thereby greatly increasing the problem
size. This defect becomes even severer in calculating the star discrepancy
since the degree of corresponding objective function is nothing but the
dimension of underlying space and is usually much greater than 1 (See Eq.
(6.15) in Section 6).

Alternatively, developing continuous approaches for discrete problems has
gotten more attention since, e.g., the Hopfield network [10], an early ODE
approach, was proposed for the $0$-$1$ quadratic programming in 1980s. The
Hopfield network has been extended to other combinatorial optimization
problems, see, e.g. [11, 23], where many useful mathematical techniques in
dynamical system have been introduced. In mostly recent work [18], an ODE
approach was proposed through a quartic penalty approximation of the Boolean
polynomial program. In line with this, here we propose to use the solutions of
the following ODE system

$\left\\{\begin{array}[]{ll}\displaystyle\frac{\mathrm{d}{y^{(j)}_{i}}}{\mathrm{d}t}=-y^{(j)}_{i}+\sigma_{i}(-\Phi^{(j)}(y);1/T),\quad
j\in[m],\ i\in[d_{j}],\\\ \displaystyle
y(0)=y_{0}\in[0,1]^{n},\end{array}\right.$ (1.2)

to approximate the solutions of MCPP, where the time-dependent vector
$y(t):[0,+\infty)\rightarrow\mathbb{R}^{n}$ is divided as $x$ in Eq. (1.1)
does, the initial data $y_{0}$ is required to satisfy the “continuous”
multiple choice constraint: $\sum_{i=1}^{d_{j}}(y_{0})_{i}^{(j)}=1$ for all
$j\in[m]$, $\Phi^{(j)}$ denotes the partial derivative of $f$ with respect to
$x^{(j)}$: $\Phi^{(j)}_{i}={\partial f}/{\partial x_{i}^{(j)}}$, $T$ is a
positive parameter called “temperature”, and
$\sigma_{i}(z;\beta):\mathbb{R}^{d}\times\mathbb{R}_{+}\rightarrow(0,1)$ gives
the softmax function defined by

$\sigma_{i}(z;\beta)=\frac{\exp\left(\beta
z_{i}\right)}{\sum_{k=1}^{d}\exp\left(\beta z_{k}\right)},$ (1.3)

with $d$ being the dimension of input argument $z$. Note in passing that
$y^{(j)}_{i}(t)\in(0,1)$ for arbitrary $t>0$, i.e.,
$y(t):[0,+\infty)\rightarrow(0,1)^{n}$. We are able to prove that the
equilibrium points of the ODE system (1.2) represent the local optimum
solutions of MCPP (1.1). Therefore, we just have to numerically integrate
(1.2) until we find an equilibrium (usually fast with fewer iterations) which
can be rounded into a local solution of MCPP. Along the way, various well-
established techniques in numerically solving ODEs can be incorporated to
improve the efficiency without sacrificing the solution quality. Actually, the
ODE system (1.2) can be regarded as a continuous version of simulated
annealing (SA) [13] by starting from a continuous-time Markov chain and
exploiting two intrinsic properties of MCPP (1.1). One is that the explicit
form of the feasible region is so straightforward that the state space of the
Markov chain can be easily determined. The other is that the objective
function of MCPP, $f(x^{(1)},x^{(2)},\dots,x^{(m)})$, is affine with respect
to $x^{(j)}$ for all $j\in[m]$ (otherwise we can modify it equivalently), with
which the definition of transition rate can be significantly simplified. It
should be noted that Eq. (1.2) may recover the ODE system used in [19] for the
unconstrained binary quadratic programming (see Example 2.1).

We apply the proposed ODE approach into two typical NP-hard problems:
MAX-$k$-CUT with $k=2,3,4,5$ and the star discrepancy calculation, and two
state-of-the-art algorithms: the multiple operator heuristic (MOH) [16] and
the improved threshold accepting (TA_improved) [6] methods are employed to as
the reference, respectively. For MAX-$k$-CUT problems, the test bed is
G-Set***available from https://web.stanford.edu/~yyye/yyye/Gset/. More than
half of the ratios between the best cut values achieved by our ODE approach
and those by MOH are above 0.99. For star discrepancy calculation, the test
set is a group of good lattice point (GLP) sets sets adopted in [24]. The
ratios between the best lower bounds achieved by our ODE approach and those by
TA_improved are at least $0.91$. More importantly, the ODE approach requires
much fewer iterations to reach the equilibrium points for each trial while the
cost of each iteration is linear to the problem size. Its iteration steps over
MOH’s for MAX-$k$-CUT are about $1/10^{5}$ to $1/10^{4}$, and about $1/10^{2}$
to $1/10$ over TA_improved’s for calculating the star discrepancy.

The paper is organized as follows. In Section 2 we show how we transform the
MCPP problem (1.1) into the ODE system (1.2) step by step starting from the
well-known assumption (2.2) of simulated annealing. In Section 3 we prove that
a local optimum of MCPP (1.1) can be indeed obtained via an equilibrium point
of the ODE system (1.2). After that, Section 4 details the procedure of
finding the equilibrium points by numerically integrating Eq. (1.2). In
Sections 5 and 6, we apply the proposed ODE approach into the MAX-$k$-CUT
problem and approximating the star discrepancy, respectively. The paper is
concluded in Section 7 with a few remarks.

##  2 The ODE approach: A continuous version of simulated annealing

Under the multiple choice constraint in Eq. (1.1), $x^{(j)}$ must be a
standard unit vector in $\mathbb{R}^{d_{j}}$. Let
$B_{d}=\\{e_{i}^{d}\\}_{i=1}^{d}\subset\mathbb{R}^{d}$ collect the standard
base of $\mathbb{R}^{d}$ where $e^{d}_{i}$ denotes the $d$-dimensional unit
vector the $i$-th entry of which equals $1$ and the rests are $0$s. The
multiple choice constraint is then equivalent to $x^{(j)}\in B_{d_{j}},\
j=1,2,\dots,m$. Therefore, the feasible region for MCPP is

$\mathcal{X}:=B_{d_{1}}\times B_{d_{2}}\times\dots\times B_{d_{m}}.$ (2.1)

We should mention that it is $d_{j}$, the cardinality of each $I_{j}$, that
determines the structure of the feasible region, rather than the particular
devision, and a compact notation, $B_{d}^{m}:=\mathcal{X}$, will be used in
Sections 5 and 6 when $d_{1}=\dots=d_{m}=d$. Then MCPP (1.1) becomes
$\min_{x\in\mathcal{X}}f(x)$ starting from which we explain how regard the ODE
system (1.2) as a continuous version of SA. The well-known SA method [13]
approaches “thermal equilibrium” by generating a Markov chain such that its
stationary distribution is the Boltzmann distribution, namely

$\Pr(X=x)\propto\exp(-E(x)/T),$ (2.2)

where $X\in\mathcal{X}$ is a random variable, and $E(x)$ denotes the energy at
state $x$. We set $E(x):=f(x)$ for MCPP, adopt a similar approach as SA but in
a continuous, ODE-type way, by defining a continuous-time Markov chain, and
then calculate the dynamics of expected value through the forward equation.

To define a random variable of a continuous-time Markov chain, we need to
determine two things: the state space of the variable and the transition rate
between states. For MCPP, the former is nothing but the feasible region (2.1).
Let $X_{t}$ with $t\geqslant 0$ be the random variable. The state change of
$X_{t}$ is allowed to happen only in a single $B_{j}$. That is, $X_{t}$
changes between $x=(x^{(1)},x^{(2)},\dots,x^{(m)})$ and $x^{\prime}$, where
there exists $j\in[m]$ and $i^{\prime}\in[d_{j}]$ such that
$x^{\prime}=(x^{(\bar{j})},e_{i^{\prime}}^{d_{j}})$. Here
$x^{(\bar{j})}:=(x^{(1)},\dots,x^{(j-1)},x^{(j+1)},\dots,x^{(m)})$ denotes the
$(n-d_{j})$-dimensional vector by deleting $x^{(j)}$, and
$(x^{(\bar{j})},\xi):=(x^{(1)},\dots,x^{(j-1)},\xi,x^{(j+1)},\dots,x^{(m)})$
the $n$-dimensional vector by replacing $x^{(j)}$ with
$\xi\in\mathbb{R}^{d_{j}}$. We are ready to determine the transition rate.
Before that, we would like to assume that $f(x^{(1)},x^{(2)},\dots,x^{(m)})$
is affine with respect to $x^{(j)}$ after fixing $x^{(\bar{j})}$ for all
$j\in[m]$:

$f(x)=x^{(j)}\cdot\Phi^{(j)}(x^{(\bar{j})})+f((x^{(\bar{j})},0)),$ (2.3)

where the dot $\cdot$ gives the standard inner product, and the second RHS
term is independent of $x^{(j)}$. Eq. (2.3) can be easily verified if noting:
(1) every $(x^{(j)}_{i})^{k}$ in $f$ can be replaced by $x^{(j)}_{i}$ because
the power of $0$ or $1$ is equal to itself; (2) every monomial with divisor
$x^{(j)}_{i}x^{(j)}_{i^{\prime}}$, $i\neq i^{\prime}$, vanishes; (3) the
partial derivative $\Phi^{(j)}(x)$ is independent of $x^{(j)}$ and thus
$\Phi^{(j)}(x^{(\bar{j})})=\Phi^{(j)}(x)$ is defined unambiguously. For
$x=(x^{(\bar{j})},e_{i}^{d_{j}})$,
$x^{\prime}=(x^{(\bar{j})},e_{i^{\prime}}^{d_{j}})$, $i\neq i^{\prime}$,
according to Eq. (2.3), we have

$f(x)=f((x^{(\bar{j})},0))+\Phi^{(j)}_{i}(x^{(\bar{j})}),\quad
f(x^{\prime})=f((x^{(\bar{j})},0))+\Phi^{(j)}_{i^{\prime}}(x^{(\bar{j})}),$

and

$-f(x^{\prime})+f(x)=-\Phi^{(j)}_{i^{\prime}}(x^{(\bar{j})})+\Phi^{(j)}_{i}(x^{(\bar{j})}).$
(2.4)

Let $q\left(x\rightarrow x^{\prime}\right)$ denote the transition rate from
state $x$ to $x^{\prime}$. It also means $(x,x^{\prime})$-entry of the
$Q$-matrix over $\mathbb{R}^{\mathcal{X}\times\mathcal{X}}$. The detailed
balance condition reads

$\exp(-f(x)/T)q\left(x\rightarrow
x^{\prime}\right)=\exp(-f(x^{\prime})/T)q\left(x^{\prime}\rightarrow
x\right),$

or

$\frac{q\left(x\rightarrow x^{\prime}\right)}{q\left(x^{\prime}\rightarrow
x\right)}=\exp((-f(x^{\prime})+f(x))/T)$ (2.5)

for non-zero entries of $q$ which implies $x\neq x^{\prime}$. To satisfy Eq.
(2.5), combining it with Eq. (2.4), we define

$q\left(x\rightarrow
x^{\prime}\right)=\begin{cases}\sigma_{i^{\prime}}\left(-\Phi^{(j)}(x^{(\bar{j})});1/T\right),&x=(x^{(\bar{j})},e_{i}^{d_{j}}),x^{\prime}=(x^{(\bar{j})},e_{i^{\prime}}^{d_{j}}),i\neq
i^{\prime},\\\ 0,&\text{otherwise.}\end{cases}$ (2.6)

Using $x^{\prime(\bar{j})}=x^{(\bar{j})}$, we have

$q\left(x^{\prime}\rightarrow
x\right)=\sigma_{i}\left(-\Phi^{(j)}(x^{\prime(\bar{j})});1/T\right)=\sigma_{i}\left(-\Phi^{(j)}(x^{(\bar{j})});1/T\right),$

and then Eq. (2.5) can be readily verified.

Next, let’s calculate the dynamics of expected value. We denote the
probability distribution at time $t$ by $p(x;t):=\Pr(X_{t}=x)$. Then it
satisfies the forward equation

$\frac{\mathrm{d}p(x;t)}{\mathrm{d}t}=\sum_{x^{\prime}\in\mathcal{X},x^{\prime}\neq
x}q\left(x^{\prime}\rightarrow
x\right)p\left(x^{\prime};t\right)-\sum_{x^{\prime}\in\mathcal{X},x^{\prime}\neq
x}q\left(x\rightarrow x^{\prime}\right)p\left(x;t\right).$ (2.7)

Let $\mu:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be an arbitrary vector function.
From Eq. (2.7), we have

$\displaystyle\frac{\mathrm{d}\left<\mu(X_{t})\right>}{\mathrm{d}t}$
$\displaystyle=\sum_{x\in\mathcal{X}}\mu(x)\frac{\mathrm{d}p(x;t)}{\mathrm{d}t}$
(2.8) $\displaystyle=\sum_{x^{\prime}\in\mathcal{X}}p(x^{\prime};t)\sum_{x\neq
x^{\prime}}\mu(x)q(x^{\prime}\rightarrow
x)-\sum_{x\in\mathcal{X}}p(x;t)\sum_{x^{\prime}\neq x}\mu(x)q(x\rightarrow
x^{\prime})$ $\displaystyle=\left<\sum_{x\neq
X_{t}}(\mu(x)-\mu(X_{t}))q(X_{t}\rightarrow x)\right>,$

where the term in the summation only counts when both $\mu(x)\neq\mu(X_{t})$
and $x\neq X_{t}$ hold. Let $m_{t}=\left<X_{t}\right>\in[0,1]^{n}$ be the
expected value of $X_{t}$. Then
$(m_{t})^{(j)}_{i}=\sum_{x\in\mathcal{X}}x^{(j)}_{i}p(x;t)$ which assigns
$\mu$ to $(\cdot)^{(j)}_{i}:x\mapsto x^{(j)}_{i}$ where there are only few
$x$s satisfies these two conditions. Recalling the definition of $q$ in Eq.
(2.6), $x\neq x^{\prime}$ iif $x$ and $x^{\prime}$ differs in only one group
$I_{k}$. In addition, we need $x^{(j)}_{i}\neq(x^{\prime})^{(j)}_{i}$ then $k$
must be $j$. Therefore, $x=((X_{t})^{(\bar{j})},e^{d_{j}}_{i^{\prime}})$ for
some $e^{d_{j}}_{i^{\prime}}\neq(X_{t})^{(j)}$. Accordingly, from Eqs. (2.8)
and (2.6), we obtain

$\displaystyle\frac{\mathrm{d}(m_{t})^{(j)}_{i}}{\mathrm{d}t}$
$\displaystyle=\left<\sum_{1\leqslant i^{\prime}\leqslant
d_{j},e^{d_{j}}_{i^{\prime}}\neq(X_{t})^{(j)}}(x^{(j)}_{i}-(X_{t})^{(j)}_{i})q(X_{t}\rightarrow
x)\right>$ $\displaystyle=\left<\sum_{1\leqslant i^{\prime}\leqslant
d_{j},e^{d_{j}}_{i^{\prime}}\neq(X_{t})^{(j)}}(x^{(j)}_{i}-(X_{t})^{(j)}_{i})\sigma_{i^{\prime}}(-\Phi^{(j)}(x);1/T)\right>$
$\displaystyle=\left<\sum_{1\leqslant i^{\prime}\leqslant
d_{j}}(x^{(j)}_{i}-(X_{t})^{(j)}_{i})\sigma_{i^{\prime}}(-\Phi^{(j)}({x^{(\bar{j})}});1/T)\right>$
$\displaystyle=\left<\sum_{1\leqslant i^{\prime}\leqslant
d_{j}}x^{(j)}_{i}\sigma_{i^{\prime}}(-\Phi^{(j)}({(X_{t})^{(\bar{j})}});1/T)\right>-\left<\sum_{1\leqslant
i^{\prime}\leqslant
d_{j}}(X_{t})^{(j)}_{i}\sigma_{i^{\prime}}(-\Phi^{(j)}({x^{(\bar{j})}});1/T)\right>$
$\displaystyle=\left<\sigma_{i}(-\Phi^{(j)}(X_{t});1/T)\right>-\left<(X_{t})^{(j)}_{i}\right>$
(2.9)
$\displaystyle\approx\sigma_{i}(-\Phi^{(j)}(m_{t});1/T)-(m_{t})^{(j)}_{i},$
(2.10)

where we have applied $x^{(\bar{j})}=(X_{t})^{(\bar{j})}$, and
$x^{(j)}_{i}\neq 0$ iif $i^{\prime}=i$ in the first RHS term of Eq. (2.9),
$\sum_{1\leqslant i^{\prime}\leqslant
d_{j}}\sigma_{i^{\prime}}(-\Phi^{(j)}(x);1/T)=1$ in the second RHS term of Eq.
(2.9), as well as a rough approximation
$\left<\sigma_{i}(-\Phi^{(j)}(X_{t});1/T)\right>\approx\sigma_{i}(-\Phi^{(j)}(m_{t});1/T)$
in the first RHS term of Eq. (2.10). The ODE system (1.2) is manifest in Eq.
(2.10) after replacing $m_{t}$ with $y(t)\in(0,1)^{n}$. That is, $y(t)$ is
also an approximation of the dynamics of the expected value of a continuous-
time Markov chain, and by searching for its efficient numerical
approximations, we are able to approximate the Markov chain efficiently.
Moreover, its equilibrium serves as the stable distribution of the Markov
chain, namely, the Boltzmann distribution (2.2). When the temperature goes to
zero, the distribution approaches the uniform distribution on the ground
states, which may give the optimum solution. Therefore, the equilibrium gives
a reasonable approximation of the optimum, the proof of which is left for
Section 3.

Finally, $X_{t}\in\mathcal{X}$ means that the multiple choice constraint
$\sum_{i=1}^{d_{j}}(X_{t})^{(j)}_{i}=1$ for all $j\in[m]$ holds with
probability $1$ and so does $m_{t}$. In fact, we claim that $y(t)$ defined in
Eq. (1.2) also satisfy such constraint when the initial data $y(0)$ does,
which can be readily seen from the sum of Eq. (1.2):

$\frac{\mathrm{d}}{\mathrm{d}t}S_{j}(t)=-S_{j}(t)+1\quad\Leftrightarrow\quad
S_{j}(t)=(S_{j}(0)-1)\mathbbm{e}^{-t}+1,$ (2.11)

where $S_{j}(t)=\sum_{i=1}^{d_{j}}y^{(j)}_{i}(t)$.

######  Example 2.1.

The ODE approach is applicable for the unconstrained situation where the
multiple choice constraint in Eq. (1.1) is absent. By introducing extra $n$
Boolean variables $x_{n+1},\dots,x_{2n}$, the unconstrained problem
$\min_{x\in\\{0,1\\}^{n}}f(x)$ can be reformulated into the following MCPP

$\displaystyle\min_{x\in\\{0,1\\}^{2n}}$ $\displaystyle
f(x_{1},x_{2},\dots,x_{n}),$ (2.12) s.t. $\displaystyle x_{i}+x_{i+n}=1,\
i=1,2,\dots,n,$

where the index set $[2n]$ are divided into $n$ disjoint subsets:
$I_{j}=\\{j,j+n\\}$, $d_{j}=2$, $j=1,2,\dots,n$, and
$\Phi^{(j)}_{1}(x)=\frac{\partial f(x)}{\partial x_{j}}$,
$\Phi^{(j)}_{2}(x)=0$. Then from Eq. (1.2), we have that $y_{j}$ (i.e.,
$y^{(j)}_{1}$) satisfies

$\frac{\mathrm{d}y_{j}}{\mathrm{d}t}=-y_{j}+\frac{1}{2}\left(\tanh\left(-\frac{1}{2T}\frac{\partial
f(y)}{\partial x_{j}}\right)+1\right),\quad j=1,2,\dots,n.$ (2.13)

As $y_{1}^{(j)}(t)+y_{2}^{(j)}(t)=y_{j}(t)+y_{j+n}(t)=1$ for $t\geqslant 0$,
we do not have to care about the dynamics of extra variable $y_{j+n}(t)$.
Implementing a linear transformation $z_{j}=2y_{j}-1$ in Eq. (2.13) yields a
ODE sysytem for $z\in[-1,1]^{n}$:

$\frac{\mathrm{d}z_{j}}{\mathrm{d}t}=-z_{j}+\tanh\left(-\frac{1}{2T}\frac{\partial
f(\frac{1}{2}(z+1))}{\partial x_{j}}\right),\quad j=1,2,\dots,n,$

which was also mentioned in [19] for the unconstrained binary quadratic
programming.

##  3 Numerical analysis

Let $\bar{y}(T)$ be an equilibrium point of the ODE system (1.2) under the
parameter $T$, namely,

$\bar{y}(T)^{(j)}=\sigma\left(-\Phi^{(j)}(\bar{y}(T));1/T\right),\
\forall\,j\in[m].$ (3.1)

We claim that $\bar{y}(T)$ approximates a local solution of Eq. (1.1) when $T$
is sufficiently small. This can be formally explained as follows. As $\beta$
approaches $+\infty$, the limit of softmax function (1.3), denoted by
$\hat{\sigma}(z)=(\hat{\sigma}_{i}(z))$, is a kind of “hard”max:

$\hat{\sigma}_{i}(z)=\lim_{\beta\rightarrow+\infty}\sigma_{i}(z;\beta)=\begin{cases}1/r,&z_{i}=\max\\{z_{1},z_{2},\dots,z_{d}\\},\\\
0,&\text{otherwise},\end{cases}$ (3.2)

where $r=|\\{i|z_{i}=\max\\{z_{1},z_{2},\dots,z_{d}\\}\\}|$ denotes the number
of maximal entries of $z$. It can be easily observed that the range of
$\hat{\sigma}(z)$, denoted by $\bar{B}_{d}$, is a discrete set, i.e.,

$\bar{B}_{d}=\left\\{\frac{1}{|A|}\chi^{A}\in\mathbb{R}^{d}\middle|\chi^{A}_{i}=\left\\{\begin{array}[]{l}1,\,i\in
A\\\ 0,\,i\notin
A,\end{array}\right.A\subseteq[d],\;\;A\neq\varnothing\right\\}.$ (3.3)

After assuming the limit $\hat{y}:=\lim_{T\rightarrow+\infty}\bar{y}(T)$
exists, it can be formally deduced from Eq. (3.1) that

$\hat{y}^{(j)}=\hat{\sigma}\left(-\Phi^{(j)}(\hat{y})\right),\
\forall\,j\in[m].$ (3.4)

That is, $\hat{y}^{(j)}$ belongs to $\bar{B}_{d_{j}}$, and thus

$\hat{y}\in\bar{\mathcal{X}}:=\bar{B}_{d_{1}}\times\bar{B}_{d_{2}}\times\dots\bar{B}_{d_{m}}.$
(3.5)

Combining Eqs. (3.2) and (3.4), we know that the nonzero entries in
$\hat{y}^{(j)}$ correspond to the minimal entries of $\Phi^{(j)}(\hat{y})$,
and thus $\hat{y}^{(j)}$ minimizes the MCPP objective function
$f((\hat{y}^{(\bar{j})},\xi))$ for $\xi\in\mathbb{R}^{d_{j}}$ after fixing
$\hat{y}^{(\bar{j})}$ in view of its affine property given in Eq. (2.3).
Actually, given Definition 3.1, we are able to prove that $\hat{y}$ is indeed
a local optimum of MCPP (1.1) in $\bar{\mathcal{X}}$ (see Proposition 3.2).

######  Definition 3.1 (local optimality).

We call $x\in\bar{\mathcal{X}}$ a local optimum of $f$ in $\bar{\mathcal{X}}$,
if $\forall\,x^{\prime}\in\bar{\mathcal{X}}$ satisfying $\exists\,j_{0}\in[m]$
such that $(x^{\prime})^{(\bar{j_{0}})}=x^{(\bar{j_{0}})}$, we have
$f(x^{\prime})\geqslant f(x)$.

######  Proposition 3.2.

The $n$-dimensional vector $\hat{y}$ defined in Eq. (3.4) is a local optimum
of MCPP (1.1) in $\bar{\mathcal{X}}$.

###### Proof.

Let $y^{\prime}\in\bar{\mathcal{X}}$ satisfy
$(y^{\prime})^{(\bar{j_{0}})}=\hat{y}^{(\bar{j_{0}})}$. From Eqs. (3.2) and
(3.4), we have

$\hat{y}^{(j_{0})}_{i}=\frac{1}{r},\;\;\Phi_{i}^{(j_{0})}(\hat{y}^{\bar{(j_{0})}})=\min_{1\leqslant
k\leqslant d_{j_{0}}}\Phi_{k}^{(j_{0})}(\hat{y}^{(\bar{j_{0}})}),\;\;\forall\
i\ \text{such that }\hat{y}^{(j_{0})}_{i}\neq 0,$ (3.6)

where $r$ equals the number of nonzero entries of $\hat{y}^{(j_{0})}$.
Combining Eqs. (2.3) and (3.6) yields

$\displaystyle f(y^{\prime})-f(\hat{y})$
$\displaystyle=(y^{\prime})^{(j_{0})}\cdot\Phi^{(j_{0})}((y^{\prime})^{(\bar{j_{0}})})-\hat{y}^{(j_{0})}\cdot\Phi^{(j_{0})}(\hat{y}^{(\bar{j}_{0})})$
$\displaystyle=(y^{\prime})^{(j_{0})}\cdot\Phi^{(j_{0})}(\hat{y}^{(\bar{j_{0}})})-\min_{1\leqslant
k\leqslant d_{j_{0}}}\Phi_{k}^{(j_{0})}(\hat{y}^{(\bar{j_{0}})})\geqslant 0,$

where the fact $(y^{\prime})^{({j_{0}})}\in\bar{B}_{d_{j_{0}}}$ is used in the
last inequality. ∎

It should be noted that Eq. (3.4) does not hold in general since
$\sigma(\cdot;\beta)$ does not converge to $\hat{\sigma}$ uniformly. However,
in numerical experiments, the equilibriums $\bar{y}(T)$, which may not be
close to a Boolean vector, are always very close to vectors in
$\bar{\mathcal{X}}$ for sufficiently small $T$. Therefore, we can reasonably
assume that $\bar{y}$ is an approximation of $\hat{y}$ and that
$\sigma(\cdot;1/T)$ is an approximation of $\hat{\sigma}$. Then Eq. (3.1)
approximates Eq. (3.4). In fact, the theoretical results exists and in
Proposition 3.3 we give the analytical condition for the closeness between
$\bar{y}$ and $\bar{\mathcal{X}}$ under which we can claim that Eq. (3.4)
holds. Thus, by Proposition 3.2, $\hat{y}$ is a local optimum. Before that, we
give some notations and definitions.

Let $L$ be the maximal Lipschitz constant of all $\Phi_{i}^{(j)}$ with
$j\in[m]$, $i\in[d_{j}]$. Namely, we have

$\left|\Phi^{(j)}_{i}(x)-\Phi_{i}^{(j)}(x^{\prime})\right|\leqslant
L\|x-x^{\prime}\|_{\infty},\quad\forall\,j\in[m],\forall\,i\in[d_{j}],\forall\,x,x^{\prime}\in[0,1]^{n}.$
(3.7)

Moreover, we define the “minimal gap” of $\Phi$ as

$\mathsf{g}:=\min\left\\{\left|\Phi_{i^{\prime}}^{(j)}(x)-\Phi_{i}^{(j)}(x)\right|\middle|\Phi_{i^{\prime}}^{(j)}(x)\neq\Phi_{i}^{(j)}(x),x\in\bar{\mathcal{X}},j\in[m],i^{\prime},i\in[d_{j}]\right\\}$
(3.8)

to measure how close can two different $\Phi^{(j)}_{i}$s in
$\bar{\mathcal{X}}$ get, and $\hat{d}:=\max\\{d_{1},d_{2},\dots,d_{m}\\}$.

######  Proposition 3.3.

Let $\bar{y}=\bar{y}(T)$ be an equilibrium of the ODE system (1.2), and assume
that there exists $\hat{y}\in\bar{\mathcal{X}}$, whose distance to $\bar{y}$,
denoted by $\varepsilon:=\|\hat{y}-\bar{y}\|_{\infty}>0$, satisfies

$\hat{d}\varepsilon<\frac{1}{2},\;\;\;\frac{\varepsilon}{\ln(1/(\hat{d}\varepsilon)-1)}<\frac{T}{2L},\;\;\;T\ln\frac{1+\hat{d}\varepsilon}{1-\hat{d}\varepsilon}+2L\varepsilon<\mathsf{g}.$
(3.9)

Then we have Eq. (3.4) holds and thus $\hat{y}$ must be a local optimum of $f$
in $\bar{\mathcal{X}}$.

###### Proof.

Let $A_{j}$ collect the nonzero entries of $\hat{y}^{(j)}$ and $r=|A_{j}|$. In
order to obtain Eq. (3.4), according to the definition of softmax function
given in Eq. (1.3), it suffices to show for every $j\in[m]$,

$\Phi_{i}^{(j)}(\hat{y})\begin{cases}=\min\limits_{1\leqslant k\leqslant
d_{j}}\Phi_{k}^{(j)}(\hat{y}),&i\in A_{j},\\\ >\min\limits_{1\leqslant
k\leqslant d_{j}}\Phi_{k}^{(j)}(\hat{y}),&i\notin A_{j}.\end{cases}$ (3.10)

And, the conditions given in Eq. (3.9) imply that $\hat{d}\varepsilon<1$,
$\ln(1/(\hat{d}\varepsilon)-1)>0$, and
$T\ln(1/(\hat{d}\varepsilon)-1)>2L\varepsilon$.

First, we prove that
$\Phi_{i_{1}}^{(j)}(\hat{y})=\Phi_{i_{2}}^{(j)}(\hat{y})$, $\forall\
i_{1},i_{2}\in A_{j}$. Using the fact that $\bar{y}$ is an equilibrium, we
arrive at

$\ln\frac{\bar{y}_{i_{1}}^{(j)}}{\bar{y}_{i_{2}}^{(j)}}=\frac{1}{T}\left(-\Phi_{i_{1}}^{(j)}(\bar{y})+\Phi_{i_{2}}^{(j)}(\bar{y})\right).$
(3.11)

On the other hand, the distance between $\hat{y}$ and $\bar{y}$ provides a
limit for every entry of $\bar{y}^{(j)}$. By
$|\hat{y}^{(j)}_{i}-\bar{y}^{(j)}_{i}|\leqslant\|\hat{y}-\bar{y}\|_{\infty}=\varepsilon$,
we get

$\bar{y}^{(j)}_{i}\in\begin{cases}\left[\frac{1}{r}-\varepsilon,\frac{1}{r}+\varepsilon\right],&i\in
A_{j},\\\ [0,\varepsilon],&i\notin A_{j}.\end{cases}$ (3.12)

Combining Eqs. (3.9), (3.11) and (3.12) together and using $r\leqslant\hat{d}$
lead to an estimate,

$\left|\Phi_{i_{1}}^{(j)}(\bar{y})-\Phi_{i_{2}}^{(j)}(\bar{y})\right|=T\left|\ln\frac{\bar{y}_{i_{1}}^{(j)}}{\bar{y}_{i_{2}}^{(j)}}\right|\leqslant
T\ln\frac{\frac{1}{r}+\varepsilon}{\frac{1}{r}-\varepsilon}\leqslant
T\ln\frac{1+\hat{d}\varepsilon}{1-\hat{d}\varepsilon},\;\;\forall\
i_{1},i_{2}\in A_{j},$

and thus we have

$\displaystyle\left|\Phi_{i_{1}}^{(j)}(\hat{y})-\Phi_{i_{2}}^{(j)}(\hat{y})\right|$
$\displaystyle\leqslant\left|\Phi_{i_{1}}^{(j)}(\hat{y})-\Phi_{i_{1}}^{(j)}(\bar{y})\right|+\left|\Phi_{i_{1}}^{(j)}(\bar{y})-\Phi_{i_{2}}^{(j)}(\bar{y})\right|+\left|\Phi_{i_{2}}^{(j)}(\hat{y})-\Phi_{i_{2}}^{(j)}(\bar{y})\right|$
$\displaystyle\leqslant
L\varepsilon+T\ln\frac{1+\hat{d}\varepsilon}{1-\hat{d}\varepsilon}+L\varepsilon<\mathsf{g},$

which directly implies
$\Phi_{i_{1}}^{(j)}(\hat{y})=\Phi_{i_{2}}^{(j)}(\hat{y})$ for all
$i_{1},i_{2}\in A_{j}$, because $\mathsf{g}$ defined in Eq. (3.8) gives the
minimal distance between two different $\Phi^{(j)}_{i}$s in
$\bar{\mathcal{X}}$.

In order to verify Eq. (3.10), the rest is to prove that
$\Phi_{i^{\prime}}^{(j)}(\hat{y})>\Phi_{i_{1}}^{(j)}(\hat{y})$ for all
$i^{\prime}\not\in A_{j}$. According to Eqs. (3.9), (3.11) and (3.12) and the
fact that $r\leqslant\hat{d}$, we have

$\Phi_{i^{\prime}}^{(j)}(\bar{y})-\Phi_{i_{1}}^{(j)}(\bar{y})=T\ln\frac{\bar{y}_{i_{1}}^{(j)}}{\bar{y}_{i^{\prime}}^{(j)}}\geqslant
T\ln\frac{\frac{1}{r}-\varepsilon}{\varepsilon}\geqslant
T\ln\left(\frac{1}{\hat{d}\varepsilon}-1\right)>2L\varepsilon,$

and then

$\displaystyle\Phi_{i^{\prime}}^{(j)}(\hat{y})-\Phi_{i_{1}}^{(j)}(\hat{y})$
$\displaystyle=\left(\Phi_{i^{\prime}}^{(j)}(\bar{y})-\Phi_{i_{1}}^{(j)}(\bar{y})\right)+\left(\Phi_{i^{\prime}}^{(j)}(\hat{y})-\Phi_{i^{\prime}}^{(j)}(\bar{y})\right)-\left(\Phi_{i_{1}}^{(j)}(\hat{y})-\Phi_{i_{1}}^{(j)}(\bar{y})\right)$
$\displaystyle>2L\varepsilon-\left|\Phi_{i^{\prime}}^{(j)}(\hat{y})-\Phi_{i^{\prime}}^{(j)}(\bar{y})\right|-\left|\Phi_{i_{1}}^{(j)}(\hat{y})-\Phi_{i_{1}}^{(j)}(\bar{y})\right|$
$\displaystyle\geqslant 2L\varepsilon-L\varepsilon-L\varepsilon=0.$

The proof is completed. ∎

######  Remark 3.4.

There exists a more concise condition for $\varepsilon$:

$\varepsilon<\min\left\\{\frac{1}{4\hat{d}},\frac{T}{2L},\frac{\mathsf{g}}{3\hat{d}T+2L}\right\\},$
(3.13)

from which the three conditions given in Eq. (3.9) can be readily derived by
direct relaxations. That is, for a given $T$, a sufficiently small
$\varepsilon$ guarantees an equilibrium point of the ODE system (1.2) gives a
local optimum of MCPP (1.1) in $\bar{\mathcal{X}}$. However, we cannot
theoretically assure that $\varepsilon$ will be small enough as $T$ approaches
$0$ at current stage.

##  4 Implementation details

The biggest benefit of transforming MCPP (1.1) into the ODE system (1.2) is
that various numerical ODE solvers and related techniques come into play for
complex combinatorial problems. From Proposition 3.3, we need to find the
equilibrium of Eq. (1.2) quickly rather than to know how we approach it. That
is, we do not have to know the dynamics evolution quite precisely. To this
end, we adopt variable time step forward Euler scheme for accelerating the
calculation. Besides, also inspired by SA, we use an infinite temperatures
sequence (mostly it is decreasing) $T_{1},T_{2},\dots,$ and numerically solve
Eq. (1.2) with $T:=T_{s},\ s=1,2,\dots,$ until we find an equilibrium
$\bar{y}(T_{s})$. Except for the first try with $T=T_{1}$, the initial value
for the ODE system (1.2) with $T=T_{s}$ can be set to be the previous
equilibrium $\bar{y}(T_{s-1})$. When the distance
$\varepsilon:=\min_{y\in\bar{\mathcal{X}}}\|y-\bar{y}(T)\|_{\infty}$ is less
than a prescribed accuracy tolerance $\varepsilon_{0}$, we stop evolving the
ODE system and turn to rounding procedure.

It should be noted that the initial temperature $T_{1}$ should be chosen
carefully to make sure that the first equilibrium $\bar{y}(T_{1})$ is not a
uninformative equilibrium such that every entry of $\bar{y}(T_{1})^{(j)}$ is
close to $1/d_{j}$, otherwise the diversity of approximate MCPP solutions may
be largely reduced, and sometimes the algorithm even produces the same
solutions for different initial inputs. To see this, let’s consider the
extreme case by fomally letting $T=+\infty$ and replacing $1/T$ by $0$ in Eq.
(3.1). We get $\bar{y}(+\infty)^{(j)}=\bm{1}_{d_{j}}/d_{j}$, where
$\bm{1}_{d_{j}}$ is a $d_{j}$-dimensional vector with every entry identical to
$1$. A similar phenomenon happens when $T_{1}$ is large enough, in which the
numerical experiments show that we can only find one equilibrium,
$\bar{y}(T_{1})$, which is close to $\bar{y}(+\infty)$. For some actual
problems such as the MAX-$k$-CUT problem, $\bar{y}(+\infty)$ is even a trivial
equilibrium for all $T>0$.

###  4.1 Variable time step

The forward Euler (FE) scheme for an autonomous system
$\frac{\mathrm{d}y}{\mathrm{d}t}=F(y)$ reads

$y^{k}=y^{k-1}+h^{k-1}F^{k-1},$ (4.1)

where $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a vector function,
$y^{k}$ denotes the numerical approximation of $y(t)$ at time $t^{k}$,
$h^{k-1}:=t^{k}-t^{k-1}$ is the time step, and $F^{k}:=F(y^{k})$. To estimate
the local truncation error of Eq. (4.1), which is dominated by $(h^{k-1})^{2}$
when $h^{k-1}$ is small enough, we adopt a method similar to the Richardson
extrapolation [20]. Let $(h^{k-1})^{2}c$ denote the error. Assuming that the
vector $c$ keeps almost unchanged for some successive time steps and
$h^{k-1}=h^{k-2}$ for even $k$, we have another approximation to $y(t^{k})$,
$\tilde{y}^{k}:=y^{k-2}+(h^{k-1}+h^{k-2})F^{k-2}=y^{k-2}+2h^{k-1}F^{k-2}$, the
error of which is $4(h^{k-1})^{2}c$, while Eq. (4.1) has the error
$(h^{k-1})^{2}c+(h^{k-2})^{2}c=2(h^{k-1})^{2}c$ staring from the same
$y^{k-2}$. That is, we may use $(\tilde{y}^{k}-y^{k})/2\approx(h^{k-1})^{2}c$
to estimate the truncation error in practice. Specifically, we adopt a
tolerance $\Theta>0$ for $\theta^{k}:=\|\tilde{y}^{k}-y^{k}\|_{2}$ and a step
adjust ratio $\rho>1$. Since the error for the $k$-th step is proportion to
the square of step size $h^{k-1}$ and can be well approximated by
$\theta^{k}/2$, we are able to maintain it within a specific range according
to $\Theta$ and $\rho$ by a two-way adjustment of $h^{k}$ as follows. Reduce
the time step $h^{k+1}=h^{k}=h^{k-1}/\rho$ when $\theta^{k}>\Theta\rho^{2}$,
which may avoid the error increase by decreasing the time step; Increase the
time step $h^{k+1}=h^{k}=\rho h^{k-1}$ when $\theta^{k}<\Theta/\rho^{2}$,
which may exploit the maximum large time step allowed by a small error. If the
errors make a small change in successive two steps, then $\theta^{k}$ will not
be far from $[\Theta/\rho^{2},\rho^{2}\Theta]$. Thus, in doing so, the errors
can be controlled around a desired accuracy of $\Theta/2$, and few of them
will be outside the interval $[\Theta/2\rho^{2},\rho^{2}\Theta/2]$.

###  4.2 Rounding procedure

Proposition 3.3 requires $\varepsilon=\|\hat{y}-\bar{y}(T)\|_{\infty}$ to be
small enough. It is natural to think of letting $\hat{y}$ reach the minimum
among $\bar{\mathcal{X}}$. However, a naive implementation of this procedure
is impractical since each $B_{d_{j}}$ contains $2^{d_{j}}-1$ elements and thus
$|\bar{\mathcal{X}}|=\prod_{j=1}^{m}(2^{d_{j}}-1)$. Instead, we only generate
an appropriate $\hat{y}$ as follows. For each $\bar{y}^{(j)}$, let
$\eta_{j}=\max_{i\in[d_{j}]}\\{\bar{y}_{i}^{(j)}\\}$ and $r_{j}=\left\lfloor
1/\eta_{j}+1/2\right\rfloor$. Then consider the largest $r_{j}$ entries of
$\bar{y}^{(j)}$ and record the indices in set $A_{j}\subseteq[d_{j}]$. This
always can be done after noting $\eta_{j}\geqslant 1/d_{j}$. Let
$\hat{y}^{(j)}=\chi^{A_{j}}/r_{j}$ for each $j\in[m]$ then
$\hat{y}^{(j)}\in\bar{B}_{d_{j}}$ as $|A_{j}|=r_{j}\leqslant d_{j}$, and we
achieve $\hat{y}\in\bar{\mathcal{X}}$. It should be pointed out that if there
exists a $r_{j}^{\prime}$-element set $A_{j}^{\prime}\subseteq[d_{j}]$ such
that $\|\bar{y}^{(j)}-\chi^{A_{j}^{\prime}}/r_{j}^{\prime}\|_{\infty}$ is
small enough, $1/\eta_{j}$ will be close to $r_{j}^{\prime}$, thus we choose
the right $r_{j}=r_{j}^{\prime}$. Moreover, the largest $r_{j}$ entries of
$\bar{y}^{(j)}$ will be close to $1/r_{j}^{\prime}$, which form
$A_{j}^{\prime}$ and are also chosen by $A_{j}$. This explains why
$\hat{y}^{(j)}$ may serve as a good approximation for $\bar{y}^{(j)}$.

At the final step, we have to get a Boolean vector $x$ from $\hat{y}$. First,
we let $x\leftarrow\hat{y}$. Then sequentially change $x^{(j)}$ into a
standard unit vector. This process should be done greedily. More precisely, we
choose and $i\in[d_{j}]$ such that minimize $\Phi_{i}^{(j)}(x^{(\bar{j})})$
and let $x^{(j)}\leftarrow e_{i}^{d_{j}}$. This will not increase $f(x)$. When
$x$ is locally optimal, the procedure ends. Thus $f(x)$ is no larger than
$f(\hat{y})$.

###  4.3 Cost analysis

The main cost of each FE step in Eq. (4.1) is calculating $F(y)$, and
according to Eq. (1.2), the cost of $F(y)$ consists of calculating $\Phi(y)$,
$m$ softmax functions with $\Phi(y)^{(j)}$ as input and $\mathcal{O}(n)$ basic
operations. The softmax functions requires another $\mathcal{O}(n)$
operations. Therefore, the computational complexity of the gradient function
may dominate the cost and it is necessary to analyze it respectively for
different problems in practice.

##  5 Application to MAX-$k$-CUT

As a classical graph optimization problem, the MAX-$k$-CUT problem wants a
$k$-division of the vertex set $V$ for a given edge-weighted graph $G=(V,E,W)$
such that the sum of weights across any two subsets is maximized. Let
$V=\\{v_{i}\\}_{i=1}^{|V|}$, $W=(w_{ij})_{|V|\times|V|}$ and $V_{1}$,
$V_{2},\dots$, $V_{k}$ denote the subsets after division. Precisely,
MAX-$k$-CUT maximizes the following cut function

$\mathrm{cut}(V_{1},V_{2},\dots,V_{k}):=\sum_{v_{i}\in V_{r},v_{j}\in
V_{s},r\neq s,i<j}w_{ij}.$ (5.1)

It reduces to the famous MAX-CUT problem for $k=2$, one of the $21$ Karp’s NP-
complete problems [12]. MAX-$k$-CUT is widely studied and there are a number
of algorithms for finding its approximating solutions, including the heuristic
methods, see e.g. [16] and reference therein, and continuous algorithms [7,
22, 3].

We will use the proposed ODE approach to solve the MAX-$k$-CUT problem and
need to transform it into MCPP at the first step. This can be readily
implemented by setting $n\leftarrow k|V|$, $m\leftarrow|V|$, $d_{j}\leftarrow
k$ for any $j\in[m]$ in Eqs. (1.1) and (2.1), and the feasible region for
MAX-$k$-CUT turns out to be $B_{k}^{|V|}$. Accordingly, the belonging of each
vertex $v_{i}$ is mapped to an element in $B_{k}$ as follows: $v_{i}\in V_{r}$
if and only if $x^{(i)}=e^{k}_{r}$; $v_{i}$ and $v_{j}$ belong to the same
subset if and only if $(x^{(i)})^{\top}x^{(j)}=1$. That is, the corresponding
MCPP objective function becomes

$f(x)=-\mathrm{cut}(V_{1},V_{2},\dots,V_{k})=-\sum_{i<j}w_{ij}\left(1-(x^{(i)})^{\top}x^{(j)}\right).$
(5.2)

Then our ODE approach can be used directly to solve the MAX-$k$-CUT problem.

Tables 1-5 present the numerical results for $k=2,3,4,5$ on G-Set which
includes $71$ randomly generated graph with $800$ to $20000$ vertices, and the
MOH method [16], a state-of-the-art algorithm for MAX-$k$-CUT, is adopted
there to as the reference. We generate 100 random initial values for each
instance and record the maximum cut values in Tables 1-4, see the columns
headed by “ODE”. Other parameters for the proposed ODE approach are $T_{1}=3$,
$T_{s}=\gamma^{s-1}T_{1},s=2,3,\dots$, $\gamma=0.95$, $\varepsilon_{0}=1\times
10^{-3}$, $\Theta=1\times 10^{-5}$ and $\rho=1.1$. For comparison, Tables 1-4
also calculate the ratio between the maximum cut value achieved by our ODE
approach in those $100$ runs and that by MOH for each instance, and find that,
for all four problems, such ratio is larger than $0.99$ for at least $38$
instances. For the remaining graphs of G-Set, the ODE method performs
diversely. For $k=2$, the ratios are at least $0.96$ and in only one instance
is less than $0.97$ (see G64 in Table 1). For $k=3$, the ratios are all above
$0.97$. However, there are $5$ instances for $k=4$ the ratios of which are in
$(0.95,0.97]$ (see G18, G19, G20, G56, G61 in Table 3), and there are $6$
ratios in this interval for $k=5$ (see G19, G20, G21, G40, G56, G61 in Table
4). The worst solution quality happens in solving the MAX-$5$-CUT problem and
the ratio decrease to $0.9472$ for G18 and it is the only one less than $0.95$
in Table 4. For $k=2,4,5$, the ratios reach $0.98$ in at least $56$ instances,
but only $46$ ratios reach $0.98$ for $k=3$. In order to further show the
overall performance in the solution quality achieved by the ODE approach, we
plot the histogram of the ratios for all $4\times 71\times 100$ trials in
Figure 1, and observe there that, in more than $97\%$ trials, the ratios are
greater than $0.95$, thereby implying that the approximate cuts produced by
the ODE approach are of comparable quality to MOH. It should be noted that,
given an initial data, the ODE approach is deterministic, and there is no much
heuristic operations. Hence, we may conclude that, the proposed ODE approach
can produce relatively good results without any ad hoc design for the
MAX-$k$-CUT problem.

Next, we are going to show the efficiency of our ODE approach. Let
$w_{\text{tot}}=\sum_{i<j}w_{ij}$ represent the total weight of all edges and
$P=(p_{ij})\in\\{0,1\\}^{k\times|V|}$ a matrix satisfying
$p_{ij}=x^{(j)}_{i}$. It is easy to check that

$f(x)=-w_{\text{tot}}+\frac{1}{2}\mathrm{Trace}(PWP^{\top})\Rightarrow\frac{\partial
f}{\partial P}=PW,$ (5.3)

where we use the denominator layout, thereby implying that the complexity of
calculating $\nabla f$ is $\mathcal{O}(k|E|)$. To obtain this order we have
used that $W$ is a sparse matrix with $2|E|$ nonzero elements and calculating
each row of $PW$ requires $\mathcal{O}(|E|)$ operations. By the observation in
Section 4.3, we deduce that the time complexity of each FE step is
$\mathcal{O}(k(|E|+|V|))$. We implement the ODE approach with MATLAB R2021a
and run it on AMD Ryzen 5800H. Table 5 reports the approximate total steps and
time needed by running the ODE approach once and shows clearly that the
computing cost is not expensive.

Table 1: Numerical cut values for MAX-CUT. Graph | ODE | MOH | Ratio | Graph | ODE | MOH | Ratio
---|---|---|---|---|---|---|---
G1 |  11623 |  11624 |  0.9999 | G37 |  7629 |  7691 |  0.9919
G2 |  11612 |  11620 |  0.9993 | G38 |  7624 |  7688 |  0.9917
G3 |  11621 |  11622 |  0.9999 | G39 |  2355 |  2408 |  0.9780
G4 |  11646 |  11646 |  1.0000 | G40 |  2342 |  2400 |  0.9758
G5 |  11626 |  11631 |  0.9996 | G41 |  2333 |  2405 |  0.9701
G6 |  2173 |  2178 |  0.9977 | G42 |  2429 |  2481 |  0.9790
G7 |  1999 |  2006 |  0.9965 | G43 |  6645 |  6660 |  0.9977
G8 |  2005 |  2005 |  1.0000 | G44 |  6640 |  6650 |  0.9985
G9 |  2045 |  2054 |  0.9956 | G45 |  6638 |  6654 |  0.9976
G10 |  1993 |  2000 |  0.9965 | G46 |  6640 |  6649 |  0.9986
G11 |  554 |  564 |  0.9823 | G47 |  6644 |  6657 |  0.9980
G12 |  548 |  556 |  0.9856 | G48 |  6000 |  6000 |  1.0000
G13 |  576 |  582 |  0.9897 | G49 |  6000 |  6000 |  1.0000
G14 |  3049 |  3064 |  0.9951 | G50 |  5880 |  5880 |  1.0000
G15 |  3032 |  3050 |  0.9941 | G51 |  3823 |  3848 |  0.9935
G16 |  3027 |  3052 |  0.9918 | G52 |  3825 |  3851 |  0.9932
G17 |  3030 |  3047 |  0.9944 | G53 |  3829 |  3850 |  0.9945
G18 |  979 |  992 |  0.9869 | G54 |  3820 |  3852 |  0.9917
G19 |  891 |  906 |  0.9834 | G55 |  10211 |  10299 |  0.9915
G20 |  931 |  941 |  0.9894 | G56 |  3954 |  4016 |  0.9846
G21 |  914 |  931 |  0.9817 | G57 |  3414 |  3494 |  0.9771
G22 |  13325 |  13359 |  0.9975 | G58 |  19096 |  19288 |  0.9900
G23 |  13337 |  13344 |  0.9995 | G59 |  5911 |  6087 |  0.9711
G24 |  13325 |  13337 |  0.9991 | G60 |  14089 |  14190 |  0.9929
G25 |  13321 |  13340 |  0.9986 | G61 |  5678 |  5798 |  0.9793
G26 |  13308 |  13328 |  0.9985 | G62 |  4766 |  4868 |  0.9790
G27 |  3325 |  3341 |  0.9952 | G63 |  26799 |  27033 |  0.9913
G28 |  3286 |  3298 |  0.9964 | G64 |  8460 |  8747 |  0.9672
G29 |  3379 |  3405 |  0.9924 | G65 |  5444 |  5560 |  0.9791
G30 |  3389 |  3413 |  0.9930 | G66 |  6208 |  6360 |  0.9761
G31 |  3299 |  3310 |  0.9967 | G67 |  6802 |  6942 |  0.9798
G32 |  1384 |  1410 |  0.9816 | G70 |  9482 |  9544 |  0.9935
G33 |  1362 |  1382 |  0.9855 | G72 |  6840 |  6998 |  0.9774
G34 |  1362 |  1384 |  0.9841 | G77 |  9712 |  9928 |  0.9782
G35 |  7628 |  7686 |  0.9925 | G81 |  13724 |  14036 |  0.9778
G36 |  7612 |  7680 |  0.9911 |  |  |  |
Table 2: Numerical cut values for MAX-$3$-CUT. Graph | ODE | MOH | Ratio | Graph | ODE | MOH | Ratio
---|---|---|---|---|---|---|---
G1 |  15158 |  15165 |  0.9995 | G37 |  9968 |  10052 |  0.9916
G2 |  15160 |  15172 |  0.9992 | G38 |  9959 |  10040 |  0.9919
G3 |  15167 |  15173 |  0.9996 | G39 |  2837 |  2903 |  0.9773
G4 |  15176 |  15184 |  0.9995 | G40 |  2808 |  2870 |  0.9784
G5 |  15187 |  15193 |  0.9996 | G41 |  2803 |  2887 |  0.9709
G6 |  2631 |  2632 |  0.9996 | G42 |  2897 |  2980 |  0.9721
G7 |  2401 |  2409 |  0.9967 | G43 |  8571 |  8573 |  0.9998
G8 |  2420 |  2428 |  0.9967 | G44 |  8543 |  8571 |  0.9967
G9 |  2464 |  2478 |  0.9944 | G45 |  8543 |  8566 |  0.9973
G10 |  2401 |  2407 |  0.9975 | G46 |  8545 |  8568 |  0.9973
G11 |  653 |  669 |  0.9761 | G47 |  8548 |  8572 |  0.9972
G12 |  645 |  660 |  0.9773 | G48 |  6000 |  6000 |  1.0000
G13 |  670 |  686 |  0.9767 | G49 |  6000 |  6000 |  1.0000
G14 |  3979 |  4012 |  0.9918 | G50 |  6000 |  6000 |  1.0000
G15 |  3943 |  3984 |  0.9897 | G51 |  4992 |  5037 |  0.9911
G16 |  3951 |  3991 |  0.9900 | G52 |  4996 |  5040 |  0.9913
G17 |  3936 |  3983 |  0.9882 | G53 |  5004 |  5039 |  0.9931
G18 |  1176 |  1207 |  0.9743 | G54 |  4993 |  5036 |  0.9915
G19 |  1050 |  1081 |  0.9713 | G55 |  12329 |  12429 |  0.9920
G20 |  1097 |  1122 |  0.9777 | G56 |  4610 |  4752 |  0.9701
G21 |  1089 |  1109 |  0.9820 | G57 |  3985 |  4083 |  0.9760
G22 |  17080 |  17167 |  0.9949 | G58 |  24977 |  25195 |  0.9913
G23 |  17130 |  17168 |  0.9978 | G59 |  7097 |  7262 |  0.9773
G24 |  17117 |  17162 |  0.9974 | G60 |  16946 |  17076 |  0.9924
G25 |  17125 |  17163 |  0.9978 | G61 |  6668 |  6853 |  0.9730
G26 |  17103 |  17154 |  0.9970 | G62 |  5566 |  5685 |  0.9791
G27 |  3963 |  4020 |  0.9858 | G63 |  35001 |  35322 |  0.9909
G28 |  3921 |  3973 |  0.9869 | G64 |  10217 |  10443 |  0.9784
G29 |  4036 |  4106 |  0.9830 | G65 |  6341 |  6490 |  0.9770
G30 |  4063 |  4119 |  0.9864 | G66 |  7241 |  7416 |  0.9764
G31 |  3963 |  4003 |  0.9900 | G67 |  7904 |  8086 |  0.9775
G32 |  1618 |  1653 |  0.9788 | G70 |  9999 |  9999 |  1.0000
G33 |  1583 |  1625 |  0.9742 | G72 |  8010 |  8192 |  0.9778
G34 |  1573 |  1607 |  0.9788 | G77 |  11332 |  11578 |  0.9788
G35 |  9961 |  10046 |  0.9915 | G81 |  15985 |  16321 |  0.9794
G36 |  9954 |  10039 |  0.9915 |  |  |  |
Table 3: Numerical cut values for MAX-$4$-CUT. Graph | ODE | MOH | Ratio | Graph | ODE | MOH | Ratio
---|---|---|---|---|---|---|---
G1 |  16789 |  16803 |  0.9992 | G37 |  11018 |  11117 |  0.9911
G2 |  16798 |  16809 |  0.9993 | G38 |  11015 |  11108 |  0.9916
G3 |  16799 |  16806 |  0.9996 | G39 |  2935 |  3006 |  0.9764
G4 |  16800 |  16814 |  0.9992 | G40 |  2898 |  2976 |  0.9738
G5 |  16798 |  16816 |  0.9989 | G41 |  2906 |  2983 |  0.9742
G6 |  2739 |  2751 |  0.9956 | G42 |  3014 |  3092 |  0.9748
G7 |  2497 |  2515 |  0.9928 | G43 |  9353 |  9376 |  0.9975
G8 |  2513 |  2525 |  0.9952 | G44 |  9358 |  9379 |  0.9978
G9 |  2574 |  2585 |  0.9957 | G45 |  9360 |  9376 |  0.9983
G10 |  2497 |  2510 |  0.9948 | G46 |  9355 |  9378 |  0.9975
G11 |  661 |  677 |  0.9764 | G47 |  9372 |  9381 |  0.9990
G12 |  652 |  664 |  0.9819 | G48 |  6000 |  6000 |  1.0000
G13 |  680 |  690 |  0.9855 | G49 |  6000 |  6000 |  1.0000
G14 |  4396 |  4440 |  0.9901 | G50 |  6000 |  6000 |  1.0000
G15 |  4354 |  4406 |  0.9882 | G51 |  5518 |  5571 |  0.9905
G16 |  4375 |  4415 |  0.9909 | G52 |  5536 |  5584 |  0.9914
G17 |  4360 |  4411 |  0.9884 | G53 |  5528 |  5574 |  0.9917
G18 |  1207 |  1261 |  0.9572 | G54 |  5527 |  5579 |  0.9907
G19 |  1080 |  1121 |  0.9634 | G55 |  12498 |  12498 |  1.0000
G20 |  1119 |  1168 |  0.9580 | G56 |  4722 |  4931 |  0.9576
G21 |  1133 |  1155 |  0.9810 | G57 |  4036 |  4112 |  0.9815
G22 |  18739 |  18776 |  0.9980 | G58 |  27650 |  27885 |  0.9916
G23 |  18736 |  18777 |  0.9978 | G59 |  7362 |  7539 |  0.9765
G24 |  18740 |  18769 |  0.9985 | G60 |  17148 |  17148 |  1.0000
G25 |  18742 |  18775 |  0.9982 | G61 |  6842 |  7110 |  0.9623
G26 |  18734 |  18767 |  0.9982 | G62 |  5642 |  5743 |  0.9824
G27 |  4159 |  4201 |  0.9900 | G63 |  38758 |  39083 |  0.9917
G28 |  4098 |  4150 |  0.9875 | G64 |  10605 |  10814 |  0.9807
G29 |  4245 |  4293 |  0.9888 | G65 |  6430 |  6534 |  0.9841
G30 |  4261 |  4305 |  0.9898 | G66 |  7352 |  7474 |  0.9837
G31 |  4131 |  4171 |  0.9904 | G67 |  8026 |  8155 |  0.9842
G32 |  1641 |  1669 |  0.9832 | G70 |  9999 |  9999 |  1.0000
G33 |  1602 |  1638 |  0.9780 | G72 |  8123 |  8264 |  0.9829
G34 |  1595 |  1616 |  0.9870 | G77 |  11493 |  11674 |  0.9845
G35 |  11017 |  11111 |  0.9915 | G81 |  16208 |  16470 |  0.9841
G36 |  11023 |  11108 |  0.9923 |  |  |  |
Table 4: Numerical cut values for MAX-$5$-CUT. Graph | ODE | MOH | Ratio | Graph | ODE | MOH | Ratio
---|---|---|---|---|---|---|---
G1 |  17695 |  17703 |  0.9995 | G37 |  11562 |  11603 |  0.9965
G2 |  17693 |  17706 |  0.9993 | G38 |  11552 |  11601 |  0.9958
G3 |  17691 |  17701 |  0.9994 | G39 |  2944 |  3022 |  0.9742
G4 |  17687 |  17709 |  0.9988 | G40 |  2884 |  2986 |  0.9658
G5 |  17697 |  17710 |  0.9993 | G41 |  2925 |  2986 |  0.9796
G6 |  2763 |  2781 |  0.9935 | G42 |  3037 |  3109 |  0.9768
G7 |  2495 |  2533 |  0.9850 | G43 |  9747 |  9770 |  0.9976
G8 |  2516 |  2535 |  0.9925 | G44 |  9742 |  9772 |  0.9969
G9 |  2574 |  2601 |  0.9896 | G45 |  9747 |  9771 |  0.9975
G10 |  2503 |  2526 |  0.9909 | G46 |  9742 |  9774 |  0.9967
G11 |  667 |  677 |  0.9852 | G47 |  9745 |  9775 |  0.9969
G12 |  654 |  662 |  0.9879 | G48 |  6000 |  6000 |  1.0000
G13 |  680 |  689 |  0.9869 | G49 |  6000 |  6000 |  1.0000
G14 |  4611 |  4639 |  0.9940 | G50 |  6000 |  6000 |  1.0000
G15 |  4580 |  4606 |  0.9944 | G51 |  5802 |  5826 |  0.9959
G16 |  4588 |  4613 |  0.9946 | G52 |  5802 |  5837 |  0.9940
G17 |  4575 |  4603 |  0.9939 | G53 |  5807 |  5829 |  0.9962
G18 |  1201 |  1268 |  0.9472 | G54 |  5806 |  5830 |  0.9959
G19 |  1083 |  1132 |  0.9567 | G55 |  12498 |  12498 |  1.0000
G20 |  1132 |  1172 |  0.9659 | G56 |  4733 |  4971 |  0.9521
G21 |  1115 |  1162 |  0.9596 | G57 |  4049 |  4111 |  0.9849
G22 |  19513 |  19553 |  0.9980 | G58 |  29019 |  29105 |  0.9970
G23 |  19513 |  19558 |  0.9977 | G59 |  7443 |  7566 |  0.9837
G24 |  19507 |  19555 |  0.9975 | G60 |  17148 |  17148 |  1.0000
G25 |  19506 |  19554 |  0.9975 | G61 |  6866 |  7188 |  0.9552
G26 |  19504 |  19552 |  0.9975 | G62 |  5650 |  5744 |  0.9836
G27 |  4178 |  4236 |  0.9863 | G63 |  40692 |  40786 |  0.9977
G28 |  4117 |  4182 |  0.9845 | G64 |  10732 |  10896 |  0.9849
G29 |  4255 |  4327 |  0.9834 | G65 |  6428 |  6540 |  0.9829
G30 |  4268 |  4340 |  0.9834 | G66 |  7364 |  7476 |  0.9850
G31 |  4157 |  4211 |  0.9872 | G67 |  8035 |  8165 |  0.9841
G32 |  1644 |  1670 |  0.9844 | G70 |  9999 |  9999 |  1.0000
G33 |  1612 |  1638 |  0.9841 | G72 |  8145 |  8266 |  0.9854
G34 |  1589 |  1615 |  0.9839 | G77 |  11516 |  11687 |  0.9854
G35 |  11547 |  11605 |  0.9950 | G81 |  16227 |  16501 |  0.9834
G36 |  11543 |  11601 |  0.9950 |  |  |  |
$0.9$$0.92$$0.94$$0.96$$0.98$$1$$0$$2{,}000$$4{,}000$$6{,}000$$8{,}000$RatioFrequency Figure 1: Quality check: Histogram of the ratios for all $4\times 71\times 100$ trials of MAX-$k$-CUT. Table 5: Approximate total steps and time in seconds on G-Set used by running the ODE approach once for MAX-$k$-CUT from a given initial data.  $k$ | $|V|$ | $|E|$ |  Number of steps | Time |  Time per step
---|---|---|---|---|---
2 | $8\times 10^{2}$ | $2\times 10^{4}$ | $3.0\times 10^{4}$ | $1.3$ | $4.3\times 10^{-5}$
$2\times 10^{3}$ | $2\times 10^{4}$ | $1.0\times 10^{5}$ | $4.0$ | $4.0\times 10^{-5}$
$1\times 10^{4}$ | $2\times 10^{4}$ | $9.0\times 10^{4}$ | $15$ | $1.7\times 10^{-4}$
$2\times 10^{4}$ | $4\times 10^{4}$ | $1.5\times 10^{5}$ | $30$ | $2.0\times 10^{-4}$
3 | $8\times 10^{2}$ | $2\times 10^{4}$ | $9.1\times 10^{4}$ | $7.8$ | $8.6\times 10^{-5}$
$2\times 10^{3}$ | $2\times 10^{4}$ | $1.7\times 10^{5}$ | $28$ | $1.6\times 10^{-4}$
$1\times 10^{4}$ | $2\times 10^{4}$ | $1.1\times 10^{5}$ | $55$ | $5.0\times 10^{-4}$
$2\times 10^{4}$ | $4\times 10^{4}$ | $1.2\times 10^{5}$ | $84$ | $7.0\times 10^{-4}$
4 | $8\times 10^{2}$ | $2\times 10^{4}$ | $1.2\times 10^{5}$ | $20$ | $1.5\times 10^{-4}$
$2\times 10^{3}$ | $2\times 10^{4}$ | $1.4\times 10^{5}$ | $28$ | $2.0\times 10^{-4}$
$1\times 10^{4}$ | $2\times 10^{4}$ | $1.0\times 10^{6}$ | $60$ | $6.0\times 10^{-4}$
$2\times 10^{4}$ | $4\times 10^{4}$ | $1.2\times 10^{5}$ | $90$ | $7.5\times 10^{-4}$
5 | $8\times 10^{2}$ | $2\times 10^{4}$ | $1.1\times 10^{5}$ | $13$ | $1.2\times 10^{-4}$
$2\times 10^{3}$ | $2\times 10^{4}$ | $1.5\times 10^{5}$ | $30$ | $2.0\times 10^{-4}$
$1\times 10^{4}$ | $2\times 10^{4}$ | $1.2\times 10^{5}$ | $70$ | $5.5\times 10^{-4}$
$2\times 10^{4}$ | $4\times 10^{4}$ | $1.3\times 10^{5}$ | $86$ | $6.6\times 10^{-4}$

##  6 Application to approximating the star discrepancy

For a positive integer $d$, vectors $a,b\in\mathbb{R}^{d}$ and relation
$\triangle\in\\{<,>,\leqslant,\geqslant\\}$, let $a\bigtriangleup b$ represent
that the relation holds for every entry, namely $a_{i}\bigtriangleup b_{i}$
for all $i\in[d]$. We define a closed interval
$[a,b]:=\left\\{u\in\mathbb{R}^{d}\middle|u\geqslant a\;\text{and }u\leqslant
b\right\\}$ for $a<b$. The half-open and open cases can be similarly
expressed. The volume for an interval denotes as
$\operatorname{vol}([a,b]):=\prod(b_{i}-a_{i})$. When $a=0$, we drop off $a$
and simplify it as $\operatorname{vol}(b):=\operatorname{vol}([a,b])$. To
represent a series of vector with indices, we let those indices appear as
superscripts of a single symbol. The subscripts are left for denoting the
entries of these vectors. Given an $N$-point set
$U=\\{u^{1},\dots,u^{N}\\}\subset[0,1)^{d}$,
$u=\left(u_{1},\dots,u_{d}\right)\in[0,1]^{d}$, let

$A(u;U)=\left|[0,u)\cap U\right|,\quad
D(u;U)=\operatorname{vol}(u)-\frac{1}{N}A(u;U).$ (6.1)

The star discrepancy of the $N$-point set $U$ is defined by

$d_{\infty}^{*}(U)=\sup_{u\in[0,1]^{d}}|D(u;U)|,$ (6.2)

which measures the uniformity of $U$, and has been widely used in high
dimensional numerical integration [14, 9] as well as in high dimensional
statistics such as number-theoretical methods [4] and density estimation [15,
21]. However, calculating $d_{\infty}^{*}(U)$ admits NP-hardness [5] and only
a few algorithms were developed. They are the exact algorithm [2], the
threshold-accepting algorithm [24] and its variation, dubbed the TA_improved
algorithm [6]. In this section, we will apply the proposed ODE approach for
MCPP to approximate $d_{\infty}^{*}(U)$. To the best of our knowledge, this is
the first attempt to use a continuous algorithm for approximating the star
discrepancy.

For $j\in[d]$, we define
$\Gamma_{j}(U)=\\{u_{j}^{i}\,\big{|}\,i=1,2,\dots,N\\}$,
$\bar{\Gamma}_{j}(U)=\Gamma_{j}(U)\cup\\{1\\}$,
$\underline{\Gamma}_{j}(U)=\Gamma_{j}(U)\cup\\{0\\}$, and
$\Gamma(U)=\Gamma_{1}(U)\times\cdots\times\Gamma_{d}(U)$ ($\bar{\Gamma}(U)$
and $\underline{\Gamma}(U)$ can be defined in a similar way). It can be
readily observed that, on each open sub-domain of $[0,1]^{d}$ divided by the
grids in $\Gamma(U)$, $A(u;U)$ keeps unchanged and thus $|D(u;U)|$ reaches its
maximum at one of the extreme points of $\operatorname{vol}(u)$, i.e., either
the lower left or upper right corner. Since $\underline{\Gamma}(U)$ and
$\bar{\Gamma}(U)$ respectively collect all the lower left and upper right
corners in $[0,1]^{d}$, an equivalent form for the star discrepancy in Eq.
(6.2) was obtained in [17]:

$d_{\infty}^{*}(U)=\max\left\\{\max_{u\in\bar{\Gamma}(U)}D(u;U),\max_{u\in\underline{\Gamma}(U)}\bar{D}(u;U)\right\\},$
(6.3)

where

$\bar{A}(u;U)=\left|[0,u]\cap
U\right|,\quad\bar{D}(u;U)=\frac{1}{N}\bar{A}(u;U)-\operatorname{vol}(u).$
(6.4)

We should mention that, $\underline{\Gamma}(U)$, the feasible region for
maximizing $\bar{D}(u;U)$, can be slightly narrowed to $\Gamma(U)$ due to the
following reasons. For each $j\in[d]$, if $0\in\Gamma_{j}(U)$, then
$\underline{\Gamma}_{j}(U)$ equals $\Gamma_{j}(U)$. Otherwise,
$\bar{A}(u;U)=0$ if $u_{j}=0$, thus $\bar{D}(u;U)=0<d_{\infty}^{*}(U)$,
indicating that we can drop out that kind of $u$ when maximizing
$\bar{D}(u;U)$. In either case, we can modify the Cartesian product expression
of $\Gamma(U)$ by using $\Gamma_{j}(U)$ instead of
$\underline{\Gamma}_{j}(U)$. Therefore, from Eq. (6.3), we achieve a more
concise form for calculating the star discrepancy,

$d_{\infty}^{*}(U)=\max\left\\{\max_{u\in\bar{\Gamma}(U)}D(u;U),\max_{u\in\Gamma(U)}\bar{D}(u;U)\right\\}.$
(6.5)

The equivalent form (6.5) transforms the star discrepancy into two
optimization problems on discrete sets with which we are able to obtain MCPP
forms.

Let’s start from the first optimization problem
$\max_{u\in\bar{\Gamma}(U)}D(u;U)$ in Eq. (6.5). Inspired by the Cartesian-
product structure of $\bar{\Gamma}(U)$, the feasible region of corresponding
MCPP instance should be $B_{N+1}^{d}$ since the cardinality of
$\bar{\Gamma}_{j}$ is no larger than $N+1$. That is, the key parameters which
shape the MCPP problem in Eqs. (1.1) and (2.1) are: $n\leftarrow(N+1)d$,
$m\leftarrow d$, $d_{j}\leftarrow N+1$ for any $j\in[m]$. However, mapping
$\bar{\Gamma}(U)$ into $B_{N+1}^{d}$ and defining an objective function over
$B_{N+1}^{d}$ that represents $D(u;U)$ need meticulous designs we are about to
state below.

For $j\in[d]$, sort $\bar{\Gamma}_{j}(U)$ into
$\bar{u}_{1j}\leqslant\bar{u}_{2j}\leqslant\cdots\leqslant\bar{u}_{Nj}<\bar{u}_{(N+1)j}=1$,
which also sort $\Gamma_{j}(U)$ in the same order for
$\Gamma_{j}(U)=\\{\bar{u}_{1j},\dots,\bar{u}_{Nj}\\}$. For $i\in[N+1]$, let
$\bar{u}_{\sigma_{ij}j}:=u_{j}^{i}$ record the order of $u_{j}^{i}$ and
$\\{\sigma_{ij}\,\big{|}\,i\in[N+1]\\}$ be the corresponding permutation of
$[N+1]$. When some elements of $\Gamma_{j}(U)$ are identical, such permutation
may not be unique, and it is viable to choose one of them arbitrarily. Let
$x\in B_{N+1}^{d}$. For each $j\in[d]$, there is a unique entry of $x^{(j)}$,
denoted by $x^{(j)}_{s_{j}}$ with $s_{j}\in[N+1]$, that equals $1$. Then we
have a surjection from $x\in B_{N+1}^{d}$ to $u\in\bar{\Gamma}(U)$

$x^{(j)}\rightarrow s_{j}\rightarrow
u_{j}=\bar{u}_{s_{j}j},\quad\forall\,j\in[d].$ (6.6)

Let

$\nu(x)=\prod_{j=1}^{d}\sum_{i=1}^{N+1}x^{(j)}_{i}\bar{u}_{ij},\;\;\;\alpha(x)=\sum_{i=1}^{N}\prod_{j=1}^{d}\sum_{k=\sigma_{ij}+1}^{N+1}x^{(j)}_{k},\;\;\;\delta(x)=\nu(x)-\frac{1}{N}\alpha(x),$
(6.7)

all of which are affine with respect to each $x^{(j)}$. In fact, $\nu(x)$,
$\alpha(x)$ and $\delta(x)$ act as $\operatorname{vol}(u)$, $A(u;U)$ and
$D(u;U)$, respectively.

######  Proposition 6.1.

Let $u\in\bar{\Gamma}(U)$ be the image of $x\in B_{N+1}^{d}$ under the
surjection (6.6). Then we have

$\displaystyle\operatorname{vol}(u)$ $\displaystyle=\nu(x),$ (6.8)
$\displaystyle A(u;U)$ $\displaystyle\leqslant\alpha(x),$ (6.9) $\displaystyle
D(u;U)$ $\displaystyle\geqslant\delta(x).$ (6.10)

Moreover, $D(u;U)$ and $\delta(x)$ admit a deeper connection,

$\max_{u\in\bar{\Gamma}(U)}D(u;U)=\max_{x\in B_{N+1}^{d}}\delta(x).$ (6.11)

###### Proof.

It is easy to verify Eq. (6.8) from the fact
$\sum_{i=1}^{N+1}x^{(j)}_{i}\bar{u}_{ij}=u_{j}$. Next let’s express
$|[0,u)\cap\\{u^{i}\\}|$ by $x$. For each $i\in[N]$, we have

$\displaystyle u^{i}\in[0,u)\Leftrightarrow\ $ $\displaystyle
u_{j}^{i}=\bar{u}_{\sigma_{ij}j}<u_{j}=\bar{u}_{s_{j}j},\;\;\;\forall\,j\in[d],$
$\displaystyle\Rightarrow\ $
$\displaystyle\sigma_{ij}<s_{j},\;\;\;\forall\,j\in[d],$ (6.12)
$\displaystyle\Leftrightarrow\ $
$\displaystyle\sum_{k=\sigma_{ij}+1}^{N+1}x^{(j)}_{k}=1,\;\;\;\forall\,j\in[d].$

Thus
$|[0,u)\cap\\{u^{i}\\}|\leqslant\prod_{j=1}^{d}\sum_{k=\sigma_{ij}+1}^{N+1}x^{(j)}_{k}$
and summing up for all $i\in[N]$ leads to Eq. (6.9). Combining Eqs. (6.1),
(6.7), (6.8) and (6.9) together, we immediately arrive at Eq. (6.10), and thus
$\max_{u\in\bar{\Gamma}(U)}D(u;U)\geqslant\max_{x\in B_{N+1}^{d}}\delta(x)$
for the mapping (6.6) from $B_{N+1}^{d}$ to $\bar{\Gamma}(U)$ is surjective.

To prove Eq. (6.11), we only have to show
$\max_{u\in\bar{\Gamma}(U)}D(u;U)\leqslant\max_{x\in B_{N+1}^{d}}\delta(x)$,
which can be obtained via the fact that, for each $u\in\bar{\Gamma}(U)$, there
exists a preimage of $u$, $\hat{x}=\hat{x}(u)\in B_{N+1}^{d}$, such that
$A(u;U)=\alpha(\hat{x}(u))$. The $\hat{x}$ is defined by letting $s_{j}$ equal
the minimal index that $u_{j}=\bar{u}_{s_{j}j}$ holds and
$\hat{x}_{s_{j}}^{(j)}=1,\ j\in[d]$. Thus
$\bar{u}_{(s_{j}-1)j}<\bar{u}_{s_{j}j}$ whenever $s_{j}>1$, otherwise
contradicting to the minimality of $s_{j}$. When $\sigma_{ij}<s_{j}$, we then
have

$\sigma_{ij}\leqslant
s_{j}-1\Rightarrow\bar{u}_{\sigma_{ij}j}\leqslant\bar{u}_{(s_{j}-1)j}<\bar{u}_{s_{j}j}.$
(6.13)

Therefore, the “$\Rightarrow$” in Eq. (6.12) becomes “$\Leftrightarrow$”,
which results in the equivalence between $|[0,u)\cap\\{u^{i}\\}|$ and
$\prod_{j=1}^{d}\sum_{k=\sigma_{ij}+1}^{N+1}\hat{x}^{(j)}_{k}$, rather than
inequality. Hence, we get $A(u;U)=\alpha(\hat{x}(u))$. Let $u^{*}$ reaches the
maximum of $D(u;U)$. Then

$\max_{u\in\bar{\Gamma}(U)}D(u;U)=D(u^{*};U)=\delta(\hat{x}(u^{*}))\leqslant\max_{x\in
B_{N+1}^{d}}\delta(x).$

∎

Eq. (6.11) gives an MCPP form for the first optimization problem in Eq. (6.5)
and the same treatment can be also applied to the second one. The key
parameters which shape the MCPP problem in Eqs. (1.1) and (2.1) become:
$n\leftarrow Nd$, $m\leftarrow d$, $d_{j}\leftarrow N$ for any $j\in[m]$, and
we need the following functions defined on $B_{N}^{d}$:

$\bar{\nu}(x^{\prime})=\prod_{j=1}^{d}\sum_{i=1}^{N}(x^{\prime})^{(j)}_{i}\bar{u}_{ij},\;\;\;\bar{\alpha}(x^{\prime})=\sum_{i=1}^{N}\prod_{j=1}^{d}\sum_{k=\sigma_{ij}}^{N}(x^{\prime})^{(j)}_{k},\;\;\;\bar{\delta}(x^{\prime})=\frac{1}{N}\bar{\alpha}(x^{\prime})-\bar{\nu}(x^{\prime}),$

to present the MCPP form

$\max_{u\in\Gamma(U)}\bar{D}(u;U)=\max_{x^{\prime}\in
B_{N}^{d}}\bar{\delta}(x^{\prime}),$ (6.14)

the proof of which is very similar to that of Eq. (6.11) and thus skipped
here. In a word, according to Eqs. (6.11) and (6.14), the star discrepancy is
now determined by two MCPP problems

$d_{\infty}^{*}(U)=\max\left\\{\max_{x\in
B_{N+1}^{d}}\delta(x),\max_{x^{\prime}\in
B_{N}^{d}}\bar{\delta}(x^{\prime})\right\\},$ (6.15)

which can be solved approximately by the proposed ODE approach in a
straightforward manner.

Compared to the quadratic objective function of MAX-$k$-CUT problem in Eq.
(5.2), both $\delta(x)$ and $\bar{\delta}(x^{\prime})$ in Eq. (6.15) are of
degree of $d$ and it may be more complicated to store and calculate these
function values as well as their gradients when $d>2$. For example,
calculating the gradients of $\delta$ and $\bar{\delta}$ may involve high
computational cost if handling it inappropriately. It forms the main cost in
using FE scheme (4.1) to solve Eq. (6.15), thereby requiring careful
optimization. More precisely, it is easy to check that for
$i\in[N+1],j\in[d]$,

$\frac{\partial\nu(x)}{\partial
x^{(j)}_{i}}=\bar{u}_{ij}\prod_{j^{\prime}\in[d]\backslash\\{j\\}}\sum_{i^{\prime}=1}^{N+1}x^{(j^{\prime})}_{i^{\prime}}\bar{u}_{i^{\prime}j^{\prime}},\;\;\;\frac{\partial\alpha(x)}{\partial
x^{(j)}_{i}}=\sum_{\begin{subarray}{c}i^{\prime}\in[N]\\\
\sigma_{i^{\prime}j}<i\end{subarray}}\prod_{j^{\prime}\in[d]\backslash\\{j\\}}\sum_{k=\sigma_{i^{\prime}j^{\prime}}+1}^{N+1}x^{(j^{\prime})}_{k}.$
(6.16)

If we calculate the sums and products in $\frac{\partial\alpha(x)}{\partial
x^{(j)}_{i}}$ straightforwardly, there will be at least $\mathcal{O}(N^{2}d)$
operations. Considering that there are $(N+1)d$ entries, the complexity of
$\nabla\delta(x)$ will be at least $\mathcal{O}(N^{3}d^{2})$, which is barely
acceptable. However, we would like to point out that the complexity can be
limited to $\mathcal{O}(Nd)$ according to the following procedure. First, it
should be observed that given $U$, Eq. (6.16) can be simplified as

$\frac{\partial\nu(x)}{\partial
x^{(j)}_{i}}=\bar{u}_{ij}\prod_{j^{\prime}\in[d]\backslash\\{j\\}}b_{j^{\prime}}(x),\;\;\;\frac{\partial\alpha(x)}{\partial
x^{(j)}_{i}}=\sum_{\begin{subarray}{c}i^{\prime}\in[N]\\\
\sigma_{i^{\prime}j}<i\end{subarray}}\prod_{j^{\prime}\in[d]\backslash\\{j\\}}z_{i^{\prime}j^{\prime}}(x),$
(6.17)

where the vector function $b:[0,1]^{n}\rightarrow\mathbb{R}^{d}$ and the
matrix function $z:[0,1]^{n}\rightarrow\mathbb{R}^{N\times d}$ are independent
with $i$ or $j$. Thus $b$ and $z$ can be calculated in advance. To calculate
$z$ efficiently, noting that
$z_{i^{\prime}j^{\prime}}(x)=\sum_{k>\sigma_{i^{\prime}j^{\prime}}}x^{(j^{\prime})}_{k}$
is a partial sum of entries of $x^{(j^{\prime})}$ in reverse order, we can
compute each column of $z$ in a group within a cost of only $\mathcal{O}(N)$,
thereby a cost of $\mathcal{O}(Nd)$ for computing $z$. To get
$g_{i^{\prime}j}(x):=\prod_{j^{\prime}\in[d]\backslash\\{j\\}}z_{i^{\prime}j^{\prime}}(x)$
fast, we replace the multiplication by
$\prod_{j^{\prime}\in[d]}z_{i^{\prime}j^{\prime}}(x)/z_{i^{\prime}j}(x)$,
whose numerator can be reused for fixed $i^{\prime}$. Therefore, the
complexity of computing $g$ from $z$ is also $\mathcal{O}(Nd)$. To get
$\nabla\alpha$ from $g$, the partial sum trick can also be applied in the
summation $\sum_{i^{\prime}:\sigma_{i^{\prime}j}<i}g_{i^{\prime}j}(x)$, albeit
$\sigma_{\cdot j}$ determining the order of entries in $g_{\cdot j}$. Hence,
the total time cost is $\mathcal{O}(Nd)$. In the same spirit, we are able to
get $\nabla\nu$ with cost of $\mathcal{O}(Nd)$. By these means, the complexity
of $\nabla\delta$ is limited to $\mathcal{O}(Nd)$. This is also true for
$\nabla\bar{\delta}$. Since $n=\mathcal{O}(Nd)$ for these two MCPP instances,
the cost of a FE step is also $\mathcal{O}(Nd)$ by the observation in Section
4.3.

The GLP sets used in [24] are adopted for test in this work, and include $30$
small sets with $N$ ranging from $28$ to $487$ and $d$ from $4$ to $6$ as well
as $6$ large sets with $N$ from around $2000$ to $5000$ and $d$ from $6$ to
$11$. The parameters for our ODE approach are $T_{1}=1\times 10^{-4}$,
$T_{s}=\gamma^{s-1}T_{1},s=2,3,\dots$, $\gamma=0.95$, $\varepsilon_{0}=1\times
10^{-3}$, $\Theta=1\times 10^{-6}Nd$ and $\rho=1.1$. Numerical values of the
star discrepancy obtained by re-running the ODE approach 100 times are listed
in Table 6, where those obtained by TA_improved are directly copied from [6].
We also list the ratios of numerical values achieved by two methods there. It
can be seen that except for $4$ out of $30$ small instances, the ODE method
gives the same value as TA_improved did, see
$(N,d)=(312,4),(376,4),(487,4),(73,6)$ in Table 6, and for these four
instances, the ratios are at least $0.98$. For the last $6$ large instances,
the ratios are no smaller than $0.91$. The histogram of the ratios for
$36\times 100$ trials is plotted in Figure 2. In more than $85\%$ trials, the
ratios are greater than $0.9$. Compared to TA_improved’s 100000 iterations per
trial used in [6] and the same $\mathcal{O}(Nd)$ complexity per iteration, our
method requires much less computational resource overall while producing
solutions with similar quality. We also record the approximate typical running
time and steps for different sizes of sets in Table 7. It can be seen there
that the time consumed per iteration is roughly proportional to $Nd$.

Table 6: Numerical values of the star discrepancy on GLP sets. Values in the columns headed by “TA” are directly copied from [6] as reference and obtained there by the TA_improved algorithm. $N$ | $d$ | ODE | TA | Ratio | $N$ | $d$ | ODE | TA | Ratio
---|---|---|---|---|---|---|---|---|---
$145$ | $4$ |  $0.0731$ |  $0.0731$ |  $1.0000$ | $28$ | $6$ |  $0.5360$ |  $0.5360$ |  $1.0000$
$255$ | $4$ |  $0.1093$ |  $0.1093$ |  $1.0000$ | $29$ | $6$ |  $0.2532$ |  $0.2532$ |  $1.0000$
$312$ | $4$ |  $0.0617$ |  $0.0618$ |  $0.9974$ | $35$ | $6$ |  $0.3431$ |  $0.3431$ |  $1.0000$
$376$ | $4$ |  $0.0752$ |  $0.0753$ |  $0.9979$ | $50$ | $6$ |  $0.3148$ |  $0.3148$ |  $1.0000$
$388$ | $4$ |  $0.1297$ |  $0.1297$ |  $1.0000$ | $61$ | $6$ |  $0.1937$ |  $0.1937$ |  $1.0000$
$442$ | $4$ |  $0.0620$ |  $0.0620$ |  $1.0000$ | $73$ | $6$ |  $0.1467$ |  $0.1485$ |  $0.9876$
$448$ | $4$ |  $0.0548$ |  $0.0548$ |  $1.0000$ | $81$ | $6$ |  $0.2500$ |  $0.2500$ |  $1.0000$
$451$ | $4$ |  $0.0271$ |  $0.0271$ |  $1.0000$ | $88$ | $6$ |  $0.2658$ |  $0.2658$ |  $1.0000$
$471$ | $4$ |  $0.0286$ |  $0.0286$ |  $1.0000$ | $90$ | $6$ |  $0.1992$ |  $0.1992$ |  $1.0000$
$487$ | $4$ |  $0.0413$ |  $0.0413$ |  $0.9995$ | $92$ | $6$ |  $0.1635$ |  $0.1635$ |  $1.0000$
$102$ | $5$ |  $0.1216$ |  $0.1216$ |  $1.0000$ | $2129$ | $6$ |  $0.0246$ |  $0.0254$ |  $0.9685$
$122$ | $5$ |  $0.0860$ |  $0.0860$ |  $1.0000$ | $3997$ | $7$ |  $0.0251$ |  $0.0254$ |  $0.9882$
$147$ | $5$ |  $0.1456$ |  $0.1456$ |  $1.0000$ | $3997$ | $8$ |  $0.0242$ |  $0.0254$ |  $0.9528$
$153$ | $5$ |  $0.1075$ |  $0.1075$ |  $1.0000$ | $3997$ | $9$ |  $0.0387$ |  $0.0387$ |  $1.0000$
$169$ | $5$ |  $0.0755$ |  $0.0755$ |  $1.0000$ | $4661$ | $10$ |  $0.0256$ |  $0.0272$ |  $0.9412$
$170$ | $5$ |  $0.0860$ |  $0.0860$ |  $1.0000$ | $4661$ | $11$ |  $0.0259$ |  $0.0283$ |  $0.9152$
$195$ | $5$ |  $0.1574$ |  $0.1574$ |  $1.0000$ |  |  |  |  |
$203$ | $5$ |  $0.1675$ |  $0.1675$ |  $1.0000$ |  |  |  |  |
$235$ | $5$ |  $0.0786$ |  $0.0786$ |  $1.0000$ |  |  |  |  |
$236$ | $5$ |  $0.0582$ |  $0.0582$ |  $1.0000$ |  |  |  |  |
$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$500$$1{,}000$$1{,}500$RatioFrequency Figure 2: Quality check: Histogram of the ratios for all $36\times 100$ trials of approximating discrepancy. Table 7: Approximate total steps and time in milliseconds used by the ODE approach for approximating the star discrepancy in Eq. (6.15). $N$ | $d$ | Time |  Number of steps |  Time per step
---|---|---|---|---
$1.0\times 10^{2}$ | $4$ | $4.0\times 10^{1}$ | $8.0\times 10^{2}$ | $5.0\times 10^{-2}$
$2.5\times 10^{2}$ | $4$ | $6.0\times 10^{1}$ | $8.0\times 10^{2}$ | $6.0\times 10^{-2}$
$5.0\times 10^{2}$ | $4$ | $1.0\times 10^{2}$ | $1.5\times 10^{3}$ | $7.0\times 10^{-2}$
$5.0\times 10^{2}$ | $6$ | $6.0\times 10^{2}$ | $5.0\times 10^{3}$ | $1.2\times 10^{-1}$
$2.0\times 10^{3}$ | $6$ | $3.2\times 10^{3}$ | $8.0\times 10^{3}$ | $4.0\times 10^{-1}$
$4.0\times 10^{3}$ | $10$ | $1.3\times 10^{4}$ | $1.3\times 10^{4}$ | $1.0\times 10^{0}$

##  7 Conclusion and discussion

We proposed an ODE approach for multiple choice polynomial programming and
demonstrated its validity via both theoretical analysis and numerical
experiments. We are conducting analysis on the existence of equilibrium points
and the realizability of conditions that ensures the local optimality, and
trying more advanced numerical techniques for evolving an ODE to its
equilibrium points. We are going to extend the proposed ODE approach to some
kinds of mixed integer programming problems with more constraints rather than
multiple choice.

## Data Availability

The datasets generated during and/or analyzed during the current study are
available from the corresponding author on reasonable request.

## Acknowledgements

This research was supported by the National Key R&D Program of China (Nos.
2020AAA0105200, 2022YFA1005102) and the National Natural Science Foundation of
China (Nos. 12288101, 11822102). SS is partially supported by Beijing Academy
of Artificial Intelligence (BAAI).

## References

  * [1] E. Boros and P. L. Hammer, Pseudo-Boolean optimization, Discrete Appl. Math., 123 (2002), pp. 155–225.
  * [2] P. Bundschuh and Y. Zhu, A method for exact calculation of the discrepancy of low-dimensional finite point sets I, Abh. Math. Semin. Univ. Hambg., 63 (1993), pp. 115–133.
  * [3] S. Burer, R. D. C. Monteiro, and Y. Zhang, Rank-two relaxation heuristics for MAX-CUT and other binary quadratic programs, SIAM J. Optim., 12 (2002), pp. 503–521.
  * [4] K.-T. Fang and Y. Wang, Number-theoretic methods in statistics, Chapman and Hall, London, 1993.
  * [5] M. Gnewuch, A. Srivastav, and C. Doerr, Finding optimal volume subintervals with k points and calculating the star discrepancy are NP-hard problems, J. Complex., 25 (2009), pp. 115–127.
  * [6] M. Gnewuch, M. Wahlström, and C. Winzen, A new randomized algorithm to approximate the star discrepancy based on threshold accepting, SIAM J. Numer. Anal., 50 (2011), pp. 781–807.
  * [7] M. X. Goemans and D. P. Williamson, Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming, J. ACM, 42 (1995), pp. 1115–1145.
  * [8] W. C. Healy, Jr, Multiple choice programming (a procedure for linear programming with zero-one variables), Oper. Res., 12 (1964), pp. 122–138.
  * [9] E. Hlawka, Funktionen von beschränkter variatiou in der theorie der gleichverteilung, Ann. di Mat. Pura ed Appl., 54 (1961), pp. 325–333.
  * [10] J. J. Hopfield, Neurons with graded response have collective computational properties like those of two-state neurons., Proc. Natl. Acad. Sci. USA, 81 (1984), pp. 3088–3092.
  * [11] A. Jagota, Approximating maximum clique with a Hopfield network, IEEE Trans. Neural Netw., 6 (1995), pp. 724–735.
  * [12] R. M. Karp, Reducibility among combinatorial problems, in Complexity of Computer Computations, R. E. Miller, J. W. Thatcher, and J. D. Bohlinger, eds., The IBM Research Symposia Series, Springer US, Boston, MA, 1972, pp. 85–103.
  * [13] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization by simulated annealing, Science, 220 (1983), pp. 671–680.
  * [14] J. F. Koksma, A general theorem from the theory of uniform distribution modulo 1, Mathematica B (Zutphen), 11 (1942), pp. 7–11.
  * [15] D. Li, K. Yang, and W. H. Wong, Density estimation via discrepancy based adaptive sequential partition, in Advances in Neural Information Processing Systems, vol. 29, Curran Associates, Inc., 2016.
  * [16] F. Ma and J.-K. Hao, A multiple search operator heuristic for the max-k-cut problem, Ann. Oper. Res., 248 (2017), pp. 365–403.
  * [17] H. Niederreiter, Discrepancy and convex programming, Ann. di Mat. Pura ed Appl., 93 (1972), pp. 89–97.
  * [18] Y.-S. Niu and R. Glowinski, Discrete dynamical system approaches for Boolean polynomial optimization, J. Sci. Comput., 92 (2022), p. 46.
  * [19] C. Peterson and J. R. Anderson, Neural networks and NP-complete optimization problems; a performance study on the graph bisection problem, Complex Systems, 2 (1988), pp. 58–59.
  * [20] L. F. Richardson and R. T. Glazebrook, The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam, Philos. Trans. Royal Soc. A, 210 (1911), pp. 307–357.
  * [21] S. Shao and Y. Xiong, SPADE: sequential-clustering particle annihilation via discrepancy estimation, May 2020.  arXiv:2005.05129.
  * [22] S. Shao, D. Zhang, and W. Zhang, A simple iterative algorithm for maxcut, Aug. 2019.  arXiv:1803.06496.
  * [23] J. Wang, An improved discrete Hopfield neural network for Max-Cut problems, Neurocomputing, 69 (2006), pp. 1665–1669.
  * [24] P. Winker and K.-T. Fang, Application of threshold-accepting to the evaluation of the discrepancy of a set of points, SIAM J. Numer. Anal., 34 (1997), pp. 2028–2042.