PET short=PET, long=Positron Emission Tomography SE short=SE, long=squeeze and
excitation

11institutetext: Biomedical Engineering, Yale University, New Haven, CT, USA
22institutetext: Radiology and Biomedical Imaging, Yale University, New Haven,
CT, USA

# Simultaneous Denoising and Motion Estimation for Low-dose Gated PET using a
Siamese Adversarial Network with Gate-to-Gate Consistency Learning

Bo Zhou 11 Yu-Jung Tsai 22 Chi Liu 1122

###### Abstract

Gating is commonly used in PET imaging to reduce respiratory motion blurring
and facilitate more sophisticated motion correction methods. In the
applications of low dose PET, however, reducing injection dose causes
increased noise and reduces signal-to-noise ratio (SNR), subsequently
corrupting the motion estimation/correction steps, causing inferior image
quality. To tackle these issues, we first propose a Siamese adversarial
network (SAN) that can efficiently recover high dose gated image volume from
low dose gated image volume. To ensure the appearance consistency between the
recovered gated volumes, we then utilize a pre-trained motion estimation
network incorporated into SAN that enables the constraint of gate-to-gate
(G2G) consistency. With high-quality recovered gated volumes, gate-to-gate
motion vectors can be simultaneously outputted from the motion estimation
network. Comprehensive evaluations on a low dose gated PET dataset of 29
subjects demonstrate that our method can effectively recover the low dose
gated PET volumes, with an average PSNR of 37.16 and SSIM of 0.97, and
simultaneously generate robust motion estimation that could benefit subsequent
motion corrections.

###### Keywords:

Low-dose Gated PET, Denoising, Motion Estimation, Motion Correction

##  1 Introduction

PET is a commonly used functional imaging modality. To obtain high quality
image, a small amount of radioactive tracer is administered to patient,
introducing radiation exposure to both patients and healthcare providers [1].
PET data acquisition typically takes several minutes. During this period,
patient’s breathing inevitably introduces blurring in the lung and abdominal
regions. Respiratory gating facilitated by external motion monitoring devices
are typically used to reduce respiratory-induced motion blurring. However,
each gated image is generated by only a fraction of detected events, leading
to high image noise in each gate. To address the noise issue, previous works
proposed motion correction approaches involving non-rigid image registration
among gated images, and utilize the motion vectors to correct motion by using
all detected events to reduce image noise [2]. In the applications of
radiation dose reduction, reduction of injection dose is the first choice but
will increase the image noise and result in low signal-to-noise ratio (SNR).
In the cases where respiratory gating is performed using low-dose data, the
image noise is further increased, potentially causing errors in motion vector
estimation, which subsequently affects the final motion correction results, as
illustrated in Figure 1. To address this challenge, we aim to simultaneously
tackle the image denoising and motion estimation problems in low-dose gated
PET data.

Figure 1: Illustration of phase gated PET acquisition with 6 gates for both
100% full count and 1.5% count levels. End-expiration gate with the least
intra-gate motion (G4) is used as reference gate. Each low dose gated volume
needs to be denoised and registered to the reference gated volume.

Previous works on denoising low-dose PET can be summarized into two
categories: conventional post-processing [3, 4, 5] and deep learning based
post-processing [6, 7, 8, 9]. Conventional post-processing techniques, such as
Gaussian filtering, is the standard technique to reduce PET image noise, but
has challenge to preserve local structure. More recently, non-local mean
filter [3] and block-matching 4D filter [4] have been proposed to denoise low-
dose PET while better preserving the structural information. Deep learning
based methods, such as deep auto-context CNN [6], 3D cGAN [7], UNet [8], and
GAN [9], were developed for recovering standard-dose PET from low-dose PET.
Compared to conventional methods, these deep learning based methods achieved
promising denoising performance on static low-dose PET. However, none of these
previous studies addressed denoising and motion estimation in low-dose
respiratory gated PET in a unified fashion.

In this work, we proposed a Siamese adversarial network (SAN) with gate-to-
gate consistency learning (G2G) to simultaneously denoise low dose gated
volumes and estimate the motion among the gates. We evaluated our method on a
challenging low dose gated PET dataset with only $1.5\%$ count level. Our
experimental results demonstrated that our proposed method can effectively
reduce the noise while preserving the structural information and improve the
accuracy of motion estimation.

##  2 Problem Formulation

Assuming a phase gated PET exam generates 6 gates with gate 4 as the reference
gate, we denote high-dose PET (HDPET) and low-dose PET (LDPET) gated volumes
as $H_{n},L_{n}\in\mathbb{R}^{h\times w\times d}$ with gate index
$n\in\\{1,2,3,4,5,6\\}$ and volume size $h\times w\times d$. The
transformation predicted between $\\{L_{4},L_{n}\\}$ is expected to be
different from the transformation predicted between $\\{H_{4},H_{n}\\}$ due to
the high noise level of LDPET. Given that the distribution of HDPET is
unknown, our goal is to recover $H_{n}$ from the degraded $L_{n}$. Previous
methods have been trying to solve the inverse problem by finding the
generative model $\mathcal{P}_{D}$ parameterized by $\theta_{D}$ such that
$\mathcal{P}_{D}(\sum_{n=1}^{6}L_{n};\theta_{D})=\bar{H}_{nmc}\approx
H_{nmc}$, where $\bar{H}_{nmc}$ is the non-gated denoised volume with no
motion correction (nmc). Since no motion estimation and corresponding motion
correction component are considered, degradation in the final image can be
expected. Therefore, we aim to tackle these issues by recovering the HDPET
from LDPET for each gate and simultaneously estimate the motion field between
gates. Specifically, we want to find single gate generative model
$\mathcal{P}_{D}$ such that
$\mathcal{P}_{D}(L_{n};\theta_{D})=\bar{H}_{n}\approx H_{n}$ where
$\bar{H}_{n}$ is the recovered HDPET for gate $n$. Then, the motion
transformation between the reference gate (assume to be gate 4 here) and gate
$n$ is estimated by
$\bar{T}_{n}=\mathcal{P}_{R}(\mathcal{P}_{D}(L_{4};\theta_{D}),\mathcal{P}_{D}(L_{n};\theta_{D});\theta_{R})\approx
T_{n}$, where $\bar{T}_{n}$ is the predicted transformation from our motion
estimator $\mathcal{P}_{R}$. In this work, our goal is to obtain the optimal
$\mathcal{P}_{D}$ and $\mathcal{P}_{R}$ for simultaneous denoising and motion
estimation.

##  3 Methods

Figure 2: Our two-stage training procedure consists of: the pre-training of
our motion estimator (R), and Siamese adversarial training of our generator.
Two shared weights generators $G$ learn mapping from LDPET to HDPET, which are
supervised by a structure recovery loss
($\mathcal{L}_{SR}=\mathcal{L}_{1}+\mathcal{L}_{SSIM}+\mathcal{L}_{adv}$), and
a transform consistency loss
($\mathcal{L}_{G2G}=\mathcal{L}_{1}+\mathcal{L}_{KL}$), respectively. Motion
estimator $R$ is pre-trained with the ground truth HDPET, and concatenated to
the generator for end-to-end optimization. Network architecture details are
listed in the supplementary.

The overall pipeline of our method is illustrated in Figure 2. It consists of
three major parts: 1) Siamese generative networks with supervision from our
structure recovery loss; 2) unsupervised motion estimation network; and 3)
gate-to-gate consistency training. The Siamese generator $G$ maps the target
gate LDPET ($L_{tgt}$) and the reference gate LDPET ($L_{ref}$) to the HDPET
space simultaneously, thus generating denoised HDPET gated volumes. The
generator $G$ is first optimized based on the structure recovery loss that
measures the dissimilarity between prediction and ground truth, yielding the
high quality denoised HDPET volumes. In the meantime, the motion estimation
network $R$ is pre-trained using the ground truth HDPET gated volumes $H$, and
concatenated to the Siamese generative networks. By replacing the input for
$R$ with the synthetic HDPET volumes $\hat{H}$ generated by $G$, the joint
network enforces gate-to-gate consistency in the transformed synthetic HDPET
for each target gate, providing additional supervision for training $G$. The
details are as follows.

Siamese Generative Network is illustrated in Figure 2. The Siamese generative
network $G$ with encoding and decoding architecture is firstly supervised by a
$\mathcal{L}_{1}$ loss, a structural similarity (SSIM) loss, and an
adversarial loss to ensure the noise reduction and structure recovery.
Specifically, we use a $\mathcal{L}_{1}$ loss to ensure the general appearance
recovery and a $\mathcal{L}_{SSIM}$ loss to ensure the fine-detailed structure
recovery. $\mathcal{L}_{1}$ loss allows noise suppression and SNR improvement,
at the expense of reduced image sharpness. On the other hand,
$\mathcal{L}_{SSIM}$ loss encourages image to have high contrast, sharpness
and resolution. Given $L_{tgt}$ and $L_{ref}$ the target and reference LDPET
gated volumes respectively, $G$ takes a pair of $[L_{tgt},L_{ref}]$ and
channel-wise concatenates each volumes with anatomical prior CT ($\rho$) to
predict $\bar{H}_{tgt}=G(L_{tgt},\rho;\theta_{G})$ and
$\bar{H}_{ref}=G(L_{ref},\rho;\theta_{G})$ simultaneously. The
$\mathcal{L}_{1}$ loss and the $\mathcal{L}_{SSIM}$ loss can be written as:

$\footnotesize\mathcal{L}_{1}=\sum_{i}||H_{i}-\bar{H}_{i}||,\quad
i\in\\{tgt,ref\\}$ (1)
$\footnotesize\mathcal{L}_{SSIM}=\sum_{i}[1-SSIM(H_{i},\bar{H}_{i})],\quad
i\in\\{tgt,ref\\}$ (2a) $\footnotesize
SSIM(x,y)=\frac{2m_{x}m_{y}+C_{1}}{m_{x}^{2}+m_{y}^{2}+C_{1}}\cdot\frac{2\sigma_{xy}+C_{2}}{\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2}}$
(2b)

where $[m_{x},m_{y}]$ and $[\sigma_{x},\sigma_{y}]$ denote mean and standard
deviation of an image pair $[x,y]$. The cross-covariance of $[x,y]$ is denoted
as $\sigma_{xy}$. $C_{1}$ and $C_{2}$ are constant parameters. The adversarial
loss from the discriminator $D$ provides an indication of discrepancy between
prediction and ground truth as both $G$ and $D$ progressively optimized. Thus,
the adversarial loss is also added to minimize the perceptual difference
between prediction and ground truth from a CNN perspective. We utilize the
adversarial loss in Wasseerstein GAN with gradient penalty (WGAN-GP) to
achieve stable adversarial training [10], which is formulated as:

$\displaystyle\footnotesize\begin{split}\mathcal{L}_{adv}=\sum_{i}\mathbb{E}[D(\bar{H}_{i})]-\mathbb{E}[D(H_{i})]+&\lambda_{gp}\mathbb{E}[(||\nabla_{\ddot{H}_{i}}D(\ddot{H}_{i})||_{2}-1)^{2}],i\in\\{tgt,ref\\}\end{split}$
(3)

where $\ddot{H}$ represents a linear combination of $\bar{H}$ and $H$ with a
weight $t$ uniformly sampled between $0$ and $1$. Thereby, $\lambda_{gp}$
controls the gradient penalty level and is set to $3$ here. The combination of
these three loss functions formulates our Structure Recovery (SR) loss as:

$\footnotesize\mathcal{L}_{SR}=\beta_{1}\mathcal{L}_{1}+\beta_{2}\mathcal{L}_{SSIM}+\beta_{3}\mathcal{L}_{adv}$
(4)

where $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ are loss weights. In our
experiments, we empirically set $\beta_{1}=1$, $\beta_{2}=1$, and
$\beta_{3}=0.2$ for balance training.

Motion estimation network $R$ aims to predict the transformation between
target and reference gated volumes. Here, we use a probabilistic generative
model [11] to predict the transformation, as illustrated in Figure 2’s left
section. Assuming $H_{ref}$ and $H_{tgt}$ are volumes that need to be
registered and the transformation between them is parameterized by a sampled
velocity field $V$, $R$ aims to find the most likely registration field by
optimizing the posterior probability $p(V|H_{ref},H_{tgt})$. Thus, the loss
function for network $R$ can be derived and written as:

$\displaystyle\footnotesize\begin{split}\mathcal{L}_{R}(H_{ref},H_{tgt})=&\frac{1}{K}\sum_{k}||H_{ref}-T\circ
H_{tgt}||+\text{KL}[q_{\theta_{R}}(V|H_{ref},H_{tgt})||p(V)]\end{split}$ (5)

where $K$ is the number of samples in each training batch, $T$ is the
transformation function parameterized by $V\sim
q_{\theta_{R}}(V|H_{ref},H_{tgt})$. The first term minimizes the L1 distance
between reference volume $H_{ref}$ and warped target volume $H_{tgt}$. The
second term ensures the distribution similarity between posterior and prior of
$V$. $\mathcal{L}_{R}$ is the transform consistency loss in Figure 2. During
the inference stage, the predicted $V$ is fed into the scaling and squaring
layer [12] to integrate $V$ over $[0,1]$, and produce the final transformation
$T$. Then, $T$ and the target volume $H_{tgt}$ are inputted into the spatial
transform layer to generate the warped target volume $T\circ H_{tgt}$.
Detailed derivation is in our supplementary.

Gate-to-Gate Consistency Learning The Siamese generator in the first part maps
$L$ to $H$ with SR loss $\mathcal{L}_{SR}$ for individual gates. However, the
appearance consistency constraint between gates is not utilized. A gate-to-
gate consistency should sustain when the gated volumes are registered. The
gate-to-gate consistency learning is achieved by feeding the synthetic pair of
HDPET volumes, $[\bar{H}_{tgt},\bar{H}_{ref}]$ generated using the Siamese
generative network $G$, into the pre-trained motion estimation network $R$
after concatenating these two networks. Therefore, the transformation
prediction process of the joint network can be described as:

$\footnotesize\bar{T}=R(\bar{H}_{ref},\bar{H}_{tgt};\theta_{R})=R(G(L_{ref},\rho;\theta_{G}),G(L_{tgt},\rho;\theta_{G});\theta_{R})$
(6)

Given the transformation $\bar{T}$, we warp the synthetic $\bar{H}_{tgt}$ and
obtain $\bar{T}\circ\bar{H}_{tgt}$. We aim to minimize the distance between
$\bar{T}\circ\bar{H}_{tgt}$ and ground truth $H_{ref}$, such that the
transformed target gated volume and reference gated volume are consistent.
Thus, the gate-to-gate transform consistency loss can be formulated as:

$\displaystyle\footnotesize\begin{split}\mathcal{L}_{G2G}=\frac{1}{K}\sum_{k}||H_{ref}-\bar{T}\circ\bar{H}_{tgt}||+\text{KL}[q_{\theta_{R}}(V|\bar{H}_{ref},\bar{H}_{tgt})||p(V)]\end{split}$
(7)

The first term encourages the gate-to-gate appearance consistency using a
$\mathcal{L}_{1}$ norm and the second term ensures the distribution similarity
between posterior and prior of $V$. $\mathcal{L}_{G2G}$ provides additional
supervision for optimizing $G$ by utilizing the inter-gate relationship. It is
the key in our Siamese network design that enables us to randomly sample pairs
of gated volume, which augments the number of available training data for each
subject to $A_{6}^{2}=30$. Therefore, the denoising and structural recovery
from LDPET to HDPET will be more reliable.

Finally, our full loss function for optimizing $G$ is
$\mathcal{L}_{tot}=\mathcal{L}_{SR}+\mathcal{L}_{G2G}$, which is trained in an
adversarial manner. $G$ and $R$ try to minimize this loss collaboratively,
while $D$ tries to maximize it. To optimize the overall network, we update
$G$, $R$, and $D$ alternatively by: optimizing $D$ with $G$ and $R$ fixed,
then optimizing $G$ with $D$ and $R$ fixed.

###  3.1 Evaluation with Human Data

We collected 29 pancreas 18F-FPDTBZ [13] PET/CT studies with respiration
gating facilitated by the Anzai system. The total acquisition time was $120$
mins for each study. We used phase gating to generate 6 gates for each study.
To eliminate the mismatch between attenuation correction (AC) map and gated
PET, instead of using CT as AC-map, we utilized the maximum likelihood
estimation of activity and attenuation (MLAA) [14] to generated AC-map for
each gated volume to ensure phase-matched attenuation correction, where CT was
used as initial estimation for MLAA iterations. The HDPET volumes were
reconstructed with 100% of the listmode data mimicking high radiation dose
data. The LDPET volumes were reconstructed with 1.5% of the listmode data with
random sampling. Each data was reconstructed into a $400\times 400\times 109$
volume with voxel size of $2.032\times 2.032\times 2.027$ $mm^{3}$. The
central $200\times 200\times 109$ voxels were kept to remove most voxels
outside the human body contour and resized to $128\times 128\times 128$. The
end expiration gate (typically Gate 4) was used as the reference gate since it
shows minimum intra-gate motion.

The dataset were split into training set of 22 studies and test set of 7
studies. The evaluation was performed on the 7 test studies with 6 gated
volumes in each study. For quantitative evaluation, the denoising results were
evaluated by comparing the synthetic HDPET volumes to the ground truth HDPET
volumes using the Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity
Index (SSIM). The motion estimation results were evaluated using the Mean
Vector Euclidean Distance (MVED) that measures the 3D Euclidean distance
between the predicted vector field and the ground truth vector field, which
was defined as the motion field predicted by $R$ using the ground truth HDPET.
For comparative study, we compared our results against the following
algorithms: Gaussian filtering (GAU), Non-local mean filtering (NLM) [3],
Block-matching 4D filtering (BM4D) [4], UNet [8, 15], and cGAN [7].

Figure 3: Sample HD and 1.5% LD PET slices with enlarged subregions using
various denoising methods for Two sample subjects. The corresponding PSNR and
SSIM are indicated at the bottom of the images. Comparison of intensity
profile is also shown on the right. -/+G2G denotes without/with gate-to-gate
consistency learning.

##  4 Results

The qualitative comparison of various denoising methods is shown in Figure 3.
As we can observe on the figure, conventional post-processing methods, such as
NLM [3] and BM4D [4], have difficulties in structural recovery when only 1.5%
of the total counts was considered. The high noise level also introduced
additional artifacts, resulting in inferior performance compared to the
standard Gaussian filtering. In contrast, deep learning based methods achieved
better performance in noise reduction and structural recovery.

Table 1 outlines the quantitative comparison of different methods on PET image
denoising. Both PSNR and SSIM were evaluated for each gated volumes
($\mathcal{G}_{n}$), along with averaged value computed on the last column.
Among them, our SAN without G2G outperforms the previous deep learning based
methods, and the addition of G2G learning that utilizes the information over
gates further improved the performance. In parallel, Figure 4 illustrates a
qualitative comparison of motion estimation based on the discussed denoising
methods. As we can see, our proposed SAN+G2G yields the most consistent motion
vectors between the estimated and ground truth motion vectors. The
quantitative comparison of motion estimation among different denoising methods
is given in Table 2. As shown in the table, our SAN+G2G was able to improve
the motion estimation accuracy by $20\%$ in average, achieving the lowest
$0.264$ in averaged MVED, compared to other studied methods. Using our
proposed method, denoised gated LDPET volumes can be generated with
corresponding motion vectors to the reference gate. We then registered all
gated volumes of LDPET, HDPET, and LDPET with SAN+G2G to the reference gate.
As an example shown in Figure 5, the proposed network is able to generate
gated PET volumes with reduced noise level and a final motion corrected image
that averaged all registered image volumes with reduced motion blurring using
low dose gated data.

Table 1: Quantitative comparison of denoising results using PSNR (dB) and SSIM ($\times 10^{2}$). Among conventional post-processing methods and deep learning based methods, the optimal results are marked in red. PSNR/SSIM | $\mathcal{G}1$ | $\mathcal{G}2$ | $\mathcal{G}3$ | $\mathcal{G}4$ | $\mathcal{G}5$ | $\mathcal{G}6$ |  Average
---|---|---|---|---|---|---|---
LDPET | 28.48/86.3 | 28.47/86.5 | 28.27/86.6 | 28.47/86.9 | 27.93/86.1 | 28.11/86.2 |  28.29/86.4
GAU | 34.92/94.6 | 34.25/94.5 | 34.18/94.5 | 34.32/94.6 | 34.36/94.4 | 34.66/94.6 |  34.45/94.6
NLM [3] | 31.71/94.5 | 31.80/94.5 | 31.39/94.4 | 31.60/94.6 | 31.03/94.3 | 31.23/94.5 |  31.46/94.5
BM4D [4] | 31.32/93.5 | 31.32/93.5 | 31.03/93.5 | 31.21/93.8 | 30.68/93.4 | 30.83/93.4 |  31.07/93.5
UNet [8] | 37.12/95.9 | 36.01/95.8 | 36.24/95.7 | 36.32/96.1 | 36.34/95.9 | 36.85/96.1 |  36.48/95.9
cGAN [7] | 37.38/96.3 | 36.21/96.2 | 36.41/96.1 | 36.43/96.1 | 36.58/96.2 | 37.02/96.2 |  36.67/96.2
SAN-G2G | 37.55/96.8 | 36.48/96.6 | 36.59/96.6 | 36.67/96.7 | 36.96/96.6 | 37.25/96.8 |  36.92/96.7
SAN+G2G | 37.81/97.1 | 36.74/96.9 | 36.77/96.9 | 36.87/96.9 | 37.35/97.0 | 37.43/97.1 |  37.16/97.0

Figure 4: Qualitative motion estimation results from different denoising
methods. Ground truth (green arrows) and predicted (magenta arrows) motion
estimation vectors are overlaid on denoised images.
Figure 5: Illustration of motion blurred images (left) and averaged image of
all gates average registered to the reference frame (middle). The green arrows
indicate where significant motion reduction is observed after applying the
proposed SAN+G2G.

Table 2: Quantitative comparison of motion estimation results evaluated in terms of MVED. G4 is the reference gate. Optimal results are marked in red. MVED | $\mathcal{G}1$ | $\mathcal{G}2$ | $\mathcal{G}3$ |  $\mathcal{G}4$ (Ref) | $\mathcal{G}5$ | $\mathcal{G}6$ |  Average
---|---|---|---|---|---|---|---
LDPET | 0.392 | 0.342 | 0.289 | - | 0.282 | 0.350 |  0.331
GAU | 0.368 | 0.372 | 0.328 | - | 0.328 | 0.352 |  0.349
NLM [3] | 0.360 | 0.363 | 0.312 | - | 0.309 | 0.342 |  0.337
BM4D [4] | 0.351 | 0.341 | 0.2997 | - | 0.287 | 0.331 |  0.322
SAN-G2G | 0.309 | 0.309 | 0.262 | - | 0.263 | 0.297 |  0.288
SAN+G2G | 0.289 | 0.285 | 0.236 | - | 0.237 | 0.274 |  0.264

##  5 Discussion and Conclusion

In this work, we propose a Siamese adversarial network with gate-to-gate
consistency learning, a novel framework for low dose gated PET denoising and
motion estimation, simultaneously. We first pre-train our motion estimation
network on the ground truth HDPET, and concatenate it to our Siamese
adversarial network that enables the gate-to-gate consistency learning for
improving the denoising performance. The denoised low-dose gated volumes are
then fed into the motion estimation network for robust motion estimation. In
our framework, the Siamese input design allows us to efficiently augment the
training data from each patient, thus can better train generalizable denoising
and motion estimation models. We demonstrated the feasibility of our method on
the tasks of PET image denoising and motion estimation with promising
performance.

The potential clinical feasibility of our work is two-fold. Firstly, as high-
noise level and motion are inevitable in the chest and abdominal low-dose PET
acquisitions, it will affect the visualization of small pathological findings,
such as lung/liver lesions. Our work is potentially useful for recovering
these small objects from noise and correcting motions to improve the
delineation of distorted objects. Secondly, the estimated motion can be
incorporated into the motion compensated PET reconstruction frameworks toward
motion-free low-dose PET reconstructions, which will also improve the
reconstruction quality by reducing the motion artifacts. We will explore these
directions in our future works.

## References

  * [1] Strauss, K.J., Kaste, S.C.:  The alara (as low as reasonably achievable) concept in pediatric interventional and fluoroscopic imaging: striving to keep radiation doses as low as possible during fluoroscopy of pediatric patients—a white paper executive summary.  Radiology 240(3) (2006) 621–622
  * [2] Catana, C.:  Motion correction options in pet/mri.  In: Seminars in nuclear medicine. Volume 45., Elsevier (2015) 212–223
  * [3] Dutta, J., Leahy, R.M., Li, Q.:  Non-local means denoising of dynamic pet images.  PloS one 8(12) (2013) e81390
  * [4] Maggioni, M., Katkovnik, V., Egiazarian, K., Foi, A.:  Nonlocal transform-domain filter for volumetric data denoising and reconstruction.  IEEE transactions on image processing 22(1) (2012) 119–133
  * [5] Mejia, J., Mederos, B., Mollineda, R.A., Maynez, L.O.:  Noise reduction in small animal pet images using a variational non-convex functional.  IEEE Transactions on Nuclear Science 63(5) (2016) 2577–2585
  * [6] Xiang, L., Qiao, Y., Nie, D., An, L., Lin, W., Wang, Q., Shen, D.:  Deep auto-context convolutional neural networks for standard-dose pet image estimation from low-dose pet/mri.  Neurocomputing 267 (2017) 406–416
  * [7] Wang, Y., Yu, B., Wang, L., Zu, C., Lalush, D.S., Lin, W., Wu, X., Zhou, J., Shen, D., Zhou, L.:  3d conditional generative adversarial networks for high-quality pet image estimation at low dose.  Neuroimage 174 (2018) 550–562
  * [8] Lu, W., Onofrey, J.A., Lu, Y., Shi, L., Ma, T., Liu, Y., Liu, C.:  An investigation of quantitative accuracy for deep learning based denoising in oncological pet.  Physics in Medicine & Biology 64(16) (2019) 165019
  * [9] Kaplan, S., Zhu, Y.M.:  Full-dose pet image estimation from low-dose pet image using deep learning: a pilot study.  Journal of digital imaging 32(5) (2019) 773–778
  * [10] Arjovsky, M., Chintala, S., Bottou, L.:  Wasserstein gan.  arXiv preprint arXiv:1701.07875 (2017)
  * [11] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.:  Unsupervised learning for fast probabilistic diffeomorphic registration.  In: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer (2018) 729–738
  * [12] Arsigny, V., Commowick, O., Pennec, X., Ayache, N.:  A log-euclidean framework for statistics on diffeomorphisms.  In: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer (2006) 924–931
  * [13] Normandin, M.D., Petersen, K.F., Ding, Y.S., Lin, S.F., Naik, S., Fowles, K., Skovronsky, D.M., Herold, K.C., McCarthy, T.J., Calle, R.A., et al.:  In vivo imaging of endogenous pancreatic $\beta$-cell mass in healthy and type 1 diabetic subjects using 18f-fluoropropyl-dihydrotetrabenazine and pet.  Journal of Nuclear Medicine 53(6) (2012) 908–916
  * [14] Rezaei, A., Michel, C., Casey, M.E., Nuyts, J.:  Simultaneous reconstruction of the activity image and registration of the ct image in tof-pet.  Physics in Medicine & Biology 61(4) (2016) 1852
  * [15] Ronneberger, O., Fischer, P., Brox, T.:  U-net: Convolutional networks for biomedical image segmentation.  In: International Conference on Medical image computing and computer-assisted intervention, Springer (2015) 234–241


# ReDUCE: Reformulation of Mixed Integer Programs using Data from Unsupervised
Clusters for Learning Efficient Strategies

Xuan Lin1, Gabriel I. Fernandez1, and Dennis W. Hong1 1X. Lin, Gabriel I.
Fernandez, and Dennis W. Hong are with the Robotics and Mechanisms Laboratory,
Department of Mechanical and Aerospace Engineering, University of California,
Los Angeles, CA 90095, USA<EMAIL_ADDRESS>
###### Abstract

Mixed integer convex and nonlinear programs, MICP and MINLP, are expressive
but require long solving times. Recent work that combines learning methods on
solver heuristics has shown potential to overcome this issue allowing for
applications on larger scale practical problems. Gathering sufficient training
data to employ these methods still present a challenge since getting data from
traditional solvers are slow and newer learning approaches still require large
amounts of data. In order to scale up and make these hybrid learning
approaches more manageable we propose ReDUCE, a method that exploits structure
within small to medium size datasets. We also introduce the bookshelf
organization problem as an MINLP as a way to measure performance of solvers
with ReDUCE. Results show that existing algorithms with ReDUCE can solve this
problem within a few seconds, a significant improvement over the original
formulation. ReDUCE is demonstrated as a high level planner for a robotic arm
for the bookshelf problem.

##  I Introduction

Optimization-based methods are useful tools for solving robotic motion
planning problems. Typical approaches such as mixed-integer convex programs
(MICPs) [1, 2], nonlinear or nonconvex programs (NLPs) [3, 4, 5] and mixed-
integer NLPs (MINLPs) [6] offer powerful tools to formulate these problems.
However, each has its own drawbacks. NLPs tend to suffer from local optimal
solutions. In practice, local optimal solutions can sometimes have bad
properties, such as inconsistent behavior as they depend on initial guesses.
Mixed-integer programs (MIPs) are a type of NP-hard problem. Branch-and-bound
is usually used to solve MIPs [7]. MIP solvers seek global optimal solutions,
therefore, having more consistent behavior than NLP solvers. For small-scale
problems, these algorithms usually find optimal solutions within a reasonable
time [2, 8]. On the contrary, MIPs can require impractically long solving
times for problems with a large number of integer variables [9]. MINLPs
incorporate both integer variables and nonlinear constraints, hence, very
expressive. Unfortunately, we lack efficient algorithms to tackle MINLPs. Many
practical problems require a solving speed of at most a few seconds. As a
result, it is difficult to implement most of the optimization schemes online
for larger-scale problems.

Figure 1: Angle and x position of solutions to the bookshelf problem’s
furthest left book are depicted by colored numbers indicating cluster
membership. Each box on the grid corresponds to a different relaxation scheme
using integer variables: (a) full relaxation of the whole space, (b)
relaxation over the space with data, (c) relaxation for specific modes from
clusters. For visualization purposes only mode 4 region is explicitly
displayed. Notice how different modes can reduce the amount of integer
variables needed.

Recently, researchers have started to investigate machine learning methods to
gather problem specific heuristics and speed up the MIP solving process.
Standard algorithms to solve MIPs such as branch-and-bound and cutting plane
methods rely on heuristics to quickly remove infeasible regions. Learning
methods can be used to acquire better heuristics. For example, [10] used graph
neural networks to learn heuristics. [11] used reinforcement learning to
discover efficient cutting planes. On the other hand, data can be collected to
learn and solve specific problems [12, 13]. In [13], the authors proposed CoCo
which collects problem features and solved integer variable strategies offline
and trains a neural network. Effectively, this becomes a classification
problem, where each strategy has a unique label. For online solving, the
neural network will propose candidate solutions, reducing it to a convex
program. The results show that CoCo can solve MIPs with around 50 integer
variables from within a second. However, with more integer variables, CoCo
suffers in two aspects. First, the number of possible strategies equals to
$2^{\lambda}$, where $\lambda$ is the number of integer variables. As a result
the number of unique strategies tend to be close to the total amount of data.
This induces overfitting. Second, as the amount of integer variables increase,
the solving speed dramatically slows. As a consequence, it is difficult to
collect enough reliable training data.

In this paper, we propose ReDUCE, an algorithm that combines previous
unsupervised learning work [9] with supervised learning, e.g., CoCo [13], to
solve larger-scale MICPs and MINLPs online. Unsupervised learning is employed
on a small amount of initial data to retrieve sub-regions, clusters, inside
the solution space. Integer variables are then assigned to each cluster. This
allows us to retrieve important regions on the solution manifold and reduce
the amount of integer variables needed. This then allows for fast generation
of much larger datasets to train supervised learners on. All datasets
generated can then be used to train a final learning model which allows
interpolation between clusters. We also introduce the bookshelf organization
problem in this paper to demonstrate ReDUCE’s capabilities. Given a bookshelf
with several books on top, an additional book needs to be placed on the shelf
with minimal disturbance on the existing books. The bookshelf problem works
well as a good benchmark because it: 1) is an MINLP that can be converted to
an MICP problem with hundreds of integer variables, 2) can easily be scaled to
push algorithms to their limit, and 3) has practical significance where data
can reasonably be gathered, such as in the logistics industry.

To summarize, our contributions are as follows:

  1. 1.

Extend supervised learning schemes, e.g., CoCo, to solve MINLPs, where
nonlinear or non-convex constraints are converted into MICPs constraints using
convex envelopes,

  2. 2.

Use unsupervised learning to formulate the mixed-integer envelope constraints,
significantly speeding up the data collection process for problems unsuitable
for MICPs, i.e., very slow solving speeds, and

  3. 3.

Formulate the bookshelf organization problem as an MINLP and solve it within
seconds with ReDUCE.

##  II Bookshelf Organization Problem Setup

Figure 2: Complete formulation of the bookshelf organization problem.

Assume a 2D bookshelf with limited width $W$ and height $H$ contains
rectangular books where book $i$ has width $W_{i}$ and height $H_{i}$ for
$i=1,...,N-1$. A new book, $i=N$, is to be inserted into the shelf. The
bookshelf contains enough books in various orientations, where in order to
insert book $N$, the other $N-1$ books may need to be moved, i.e., optimize
for minimal movement of $N-1$ books. This problem is useful in the logistics
industry, such as robots filling shelves.

Fig. 2 shows the constraints, variables and objective function for the
bookshelf problem. The variables that characterize book $i$ are: position
$\textbf{x}_{i}=[x_{i},y_{i}]$ and angle $\theta_{i}$ about its centroid.
$\theta_{i}=0$ when a book stands upright. The rotation matrix is:
$\textbf{R}_{i}=[cos(\theta_{i}),\ -sin(\theta_{i});\ sin(\theta_{i}),\
cos(\theta_{i})]$. Let the 4 vertices of book $i$ be $\textbf{v}_{i,k}$,
$k=1,2,3,4$. The constraint A in Fig. 2 shows the linear relationship between
$\textbf{x}_{i}$ and $\textbf{v}_{i,k}$, where $\textbf{h}_{i,k}$ is the
constant offset vector from its centroid to vertices. Constraint B enforces
that all vertices of all books stay within the bookshelf, a linear constraint.
Constraint C enforces the orthogonality of the rotation matrix, a bilinear
(non-convex) constraint. Constraint D enforces that the angle $\theta_{i}$
stays within $[-90\degree,90\degree]$, storing books right side up.

To ensure that the final book positions and orientations do not overlap with
each other, separating plane constraints are enforced. For convex shapes, the
two shapes do not overlap with each other if and only if there exists a
separating hyperplane $\textbf{a}^{T}\textbf{x}=b$ in between [14]. That is,
for any point $\textbf{p}_{1}$ inside shape 1 then
$\textbf{a}^{T}\textbf{p}_{1}\leq b$, and for any point $\textbf{p}_{2}$
inside shape 2 then $\textbf{a}^{T}\textbf{p}_{2}\geq b$. This is represented
by constraint E. Constraint F enforces a to be a normal vector. Both E and F
are bilinear constraints.

Finally, we need to assign a state to each book $i$. For each book, it can be
standing straight up, laying down on its left or right, or leaning towards
left or right against some other book, as shown in the far right column in
Fig. 2. For each book $i$, we assign a set of integer variables ${z_{i}}$. If
book $i$ stands upright ($z_{i,2}=1$) or lays flat on its left ($z_{i,1}=1$)
or right ($z_{i,3}=1$), constraints I1 or J1 or H1 are enforced, respectively.
If book $i$ leans against another book on the left or right, constraints in K
and L are enforced, respectively. To this end checks need to indicate the
contact between books. By looking at the right column in Fig. 2, we can
reasonably assume that the separating plane $\textbf{a}^{T}\textbf{x}=b$
always crosses vertex 1 of the book on the left and vertex 4 of the book on
the right. This is represented by bilinear constraints, K1 and L1. In
addition, the books need to remain stable given gravity. Constraints K2 and L2
enforce that a book is stable if its $x$ position stays between the supporting
point of itself (vertex 2 if leaning rightward and vertex 3 if leaning
leftward) and the $x$ position of the book that it is leaning onto. Lastly,
constraint K3 and L3 enforce that the books have contact with the _ground_.
For practical reasons, we assume that books cannot stack onto each other, i.e,
each book has to touch the _ground_ of bookshelf at at least one point. We
note that constraints in H, I, J, K, and L can be easily formulated as MICP
constraints using big-M formulation [15], such that they are enforced only if
the associated integer variable $z=1$. Also, it can easily be extended to
allow stacking for our problem. Any contact conditions between pairs of books
may also be added into this problem as long as it can be formulated as mixed-
integer convex constraint. Overall, this is a problem with integer variables,
${z_{i}}$, and non-convex constraints C, E, F, K1, and L1, hence, an MINLP
problem.

Practically, this problem presents challenges for retrieving high quality
solutions. If robots were used to store books, the permissible solving time is
several seconds, and less optimal solutions means longer realization times.
For example, in Fig. 4 a non-optimal insertion induces multiple additional
robot motions that dramatically increase the chance of failure. There are
several potential approaches to resolving this issue: fix one set of nonlinear
variables and solve MICP [16], convert the nonlinear constraints into piece-
wise linear constraints and formulate them into an MICP [1], or directly
applied MINLP solvers such as BONMIN [17]. As expected, these approaches
struggle to meet the requirements. In this paper, we implement ReDUCE to
satisfy them.

##  III Learning Algorithm

Previous work demonstrated the potential of learning mixed-integer strategies
offline and then using the learned model to sample candidate solutions and
solve convex programs online [12, 13, 18]. In [9], the authors proposed an
unsupervised learning method to identify important regions in the solution
space. This approach effectively reformulated the whole MIP into multiple
problems with a reduced number of integer variables. ReDUCE further builds
upon this notion and combines those two approaches.

Assume that we are given a set of problems parametrized by $\Theta$ that is
drawn from a distribution $D(\Theta)$. For each $\Theta$, we seek a solution
to the optimization problem:

$\displaystyle\underset{\textbf{x},\ \textbf{z}}{\text{minimize}}\
f_{obj}(\textbf{x},\textbf{z};\Theta)$ (1) s. t. $\displaystyle
f_{i}(\textbf{x},\textbf{z};\Theta)\leq 0,\ \ i=1,...,m_{f}$ $\displaystyle
b_{j}(\textbf{x},\textbf{z};\Theta)\leq 0,\ \ j=1,...,m_{b}$

Where x denotes continuous variables and z binary variables with
$z_{i}\in\\{0,1\\}$ for $i=1,...,dim(\textbf{z})$. Constraints $f_{i}$ are
mixed-integer convex, meaning if the binary variables z are relaxed into
continuous variables $\textbf{z}\in[0,1]$, $f_{i}$ becomes convex. Constraints
$b_{j}$ are mixed-integer bilinear, meaning that relaxing the binary variables
gives bilinear constraints. Without loss of generality, x and z are assumed to
be involved in each constraint. We omit equality constraints in (1) as they
can be turned into two inequality constraints from opposite directions.
Similar to [9], the solution set with respect to fixed z,
$X(\Theta_{i};\textbf{z})$, is defined to be the manifold containing all
feasible x given z. If z is infeasible, $X(\Theta_{i};\textbf{z})=\emptyset$.
The full solution set, $X(\Theta_{i})$, for the problem embedded inside the
solution space, $\mathbb{S}\in\mathbb{R}^{dim(\textbf{x})}$, is defined by
$\cup_{\textbf{z}}X(\Theta_{i};\textbf{z})\;\forall\,\textbf{z}$.

In this paper, we take the approach to convert bilinear constraints into
mixed-integer linear constraints by gridding the solution space $\mathbb{S}$
and approximating the constraints locally inside grids with McCormick
envelopes similar to [19]. A McCormick envelope relaxation of one bilinear
constraint $w=xy$ [20] is the best linear approximation defined over a pair of
lower and upper bounds $[x^{L},x^{U}]$ and $[y^{L},y^{U}]$. Therefore, we
first assign grids $G(\textbf{x})$ to $\mathbb{S}$. Let $\\{g_{l}\\}$,
$l=1,...,L$ be one cell in the grid with upper and lower bounds
$[\textbf{x}^{L},\textbf{x}^{U}]$. We introduce additional integer variables
n, where $n_{i}\in\\{0,1\\}$. Each unique value of n corresponds to one cell
in the grid within which McCormick envelope relaxations are applied. The
constraints $b_{j}(\textbf{x},\textbf{z};\Theta)\leq 0$, $j=1,...,m_{b}$ are
converted into $L$ constraints:
$E_{b_{j},l}(\textbf{x},\textbf{z},\textbf{n};\Theta)\leq 0$ for
$j=1,...,m_{b}$ and $l=1,...,L$, turning (1) into an MICP:

$\begin{split}&\underset{\textbf{x},\ \textbf{z},\
\textbf{n}}{\text{minimize}}\
f_{obj}(\textbf{x},\textbf{z},\textbf{n};\Theta)\\\ \text{s.
t.}\qquad\;f_{i}(\textbf{x},\textbf{z};\Theta)\leq 0,&\ \ i=1,...,m_{f}\\\
E_{b_{j},l}(\textbf{x},\textbf{z},\textbf{n};\Theta)&\leq 0,\ \
j=1,...,m_{b},\ \ l=1,...,L\qquad\end{split}$ (2)

To speed up the solving speed for MICPs, [18, 13, 21] define integer
strategies to be tuples of
$\mathcal{I}(\Theta_{i})=(\textbf{z}^{*},\textbf{n}^{*},\mathcal{T}(\Theta_{i}))$,
where $(\textbf{x}^{*},\textbf{z}^{*},\textbf{n}^{*})$ is an optimizer for
problem (2), $\mathcal{T}(\Theta_{i})=\\{i\in{1,...,m_{f}},\ l\in
1,...,L|f_{i}(\textbf{x}^{*},\textbf{z}^{*};\Theta)=0,\
E_{b_{j},l}(\textbf{x}^{*},\textbf{z}^{*},\textbf{n}^{*};\Theta)=0\\}$ is the
set of active inequality constraints. Given an optimal integer strategy
$\mathcal{I}(\Theta_{i})$, solution to (2) can be retrieved through solving a
convex optimization problem. If this approach returns a small amount of
integer variables, it may perform well. However, approximating nonlinear
constraints with mixed-integer convex constraints usually gives large number
of integer variables when high approximation accuracy is desired [19]. As a
result, the solving time grows fast. This makes it difficult to collect data
to train supervised learners.

This paper uses clustering method on a relatively smaller amount of pre-solved
data to identify important regions on $X(\Theta_{i})$, and implement McCormick
envelope relaxations around those regions which significantly improves the MIP
solving speed. The nominal dimension of the solution space is
$dim(\textbf{x})$. However, due to the existence of constraints, the actual
solution set of $X(\Theta)$ is of much lower dimension. Relaxing the complete
space of x is unnecessary as many grid cells are infeasible or non-optimal and
the MIP solver should not need to spend time in exploring those regions.
Furthermore, we almost never need the complete $X(\Theta)$ manifold in
practice. Many robotic systems operate under different modes corresponding to
regions on $X(\Theta)$ where the optimal solutions $\textbf{x}^{*}$ populate.
With data unsupervised learning can identify those regions. Within each
region, we only need a smaller amount of integer variables
$\textbf{z},\textbf{n}$. Consequently, the solving speed can be significantly
improved.

The detailed steps of ReDUCE are shown in Algorithm 1. To begin, ReDUCE
requires a relatively small amount of pre-solved data
$(\Theta_{k},\textbf{x}_{k})$, $k=1,...,K$ termed kick off data throughout.
This data may be collected by solving the non-reduced formulation (2). If (2)
represents a practical problem, the data may come from simulation or human
demonstration on real hardware. We also need to pre-assign grids
$G(\textbf{x})$ to the solution space $\mathbb{S}$. The size of grids depends
on the approximation accuracy requirement for bilinear constraints. ReDUCE
begins by performing unsupervised learning on the kick off data to retrieve
clusters that indicates regions on $X(\theta)$. Density-Based Spatial
Clustering of Applications with Noise (DBSCAN) [22] was used to cluster on x
which gives clusters ${\textbf{x}_{1}},{\textbf{x}_{2}},...,{\textbf{x}_{C}}$,
where $C$ is the number of clusters. We then trace
$\Theta\rightarrow\textbf{x}$ map backwards to create clusters in $\Theta$
space, i.e.
$(\\{\Theta\\}_{1},\\{\textbf{x}\\}_{1}),...,(\\{\Theta\\}_{C},\\{\textbf{x}\\}_{C})$.
Next, a supervised classifier is trained to classify $(\Theta,c)$. We use
random forest which requires relatively smaller amounts of data to train than
deep learning methods. Thus, the amount of _kick off_ data, $K$, can be
relatively small.

The main difficulty to generate training data for (2) is its slow solving
speed due to large amount of integer variables. If
$\mathbb{R}^{dim(\textbf{x})}$ is segmented into smaller regions, e.g.
clusters, the required integer variables for each cluster can be reduced. This
may be seen in Fig. 1 which is an instance of 2-dimensional (dim) $X(\Theta)$
from a bookshelf experiment described in Sec. IV-A. In this paper, we use a
$log_{2}N$ formulation which means $N$ grids will be represented by $\lceil
log_{2}N\rceil$ integer variables. For example, $17\sim 32$ grids can be
represented by 5 integer variables. In Fig. 1, the complete space has 27 grids
which can be represented with at least 5 integer variables. However, mode 4
occupies only 6 grids which is represented with 3 integer variables. One
advantage of using DBSCAN is that it uses a threshold to decide the boundaries
of each cluster and identify outliers. Since DBSCAN is based off of densities,
it is able to predict outliers from the data. Outliers can be removed from
training data wherein the classifier can be used to make a prediction. Based
on that prediction, it is possible that some outliers may have membership to
multiple clusters. Also outliers may indicate insufficient sample size near
that particular region. This may guide practitioners where to collect more
data.

With the trained classifier, we can classify a much greater amount of problems
$\Theta\sim D(\Theta)$ and solve them quickly within each clustered region to
collect more training data, indicated in line 4-14 of Algorithm 1. For all
$\\{\textbf{x}\\}_{c}$ in cluster $c$, we find the grid cells that they occupy
and re-assign integer variables $\textbf{z}_{c}$ and $\textbf{n}_{c}$ to each
cluster (line 8-10). Eventually all data are put together for training, so the
definition of integer variables needs to be consistent across clusters. Thus,
we recover the original integer variables z and n (line 13). With $\Theta$,
$\textbf{z}_{c}$ and $\textbf{n}_{c}$ for cluster $c$, problem (2) can be
reformulated and solved. The feature $\Theta$, solution $\textbf{x}^{*}$,
recovered integer variables $\textbf{z}^{*}$ and $\textbf{n}^{*}$ are added to
the dataset. Finally, the dataset is used to train a strategy learner, e.g.
CoCo.

The idea of solving smaller sub-MIPs bears similarities with algorithms such
as RENS [23] and Neural Diving [10] where a subset of integer variables are
fixed from linear program solutions or a learned model, while the others are
solved. Our method, however, segments the problem through an unsupervised
clustering approach with the goal of generating enough training data for
supervised learning.

Algorithm 1 $\operatorname{ReDUCE}$

Input |  Kick off data $\\{(\Theta_{k},\;\textbf{x}_{k})\\},\;k=1,...,K$,
---|---
|  Grid $G(\textbf{x})\;\text{with cells}\;\\{g_{l}\\},\;l=1,...,L$,
|  DBSCAN clustering threshold  $\epsilon$
|  $d_{c}$ desired amount of samples on label $c,\;c=1,...,C$

1:
DBSCAN$(\\{\textbf{x}\\},\epsilon)\longrightarrow\\{\textbf{x}\\}_{1},...,\\{\textbf{x}\\}_{C}$

2:
$\\{\\{\Theta\\}_{c}:\textbf{x}\in(\Theta,\textbf{x})\land\\{\textbf{x}\\}_{c}\;\forall\;c=1,...,C\\}$

3: Train Random Forest: RF$(\\{\Theta\\})\rightarrow\\{c\\}$

4: while Samples $|\\{\Theta\\}_{c}|<d_{c}$ do

5: Sample $\Theta\sim D(\Theta)$, Classify RF$(\Theta)$ and Store $(\Theta,c)$

6: for $c=1,...,C$ do

7: Initialize grid cell list $\\{g_{c}\\}$ and dataset $S_{c}$

8: for each point in $\\{\textbf{x}\\}_{c}$ do

9: Find corresponding grid $g_{l}$ and add $g_{l}$ to $\\{g_{c}\\}$

10: Assign integer variables $\textbf{z}_{c},\textbf{n}_{c}$ to $\\{g_{c}\\}$

11: Formulate $\mathcal{P}$ from (2) with $\\{\Theta\\}_{c}$, replacing z, n
by $\textbf{z}_{c},\textbf{n}_{c}$

12: if $\mathcal{P}$ has solution
$\textbf{z}_{c}^{*},\textbf{n}_{c}^{*},\textbf{x}^{*}$ then

13: Recover original $\textbf{z}^{*}$, $\textbf{n}^{*}$ from
$\textbf{z}_{c}^{*},\textbf{n}_{c}^{*}$

14: Add $\Theta$, $\textbf{z}^{*}$, $\textbf{n}^{*}$, $\textbf{x}^{*}$ to
$S_{c}$

15: Initialize Neural Network or desired model: $f(\cdot)$

16: Use $S_{c}$ to train:
$f(\Theta)\rightarrow(\textbf{z}^{*},\textbf{n}^{*})$

17: return $f(\cdot)$

##  IV Experiment

###  IV-A Experiment Setup

We use ReDUCE to solve the bookshelf organization problem. We place 3 books
inside the shelf where 1 additional book is to be inserted. Grids are assigned
to the variables involved in the non-convex constraint C, E, F, K1, L1:
$\textbf{R}_{i}(\theta_{i})$, $\textbf{a}_{j}$ and $\textbf{v}_{i,k}$. These
variables span a 48 dimensional space. The rotation angles $\theta_{i}$, which
includes $\textbf{R}_{i}$, are gridded at a $\frac{\pi}{8}$ interval. Elements
in $\textbf{a}_{j}$ are gridded on 0.25 intervals. Elements in
$\textbf{v}_{i,k}$ are gridded at intervals $\frac{1}{4}$ the shelf width $W$
and height $H$. Since books are not allowed to stack on each other, there is
an order of books from left to right. We order the feature and solution vector
according to this order. We use an MICP formulation that has $log_{2}N$
integer variables (explained in detail in the appendix) where $N$ is the
number of grids. This results in 130 integer variables in total. The feature
vector includes the centroid positions, angles, heights and widths of stored
books and height and width of the book to be inserted. The feature dimension
is 17.

The _kick off_ data was collected using a 2-dim simulated environment of books
on a shelf. Initially, 4 randomly sized books are arbitrarily placed on the
shelf, and then 1 is randomly removed and regarded as the book to be inserted.
Contrary to the sequence, the initial state with 4 books represents one
feasible (not necessarily optimal) solution to the problem of placing a book
on a shelf with 3 existing books. Since this problem can be viewed as high-
level planning for robotic systems, the simulated data is sufficient. For
applications outside of the scope of this paper real world data may be
preferable in this pipeline.

###  IV-B Unsupervised Learning

| ReDUCE | CoCo+Re | RF+Re | Baseline+Re | MIP+Re
---|---|---|---|---|---
| Cluster
---
blue!10 0
blue!10
blue!10
1
blue!10 2
blue!10
blue!10
3
| N/total | % | Int
---|---|---
blue!10 997/1000 |  blue!10 99.7% |  blue!1077
blue!10 2984/2999 |  blue!10 99.4% |  blue!1077
blue!10 4947/5000 |  blue!10 98.9% |  blue!1077
698/699 |  99.9% | 77
1996/1999 |  99.8% | 77
5976/5999 |  99.6% | 77
blue!10 599/600 |  blue!10 99.8% |  blue!1066
blue!10 2970/2999 |  blue!10 99.0% |  blue!1066
blue!10 5476/5599 |  blue!10 97.8% |  blue!1066
339/399 |  85.0% | 46
681/900 |  75.6% | 46
| S% | Det | Avg | Max
---|---|---|---
blue!10 93.4% |  blue!10 210 |  blue!10 1.39 |  blue!10 7.39
blue!10 97.0% |  blue!10 144 |  blue!10 1.25 |  blue!10 7.35
blue!10 98.0% |  blue!10 119 |  blue!10 1.05 |  blue!10 7.23
91.2% | 148 | 1.42 | 7.37
96.2% | 93 | 1.30 | 7.58
99.0% | 73 | 1.24 | 7.12
blue!10 91.0% |  blue!10 95 |  blue!10 1.36 |  blue!10 6.70
blue!10 96.8% |  blue!10 60 |  blue!10 1.13 |  blue!10 7.37
blue!10 98.4% |  blue!10 47 |  blue!10 1.07 |  blue!10 7.47
94.0% | 28 | 1.18 | 7.03
97.2% | 20 | 0.95 | 6.92
| S% | Det | Avg | Max
---|---|---|---
blue!10 93.2% |  blue!10 348 |  blue!10 1.42 |  blue!10 7.58
blue!10 97.0% |  blue!10 253 |  blue!10 1.31 |  blue!10 7.46
blue!10 96.0% |  blue!10 194 |  blue!10 1.17 |  blue!10 6.75
90.0% | 267 | 1.35 | 7.10
92.6% | 204 | 1.27 | 7.21
96.8% | 169 | 1.12 | 7.57
blue!10 90.6% |  blue!10 122 |  blue!10 1.33 |  blue!10 7.51
blue!10 93.8% |  blue!10 88 |  blue!10 1.00 |  blue!10 7.42
blue!10 95.2% |  blue!10 84 |  blue!10 1.04 |  blue!10 6.97
92.6% | 37 | 0.94 | 7.04
97.8% | 29 | 0.73 | 6.94
| S% | Det
---|---
blue!10 65.4% |  blue!10 856
blue!10 63.7% |  blue!10 862
blue!10 65.5% |  blue!10 933
57.4% | 608
60.0% | 554
60.5% | 543
blue!10 59.5% |  blue!10 281
blue!10 57.0% |  blue!10 293
blue!10 57.7% |  blue!10 251
63.2% | 162
57.5% | 132
| Avg | Max
---|---
blue!10 4.52 |  blue!10 37.86
blue!10 4.52 |  blue!10 37.86
blue!10 4.52 |  blue!10 37.86
2.27 | 16.61
2.27 | 16.61
2.27 | 16.61
blue!10 0.92 |  blue!10 5.46
blue!10 0.92 |  blue!10 5.46
blue!10 0.92 |  blue!10 5.46
0.21 | 0.78
0.21 | 0.78
TABLE I: Comparison of different algorithms with ReDUCE (+Re) for CoCo
(CoCo+Re), Random Forest (RF+Re), random sampling (Baseline+Re), and MIP
(MIP+Re). Success rate and deterioration on the objective is denoted as S% and
Det, respectively. Average solving time and maximum solving time are denoted
as Avg and Max, respectively. Int stands for the number of integer variables
within a cluster. With more data the algorithms generally tend to improve
performance.

We randomly sample 4,000 bookshelves and implement DBSCAN to realize 100
clusters. Fig. 1 shows the first 2 dimensions of the projected solution set
$X(\Theta)$ and 6 modes with distinct labels and colors. Fig. 3 shows the
solution set packed into tighter groups compared to the same sample of
features using t-distributed Stochastic Neighbor Embedding (t-SNE) [24] as a
projection of the high dimensional space to 2-dim. The upper graph in Fig. 3
shows the clusters in the solution space, while the lower one depicts the
corresponding clusters in the feature space according to $(\Theta,\textbf{x})$
from the kick off data. The colors in Fig. 3 denote the different clusters. We
can tell obvious separations in the solution space. Although the clusters are
more intertwined in the feature space because of the complexity of
$\Theta\longrightarrow\textbf{x}$ mapping, there are still distinct regions
where certain colors are more dominant. We trained a random forest (RF)
classifier on the features $\Theta$. The classification accuracy reaches 97%
indicating that RF$(\\{\Theta\\})\longrightarrow\\{\textbf{x}\\}$ is able to
achieve a reasonable mapping.

Figure 3: _Top_ : projection onto a 2-dim manifold using t-SNE to depict
clustering of a 48-dim solution space. Each color and label specify a solution
clustering. Under this projection the solution appears very much structured.
_Bottom_ : projection of the 17-dim feature space corresponding to the
solution using t-SNE. The clusters have some structure over the feature space
with certain labels only being found in certain regions. For visualization
purposes only the top 8 clusters out of 100 are being displayed.

###  IV-C Supervised Learning

| Success (%) | Det | Avg Time | Max Time
---|---|---|---|---
CoCo+Re | 99.2% | 140 | 1.21 sec | 13.47 sec
MIP+Re | 100 % | 0 | 2.36 sec | 48 sec
MIP | n/a | n/a |  $>$ 851 sec | n/a
MINLP | n/a | n/a |  $>$ 10 min | n/a
TABLE II: CoCo+Re (+Re denotes using ReDUCE) returned the top 50 candidate
solutions over 100 clusters. MIP+Re used those same clusters to solve the
problem on Gurobi. MIP also used Gurobi. MINLP ran using BONMIN as the solver.
MIP and MINLP solving time for the full set of samples have exceeded
reasonable limitations beyond practical purposes.  Figure 4: Solved scenes of
bookshelves. _Top Row_ : Instance solved by CoCo with ReDUCE. _Bottom Row_ :
Instance solved by RF with ReDUCE. Column (a), (b), (c) corresponds to cluster
0, 1, 2, respectively. Each cell contains a before and after scene with the
upper diagram showing the original bookshelf with orange rectangles
representing stored books and the lower diagram showing the solved bookshelf
where the blue rectangle represents the inserted book. For column (a), we
demonstrate a scene where RF gives a much worse result than CoCo such that
multiple books are moved. However, CoCo and RF usually give similar results.

With the RF classifier trained on 100 clusters, we can quickly sample
different bookshelf problems and solve the reduced MIP problem within seconds
to collect more data for supervised learning. To verify the benefits with
ReDUCE, we first run an experiment that observes an increase in performance
with increasing amounts of data over several clusters each containing
different number of integer variables (denoted as Int in Table I). We pick 3
different clusters with number of integer variables 77, 66, and 46,
respectively. Four different methods to solve the reduced problem within a
cluster are used and tested on a fixed testing set of 500 data points. The
results are shown in Table I. The column titled ReDUCE shows the total number
number of unique integer strategies (N) over training data (total), and the
uniqueness percentage (%) is the resulting quotient. For the column titled
CoCo+Re, we train a neural network with one hidden layer of size 10,000. The
input size is equal to dimension of feature space (17) and the output layer is
equal to the number of unique strategies within training data (N). The neural
network then samples 30 candidate strategies online according to the rank of
the softmax scores of the network and solves the associated convex
optimization problems. If one strategy gives a feasible solution, we terminate
the process and record the solving time and optimal cost. Otherwise,
infeasibility is recorded. Since the maximum number of convex problems
considered is 30 to be solved online, the problem setup time is non-
negligible. To avoid additional overhead we only setup the problem once and
keep on modifying the integer constraints for more instances. Thus, all
solving times include the problem setup time. The solving process is done on a
Core i7-7800X 3.5GHz $\times$ 12 machine with Gurobi. For column RF+Re, we
instead train a random forest with 150 decision trees and get the top 30 most
voted strategies for candidate solutions. For column titled Baseline+Re, we
simply randomly sample 30 unique strategies from the training strategies. For
column MIP+Re, we solve the MIP with reduced number of integer variables,
instead of using supervised learning or sampling method. For all sampling
methods, we record the rate of feasibility (S%), the deterioration (Det) of
optimal cost compared to the optimal cost from MIP+Re column. That is, for the
MIP+Re column, its feasibility rates are all 100%, and its optimal
deteriorations are 0. Solving times for baseline is omitted as they are
significantly longer than learning based methods. All times under Avg (average
solving time) and Max (maximum solving time) columns are given in seconds.

As the results show, the performances of learning based methods in general
improve with more data. As ReDUCE improves the solving speed to collect larger
amounts of data, we can further increase the performances with more data. The
random sampling baseline (Baseline+Re) has significantly worse performance
than learning methods. This column ensures that the clusters are not too small
making the reduced problem too simple. For larger clusters with more integer
variables (e.g. cluster 0), the learning methods demonstrate an increase in
solving speed over MIP. For smaller clusters (e.g. cluster 2), MIP has faster
solving speeds. There is a tradeoff between problem scale and solving methods
depending on the cluster size.

For a second set of experiments, we collect data for all 100 clusters, in
total 17,000, and train a neural network with the same size as above. The
network is then tested on a 1,000 testing set with 100 clusters blended
together. The network then samples 50 strategies. The percentage of
feasibility (Success (%)), the optimal cost deterioration (Det), the average
(Avg Time) and max (Max Time) solving time are shown in Table II. For
comparison, we record the average and maximal solving time for MIP with ReDUCE
(all clusters blended together). We also record the solving time for MIP
without ReDUCE, the original formulation, and the solving time for formulation
(1) solved with BONMIN, an MINLP solver. CoCo solves faster than MIP for
larger clusters which is associated with more data from the kick off data.
When sampling $D(\Theta)$, we get more samples for larger clusters. Therefore,
solving time is improved when averaged. For non-reduced MIP, Table II shows
the average solving time over 10 samples where the solving process is
interrupted if it exceeds 1,000 sec. Thus, the average solving time is at
least 851 sec. It is clear that without ReDUCE, it is intractable to gather
the amount of data required to train a learner of decent performance. Similar
speed is seen for the MINLP solver. All the code is available at:
https://github.com/RoMeLaUCLA/ReDUCE.git.

###  IV-D Hardware Experiment

We implement our optimization results on hardware with a 6 degree of freedom
manipulator [25] to insert a book onto a shelf. To automate the system, a
trajectory planner is required to generate the insertion motion, which is
common in practice. As the focus of this paper is on high level planning, we
simplify this part by manually selecting waypoints. The hardware
implementation is shown in the attached video.

##  V Conclusion, Discussion and Future Work

This paper proposes ReDUCE, an algorithm that relies on unsupervised learning
on relatively smaller dataset to reduce the number of integer variables and
generate large amounts of strategies for supervised learners. The results
demonstrated improvements in solving speed over the original MIP formulations
allowing for practical implementations.

As this algorithm requires kick off data, considering data collection is an
important aspect. Beyond collecting data from simulations or hardware
demonstrations, we consider automatically generating formulations for known
problems and transfer knowledge to a target domain as an interesting next
step. In real life, we often encounter problems with similar solutions in
different contexts. Developing tools that automatically generate MIP
formulations from other domains may help generate kick off data.

It is noticed from the test results that higher feasibility rates sometime
lead to lower optimality. This may be due to the optimal strategy for the
testing problem not included in the training data. One future work is using
generative learning approaches to combine features of different training
instances to generate new strategies. Other future work include increasing the
amount of data to further boost performance, increasing the amount of books to
test the limit of ReDUCE, implementing parallel computing to increase online
solving speed, and comparing the performance of ReDUCE with reinforcement
learning algorithms.

## Acknowledgements

The authors would like to thank Abhishek Cauligi and Professor Bartolomeo
Stellato for helpful discussions, and Yusuke Tanaka and Hyunwoo Nam for the
assisting with hardware implementation.

## References

  * [1] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in _2014 IEEE-RAS international conference on humanoid robots_. IEEE, 2014, pp. 279–286.
  * [2] X. Lin, J. Zhang, J. Shen, G. Fernandez, and D. W. Hong, “Optimization based motion planning for multi-limbed vertical climbing robots,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 1918–1925.
  * [3] H. Dai, A. Valenzuela, and R. Tedrake, “Whole-body motion planning with centroidal dynamics and full kinematics,” in _2014 IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2014, pp. 295–302.
  * [4] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, 2018.
  * [5] Y. Shirai, X. Lin, Y. Tanaka, A. Mehta, and D. Hong, “Risk-aware motion planning for a limbed robot with stochastic gripping forces using nonlinear programming,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 4994–5001, 2020.
  * [6] M. Soler, A. Olivares, P. Bonami, and E. Staffetti, “En-route optimal flight planning constrained to pass through waypoints using minlp,” 2011.
  * [7] S. Boyd and J. Mattingley, “Branch and bound methods,” _Notes for EE364b, Stanford University_ , pp. 2006–07, 2007.
  * [8] J. Tordesillas, B. T. Lopez, and J. P. How, “Faster: Fast and safe trajectory planner for flights in unknown environments,” in _2019 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2019, pp. 1934–1940.
  * [9] X. Lin, M. S. Ahn, and D. Hong, “Designing multi-stage coupled convex programming with data-driven mccormick envelope relaxations for motion planning,” in _2021 IEEE/RSJ International Conference on Robotics and Automation (ICRA)_. IEEE, 2021.
  * [10] V. Nair, S. Bartunov, F. Gimeno, I. von Glehn, P. Lichocki, I. Lobov, B. O’Donoghue, N. Sonnerat, C. Tjandraatmadja, P. Wang _et al._ , “Solving mixed integer programs using neural networks,” _arXiv preprint arXiv:2012.13349_ , 2020.
  * [11] Y. Tang, S. Agrawal, and Y. Faenza, “Reinforcement learning for integer programming: Learning to cut,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 9367–9376.
  * [12] J.-J. Zhu and G. Martius, “Fast non-parametric learning to accelerate mixed-integer programming for online hybrid model predictive control,” _arXiv preprint arXiv:1911.09214_ , 2019.
  * [13] A. Cauligi, P. Culbertson, E. Schmerling, M. Schwager, B. Stellato, and M. Pavone, “Coco: Online mixed-integer control via supervised learning,” _arXiv preprint arXiv:2107.08143_ , 2021.
  * [14] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [15] J. P. Vielma, “Mixed integer linear programming formulation techniques,” _Siam Review_ , vol. 57, no. 1, pp. 3–57, 2015.
  * [16] M. Posa, M. Tobenkin, and R. Tedrake, “Stability analysis and control of rigid-body systems with impacts and friction,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 6, pp. 1423–1437, 2015.
  * [17] P. Bonami and J. Lee, “Bonmin user’s manual,” _Numer Math_ , vol. 4, pp. 1–32, 2007.
  * [18] D. Bertsimas and B. Stellato, “The voice of optimization,” _Machine Learning_ , vol. 110, no. 2, pp. 249–277, 2021.
  * [19] H. Dai, G. Izatt, and R. Tedrake, “Global inverse kinematics via mixed-integer convex optimization,” _The International Journal of Robotics Research_ , vol. 38, no. 12-13, pp. 1420–1441, 2019.
  * [20] P. M. Castro, “Tightening piecewise mccormick relaxations for bilinear problems,” _Computers & Chemical Engineering_, vol. 72, pp. 300–311, 2015\.
  * [21] D. Bertsimas and B. Stellato, “Online mixed-integer optimization in milliseconds,” _arXiv preprint arXiv:1907.02206_ , 2019.
  * [22] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “Density-based spatial clustering of applications with noise,” in _Int. Conf. Knowledge Discovery and Data Mining_ , vol. 240, 1996, p. 6.
  * [23] T. Berthold, “Rens,” _Mathematical Programming Computation_ , vol. 6, no. 1, pp. 33–54, 2014.
  * [24] L. Van der Maaten and G. Hinton, “Visualizing data using t-sne.” _Journal of machine learning research_ , vol. 9, no. 11, 2008.
  * [25] D. Noh, Y. Liu, F. Rafeedi, H. Nam, K. Gillespie, J.-s. Yi, T. Zhu, Q. Xu, and D. Hong, “Minimal degree of freedom dual-arm manipulation platform with coupling body joint for diverse cooking tasks,” in _2020 17th International Conference on Ubiquitous Robots (UR)_. IEEE, 2020, pp. 225–232.
  * [26] P. Belotti, L. Liberti, A. Lodi, G. Nannicini, A. Tramontani _et al._ , “Disjunctive inequalities: applications and extensions,” _Wiley Encyclopedia of Operations Research and Management Science_ , vol. 2, pp. 1441–1450, 2011.
  * [27] J. P. Vielma and G. L. Nemhauser, “Modeling disjunctive constraints with a logarithmic number of binary variables and constraints,” _Mathematical Programming_ , vol. 128, no. 1, pp. 49–72, 2011.
  * [28] D. Avis and K. Fukuda, “A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra,” _Discrete & Computational Geometry_, vol. 8, no. 3, pp. 295–313, 1992.

We explain the formulation that is used to solve the MIPs in Sec. IV-A.

For a non-convex constraint, we segment it into multiple regions and locally
approximate or relax them into convex constraints. In this paper, we relax the
bilinear constraints locally into convex polytopes (McCormick envelopes). Each
polytope is associated with a unique combination of integer variable values,
hence, mixed-integer convex constraints. Assume the number of regions used is
$N$. Depending on the number of integer variables used, the formulation can
generally be divided into 2 categories: 1) If the number of integer variables
is $N$, we call it $N$ formulation, e.g., the convex hull formulation [26]. 2)
If the number of integer variables is $log_{2}N$, we call it $log_{2}N$
formulation, e.g., [27]. [27] presents a $log_{2}N$ formulation to model the
special ordered sets of type 2 (sos2). However, there are several limitations
of this formulation. For example, the segmented regions need to be connected
to be a valid sos2 constraint, i.e., the two consecutive couple of non-zero
entries can be any consecutive couple in the set. In this paper, we use MIP to
model convex polytopes of arbitrary locations. This can increase the
complexity for sos2 techniques as the polytopes can be disjunctive. The
disjunctive constraints can be handled with convex hull formulations at a
price of introducing more integer variables which may result in slower solving
speeds.

In this appendix, we demonstrate an intuitive but general $log_{2}N$
formulation to model combinations of convex polytopes at any locations which
serves as the base MIP formulation for the bilinear constraints in our paper.
Assume the variable x is enforced to be within one of the $N$ convex
polytopes, denoted by $\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$,
$i=1,...,N$. We introduce $m=log_{2}N$ binary variables $z_{1},...,z_{m}$,
$z_{i}\in\\{0,1\\}$. Each combination of unique values of binary variables can
be assigned to a convex polytope. Let the assignment be:

$\textbf{z}=\bar{\textbf{z}}_{i}\Rightarrow\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$
(3)

Where $\bar{\textbf{z}}_{i}=[\bar{z}_{i,1},...,\bar{z}_{i,m}]$ are constant
binary values associated with polytope $i$. Note $\bar{z}_{i}\neq\bar{z}_{j}$
if $i\neq j$. In other words, we require that when
$\textbf{z}=\bar{\textbf{z}}_{i}$, x stays within the polytope
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$; otherwise, the constraint is
unenforced.

Denote the vertices of the polytopes by
$\textbf{v}_{i,1},...,\textbf{v}_{i,n_{i}}$, $i=1,...,N$, where $n_{i}$ is the
number of vertices associated with polytope $i$. Each vertex
$\textbf{v}_{i,j}$ is assigned a continuous non-negative variable
$\lambda_{i,j}\in[0,1]$. In general, one can run a mathematical program (e.g.
[28]) to get vertices from the nondegenerate system of inequalities
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$. As a result, the assignment
becomes:

$\textbf{z}=\bar{\textbf{z}}_{i}\Rightarrow\ \ \begin{aligned}
&\textbf{x}=\sum_{j=1}^{n_{i}}\lambda_{i,j}\textbf{v}_{i,j}\\\
&\sum_{j=1}^{n_{i}}\lambda_{i,j}=1,\ \ \lambda_{i,j}\in[0,1]\end{aligned}$ (4)

The formulation can be written as:

$\displaystyle(5.a)\quad\textbf{x}=\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\lambda_{i,j}\textbf{v}_{i,j}$
(5) $\displaystyle(5.b)\quad\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\lambda_{i,j}=1,\
\ \lambda_{i,j}\in[0,1]$ $\displaystyle(5.c)\quad\sum_{k=1,...,N}^{k\neq
i}\sum_{j=1}^{n_{k}}\lambda_{k,j}\leq\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|$

The set of constraint in (5) enforces (4) for all polytopes, as
$\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|=0$ only when
$\textbf{z}=\bar{\textbf{z}}_{i}$, enforcing that all $\lambda$’s that are not
associated with polytope $i$ to be zero. If
$\textbf{z}\neq\bar{\textbf{z}}_{i}$, $\sum_{l=1}^{m}|z_{l}-\bar{z}_{i,l}|\geq
1$ and constraint (5.c) is looser than constraint (5.b), hence, trivial.

Note that formulation (5) works for any convex polytope that can be written as
$\textbf{A}_{i}\textbf{x}\leq\textbf{b}_{i}$. In this paper, the polytopes are
McCormick envelope constraints which is a special case.