# QAL-BP: An Augmented Lagrangian Quantum Approach for Bin Packing Problem

Lorenzo Cellini  Antonio Macaluso  Agents and Simulated Reality Department,
German Research Center for Artificial Intelligence (DFKI), Saarbruecken,
Germany  Michele Lombardi  Department of Computer Science and Engineering,
University of Bologna, Bologna, Italy

###### Abstract

The bin packing is a well-known NP-Hard problem in the domain of artificial
intelligence, posing significant challenges in finding efficient solutions.
Conversely, recent advancements in quantum technologies have shown promising
potential for achieving substantial computational speedup, particularly in
certain problem classes, such as combinatorial optimization.

In this study, we introduce QAL-BP, a novel Quadratic Unconstrained Binary
Optimization (QUBO) formulation designed specifically for bin packing and
suitable for quantum computation. QAL-BP utilizes the augmented Lagrangian
method to incorporate the bin packing constraints into the objective function
while also facilitating an analytical estimation of heuristic, but empirically
robust, penalty multipliers. This approach leads to a more versatile and
generalizable model that eliminates the need for empirically calculating
instance-dependent Lagrangian coefficients, a requirement commonly encountered
in alternative QUBO formulations for similar problems. To assess the
effectiveness of our proposed approach, we conduct experiments on a set of
bin-packing instances using a real Quantum Annealing device. Additionally, we
compare the results with those obtained from two different classical solvers,
namely simulated annealing and Gurobi. The experimental findings not only
confirm the correctness of the proposed formulation but also demonstrate the
potential of quantum computation in effectively solving the bin-packing
problem, particularly as more reliable quantum technology becomes available.

## Introduction

Bin packing is a well-established [1] combinatorial optimization problem with
wide-ranging applications in domains such as logistics, resource allocation,
and scheduling. Its primary objective is to minimize the number of fixed-
capacity bins required to pack a set of items, each with a distinct size.

Despite extensive research efforts, bin packing remains a challenging problem
due to the exponential growth of solution possibilities as the number of items
and bins increases [2]. On the other hand, Quantum Computing has recently
emerged as a promising alternative to solving various AI problems, including
coalition formation in multi-agent systems [3, 4, 5] and supervised learning
[6, 7], though a real practical quantum advantage has yet to be found
considering near-term quantum technology. The standard approach in quantum
computing involves reformulating the original problem as a Quadratic
Unconstrained Binary Optimization (QUBO) problem and employing quantum
annealers (QAs) or parametrized quantum circuits, such as QAOA[8], to find the
optimal solution. These approaches possess distinctive strengths and
weaknesses. QAOA, for example, enables theoretical solutions to any QUBO
problem with arbitrary precision by increasing the depth of the associated
quantum circuit. Nonetheless, QAs are specifically designed to identify the
lower energy state of an Ising Hamiltonian representing the original QUBO
problem and are better suited for tackling larger problems in terms of the
number of QUBO variables.

In case of constraint optimization problems, the main drawback of the
reformulation as a QUBO consists of associating a penalty term to the
constraints and including them in the objective function. This approach
requires empirical estimation of the penalty terms, which translates to
running the QUBO solver (QA or QAOA) multiple times before achieving a
feasible solution and poses several limitations especially when considering
large problem instances and the imperfections of near-term quantum technology.

This paper presents a novel method for solving the bin packing problem using
quantum computation. Specifically, we introduce an analytical heuristic
approach for estimating penalty multipliers based on the augmented Lagrangian
framework [9] which allows obtaining a complete QUBO formulation without
requiring empirical, instance-based parametrization. To demonstrate the
effectiveness of our proposed method, we conduct experiments using a real
quantum annealer and compare the results with two different state-of-the-art
classical baselines.

## Problem Formulation

The bin packing problem (BPP) is a classic optimization task that involves
packing objects of different sizes into containers, or bins, with a limited
capacity. The goal is to minimize the number of bins needed to pack all the
objects. A mathematical formulation of the bin packing problem can be
expressed as follows: given a set of $n$ items of given integer size (or
weight) $w_{j}\ (j=1,\ldots,n)$ the goal is to pack them into the minimum
number of identical bins of integer capacity $C$.

Let $m$ be any upper bound on the solution value and let introduce
$y_{i},x_{ij}$ two sets of binary variables such that: $y_{i}(i=1,\ldots,m)$
takes the value $1$ if and only if bin $i$ is used in the solution and
$x_{ij}(i=1,\ldots,m;j=1,\ldots,n)$ takes the value $1$ if and only if item
$j$ is packed into bin $i$. A commonly adopted Integer Linear Programming
(ILP) problem formulation is the following[10]:

$\displaystyle\operatorname*{arg\,min}_{x,y}\ $
$\displaystyle\sum_{i=1}^{m}y_{i}$ (1) s.t.
$\displaystyle\sum_{i=1}^{m}x_{ij}=1$ $\displaystyle\forall j=1,\dots,n$ (2)
$\displaystyle\sum_{j=1}^{n}w_{j}x_{ij}\leq Cy_{i}$ $\displaystyle\forall
i=1,\dots,m$ (3) $\displaystyle x_{ij}\in\\{0,1\\}$ $\displaystyle\forall
i=1,\dots,m,\ \forall j=1,\dots,n$ $\displaystyle y_{i}\in\\{0,1\\}$
$\displaystyle\forall i=1,\dots,m$

A practical variant of significant interest is the _online_ bin packing
problem. In this scenario, items of varying sizes are observed sequentially,
and the decision maker must determine whether to select and pack the currently
observed item or let it pass. Each decision is made without the ability to
recall previous decisions. In contrast, the _offline_ bin packing problem
allows for rearranging the items in an attempt to achieve a better packing
arrangement when additional items arrive. However, this approach necessitates
additional storage to hold the items that need to be rearranged.

## Related Works

Classical algorithms for solving the BPP rely on Linear Programming
relaxations and dynamic programming [11, 12]. However, as the number of items
increases, the problem becomes intractable, and even for medium-sized
instances, the optimal solution cannot be computed within a reasonable time
frame. For this reason, several approximation algorithms and heuristics
approaches can be adopted, such as simulated annealing[13, 14], Tabu
search[15], population-based algorithms[16], evolutionary and genetic
heuristics[17, 18, 19, 20, 21] with hyper-heuristics[22, 23, 24, 25], variable
neighborhood search meta-heuristics[26] and ad-hoc crafted heuristics[27, 28,
29].

In addition various alternative reformulations of the BPP have been proposed
to improve the computational performance, such as pseudo-polynomial models
[30, 31, 32]. While these approaches offer more efficient problem formulations
and enable the implementation of solutions that do not scale exponentially
with the input size, they suffer from the drawback of the number of variables
depending on both the number of items and the bin capacity.

More recently, the adoption of quantum computing has been explored for solving
BPP. Existing quantum solutions involve reformulating the original problem as
a QUBO problem and leveraging gate-based quantum computers or quantum
annealing.

At the time of writing this paper, two end-to-end QUBO models have been
proposed for BPP, namely the Pseudo-Polynomial formulation [33] and the
Unbalanced Penalization approach [34]. Alternatively, another existing
approach [35] addresses the BPP through a hybrid approach, using quantum
annealing to solve the sub-problem of filling a single bin, with the chance of
reaching a sub-optimal solution.

##### Pseudo-Polynomial QUBO formulation

The Pseudo-Polynomial QUBO formulation for the BPP [33] is defined by three
sets of binary variables. These variables are employed to represent the
placement of weights in bins, indicate whether bins are empty or not, and
specify the filling levels of the bins. The corresponding Hamiltonian of the
QUBO formulation consists of two weighted components, necessitating empirical
estimation for the weights for each problem instance to avoid unfeasible
solutions. The primary objective of the Hamiltonian is to minimize the number
of used bins, which constitutes the classical objective function. The
constraints Hamiltonian comprises three components. Firstly, it enforces the
condition that each bin must be filled up to a unique level while ensuring
that unused bins remain unfilled. Secondly, it guarantees that every item is
allocated to a bin. Lastly, the third component penalizes configurations in
which bins are overfilled, thereby violating the capacity constraint.
Additionally, an extra term is introduced to account for considerations
regarding only non-empty bins.

A significant drawback of this formulation is its limited scalability
concerning the number of binary variables. Although the introduction of slack
variables enables the formulation to be pseudo-polynomial, it results in the
addition of $nC$ binary variables, rendering the formulation dependent on the
specific problem instance due to the influence of bin capacity on the variable
count. This scalability issue becomes particularly pertinent in the context of
modern Quantum Processing Units (QPUs), which face restrictions in handling
relatively small problem instances owing to qubit topology and connectivity
constraints [36]. Consequently, even for small problem instances, a pseudo-
polynomial Hamiltonian may become intractable when implemented on such quantum
platforms.

Moreover, achieving a well-balanced weight assignment for each term in the
Hamiltonian is of paramount importance to effectively minimize the number of
used bins while satisfactorily adhering to the defined constraints. This
necessitates running the same problem instance multiple times using various
hyperparameter sets in order to optimize the formulation’s performance.

##### Unbalanced penalization formulation

An alternative QUBO model for the BPP[34] introduces an inequality constraint
$g(x)=\sum_{i}l_{i}x_{i}-C\leq 0$ whose violation can be penalized using the
exponential function $e^{g(x)}$. To ensure a valid QUBO model, the exponential
function is expanded up to its second-order Taylor’s term, resulting in the
approximation $e^{g(x)}\approx 1+g(x)+\frac{1}{2}g(x)^{2}$.

Despite its improved efficiency compared to the pseudo-polynomial QUBO
formulation, this model has several limitations. Firstly, it necessitates the
estimation of lambda parameters specific to each problem instance. This
implies running the quantum algorithm multiple times to obtain feasible
solutions for a single problem instance.

Secondly, the model’s performance is evaluated on a limited set of problem
instances, raising concerns about its generalizability to other instances. The
model’s scalability across instances with varying numbers of items is not
demonstrated; the evaluation is restricted to randomly generated instances
with the same number of items.

Furthermore, the experimental testing of the solution relies on QAOA, which
poses restrictions on the number of problem variables due to the challenges
associated with simulating even small-scale quantum systems.

## Contribution

In this work, we presented QAL-BP (_Quantum Augmented Lagrangian method for
Bin Packing_), a novel QUBO formulation for the BPP based on the Augmented
Lagrangian method. QAL-BP is an end-to-end QUBO for the BPP that enables
efficient scaling of logical qubits and the analytical estimation of the
Lagrangian penalty terms. Specifically, we establish a connection between QUBO
models and augmented Lagrangian methods, leveraging advancements in both
fields and fostering potential future synergies.

The proposed formulation offers several advantages. Firstly, it exhibits
independence of the number of variables from the bin capacity. This eliminates
the need for introducing slack variables, which typically increase the number
of logical qubits and render quantum solutions infeasible for execution on
real quantum hardware. Secondly, by employing an augmented Lagrangian
formulation, we can analytically determine the Lagrangian penalty terms for a
specific class of instances without the requirement of running the quantum
solution multiple times, as is often necessary for alternative quantum
approaches. Through experiments conducted on a real quantum annealing device,
we demonstrate the effectiveness of our proposed approach. Thirdly, we proceed
to compare the performance of QAL-BP with state-of-the-art classical
approaches. The results demonstrate that QAL-BP consistently yields feasible
solutions, and in most cases, it leads to the global minimum.

To the best of our knowledge, this marks the first instance of an end-to-end
analytical quantum solution for the BPP that has been rigorously tested across
a diverse set of problem instances, displaying superior performance in
comparison to existing quantum solutions. Furthermore, our results indicate
promising potential concerning state-of-the-art classical solvers,
particularly when more reliable quantum devices become available.

## Methods

### QAL-BP: Quantum Augmented Lagrangian method for Bin Packing

Augmented Lagrangian methods are a class of algorithms used to solve
constrained optimization problems by incorporating the constraints into the
objective function through penalty terms. Consider a constrained minimization
problem of the form:

$\displaystyle\begin{split}\text{min}\ &f(\textbf{x})\\\ \text{s.t.
}&c_{i}(\textbf{x})=\textbf{b}\ \forall i\in\mathcal{D}\end{split}$ (4)

where $x$ is a candidate solution, $c_{i}(\textbf{x})=\textbf{b}$ are a set of
equality constraints and $\mathcal{D}$ is the set of constraints, the
augmented Lagrangian method consists of defining an unconstrained problem of
the form:

$\text{min}\
\Phi(\textbf{x})=f(\textbf{x})+\sum_{i\in\mathcal{D}}\rho_{i}(c_{i}(\textbf{x})-\textbf{b})^{2}+\sum_{i\in\mathcal{D}}\lambda_{i}(c_{i}(\textbf{x})-\textbf{b})$
(5)

where $\rho_{i}$, $\lambda_{i}$ for {i = 1, …, $|\mathcal{D}$|} are the
lagrangian multipliers.

In practice, when using the augmented Lagrangian approach, it is common to
introduce additional constraints that do not alter the set of feasible
solutions but aid in faster convergence of the solver (REF).

For the BPP, we introduce the following additional constraints:

$\sum_{j=1}^{n}y_{i}x_{ij}=\sum_{j=1}^{n}x_{ij}\quad\forall i=1,\dots,m$ (6)

This constraint implies that if a bin $i$ is not used ($y_{i}=0$), it cannot
contain any items ($\sum_{j=1}^{n}x_{ij}=0$).

Therefore, the Quantum Augmented Lagrangian method for Bin Packing (QAL-BP)
embeds the constraints using the augmented Lagrangian method as follows:

$\displaystyle\text{argmin}_{x,y}\ $ $\displaystyle\delta\sum_{i=1}^{m}y_{i}$
(7)
$\displaystyle+\sum_{i=1}^{m}\lambda_{i}\left(\sum_{j=1}^{n}w_{j}x_{ij}-c_{i}y_{i}\right)$
(8)
$\displaystyle+\sum_{i=1}^{m}\rho_{i}\left(\sum_{j=1}^{n}w_{j}x_{ij}-c_{i}y_{i}\right)^{2}$
(9)
$\displaystyle+\theta\sum_{j=1}^{n}\left(\sum_{i=1}^{m}x_{ij}-1\right)^{2}$
(10)
$\displaystyle+\gamma\sum_{i=1}^{m}\left(1-y_{i}\right)\sum_{j=1}^{n}x_{ij}$
(11) $\displaystyle\text{s.t. }x_{ij}$ $\displaystyle\in\\{0,1\\}\qquad\forall
i=1\dots,m,\forall j=1\dots,n$ $\displaystyle y_{i}$
$\displaystyle\in\\{0,1\\}\qquad\forall i=1\dots,m$ $\displaystyle\delta,$
$\displaystyle\lambda,\rho,\theta,\gamma\geq 0$

The penalties (8) and (9) represent the augmented Lagrangian expansion of (3).
These terms impose a penalty of $\lambda_{i}s_{i}+\rho_{i}s_{i}^{2}$ for
infeasible configurations while providing a negative reward to the solver for
feasible configurations. Accurate estimation of $\lambda_{i}$ and $\rho_{i}$
values is crucial for correctly modeling the solution space.

Similarly, in Eq. (10), $\theta$ represents a penalty for not placing an item
$j$, and penalizes $(k-1)\theta$ when item $j$ is placed $k$ times. Notably,
this penalty term is not an augmented Lagrangian expansion of (2), but rather
a pure squared penalty. This is because we do not wish to reward the solver
when an item is not placed at all.

Finally, the term (11) represents the penalty associated with the redundant
constraints (6) which imposes a penalty of $\sum_{j\in\mathcal{J}}x_{ij}$ when
the set of items $\mathcal{J}$ is assigned to bin $i$ without setting the
corresponding $y_{i}$ to 1.

It is important to note that the standard augmented Lagrangian approach
typically transforms inequality constraints into equality constraints using
slack variables, which are then incorporated into the Lagrangian as shown in
(5). However, the proposed QUBO formulation in this study does not involve
slack variables but directly utilizes the capacity constants $c_{i}$. This
aspect provides a significant advantage over the pseudo-polynomial
approach[33].

### Penalties estimation

When incorporating constraints into the objective function, the estimation of
penalty multipliers is typically carried out by testing a large set of
parameters which requires running the algorithm multiple times with different
parameter values to find a feasible solution for the specific problem instance
and increase dramatically the time cost to find good solutions.

Given the QAL-BP formulation and its corresponding set of constraints, we
propose an analytical estimation of the penalty multipliers. The conditions
are designed based on approximate worst-case reasoning, aiming to yield
optimal or slightly sub-optimal solutions for most instances of the BPP. The
following is a set of heuristic conditions that consider each penalty
multiplier individually, with the exception of the pair $\lambda,\rho$.
Considering the $i$-th bin, the correspondent augmented Lagrangian term is
given by $\lambda_{i}(s_{i}-c_{i}y_{i})+\rho_{i}(s_{i}-c_{i}y_{i})^{2}$. When
$y_{i}=0$, i.e, …using the smallest bin usage amount should be at least as
expensive as using the bin:

$\lambda_{i}(w_{min}-0)+\rho_{i}(w_{min}-0)^{2}\geq 1,$ (12)

where $w_{min}$ is the smallest item weight, i.e. $w_{min}=\min\\{w_{j}\\}$.
If this condition is satisfied, using more capacity makes the solver choose to
set $y_{i}=1$. On the contrary, if $y_{i}=0$ and no capacity then the
condition in eq. (12) is true. Let’s consider now the case where $y_{i}=1$,
meaning that the bin $y_{i}$ is included in the candidate solution. In this
case, exceeding the capacity by any amount should be at least as expensive as
using one more bin, i.e.:

$\lambda_{i}(c_{i}+w_{min}-c_{i})+\rho_{i}(c_{i}+w_{min}-c_{i})^{2}\geq 1$
(13)

which conveniently is the same condition as Eq. (12). Lastly, it is necessary
to identify a solution space that contains only feasible solutions. In this
case, the Lagrangian term needs to provide a positive reward (i.e., negative
cost) if the constraint is satisfied. It is also necessary that such a reward
be small enough that it does not provide an incentive for using another bin,
i.e.:

$\lambda_{i}(-\frac{c_{i}}{2})+\rho_{i}(-\frac{c_{i}}{2})^{2}\geq 0$ (14)

Basically, here we are fitting a quadratic function using the conditions (12)
and (14) to approximate the solution space excluding unfeasible solutions. We
can therefore obtain values for $\lambda_{i}$ and $\rho_{i}$ by stating all
conditions for their least restrictive values:

$\displaystyle w_{min}\lambda_{i}+w_{min}^{2}\rho_{i}$ $\displaystyle=1;$ (15)
$\displaystyle-\frac{c_{i}}{2}\lambda_{i}+\frac{c_{i}}{4}^{2}\rho_{i}$
$\displaystyle=0,$ (16)

which leads to an analytical formulation of the form:

$\displaystyle\lambda_{i}$
$\displaystyle=\frac{c_{i}}{w_{min}\left(2w_{min}+c_{i}\right)}$ (17)
$\displaystyle\rho_{i}$
$\displaystyle=\lambda_{i}\frac{2}{c_{i}}=\frac{2}{w_{min}\left(2w_{min}+c_{i}\right)}$
(18)

The next step is to calibrate $\theta_{j}$. The abstract Lagrangian term
associated with item $j$ is $\theta_{j}(p_{j}-1)^{2}$, where
$p_{j}\in\mathbb{N}$ is the number of times item $j$ has been assigned to a
bin. At this stage, only the assignment is considered and not capacities, it
is possible to ignore the index $j$ when define the $\theta$ parameter.
Moreover, we want to force the solver to assign all items to at maximum one
bin, so the penalty should increase when $p_{j}\neq 1$. Thus, in case
$p_{j}=1$, the following condition holds:

$\displaystyle\theta(1-1)^{2}\geq-1.$ (19)

In case $p_{j}=0$, i.e., item $j$ is not assigned to any bin, then the
associated penalty needs to be greater than the cost of opening a new bin,
i.e.,

$\displaystyle\theta\geq 2.$ (20)

Also, the parameter $\gamma$ needs to be set. The abstract Lagrangian term
associated with the $\gamma$ term is $\gamma(1-y_{i})k_{i}$, which comes into
play only when $y_{i}=0$ and $k_{i}\neq 0$ by adding $k_{i}$ times the penalty
$\gamma$. We want the minimum penalty to be at least equal to the cost of
opening a new bin, i.e.,

$\displaystyle\gamma\geq 1.$ (21)

The final parameter to be estimated is $\delta$. Although it is not strictly a
penalty term since it serves as a multiplier of the objective function,
including it is beneficial for controlling the other parameters. One purpose
of $\delta$ is to prevent issues that may arise when working with extremely
small numbers.

Another reason for incorporating this multiplier is to address the undesirable
behavior of the model in certain item configurations. In some cases, the model
may favor configurations where one or more bins are slightly overfilled due to
the high cost associated with opening a new bin. Naturally, this behavior is
contingent on the specific combination of instance weights relative to the bin
capacity, as well as their number.

To rectify this behavior, the following requirements must be met:

$\displaystyle\delta\leq\lambda s_{min}+\rho s_{min}^{2},$ (22)

where $s_{min}$ is the minimum capacity that can be exceeded. Thus, the cost
of opening a new bin must be less than the cost of overfilling an already open
bin, of the smallest possible amount, i.e., $s_{min}\geq 1$.

### Model analysis

When solving QUBO problems with quantum computing, the number of binary
variables of the problem formulation corresponds to the number of logical
qubits to use in the quantum computer. Therefore, having an efficient
formulation that minimizes the number of variables without limiting the range
of possible solutions is crucial for the adoption of quantum approaches. In
terms of variable count, the QAL-BP model is more efficient than the pseudo-
polynomial one [33], and equivalent, with respect to the _Unbalanced
penalization formulation_[34]. Specifically, for a given problem instance
$BPP(n,C)$, the total number of variables is equal to the number of bins $m$
plus $n\times m$ decision variables representing the assignment of a specific
item to a specific bin. Thus, in the worst case, $m=n$, resulting in $n(n+1)$
binary variables.

Furthermore, from a methodological point of view, the QAL-BP approach offers a
twofold advantage. Firstly, the number of binary variables is not affected by
bin capacities and item weights, as observed in the pseudo-polynomial
formulation[37]. Secondly, the reduced number of variables enables the
execution of the QUBO problem using a smaller set of logical qubits, rendering
it suitable for current QPUs.

Fig. 1 shows a comparison of the performance, in terms of the number of
variables, between the QAL-BP and the pseudo-polynomial formulation.

Figure 1: Comparative analysis of variable growth in the Pseudo-Polynomial and
Augmented Lagrangian models concerning the number of items and bin capacity.
Three distinct values of bin capacity (C) are explored. The continuous dark
red line represents the upper limit for QUBO problems represented by fully
connected graph that can be mapped in the D-Wave Advantage Quantum Processing
Unit (QPU) equipped with 5640 qubits.

## Evaluation

#### Experimental Settings

##### Data

The experiments are performed on a set of eight classes of randomly generated
instances, ranging from $3$ to $10$ items, with corresponding weights ranging
from $4$ to $10$ and fixed bin capacity equal to $10$. This experimental
choice allows finding the best compromise between exploring problem instances
of different sizes while taking into account the limitations of modern QPUs
which are restricted to non-sparse QUBO problems with up to $180$ binary
variables[38]. In particular, five different problem instances are generated
for given a fixed number of items, then varying the weights, which are solved
using different approaches. The generated instances are shown in Table 1.

Seed |  N. items |  Weights |  Lower bound |  Instance name
---|---|---|---|---
23 |  3 |  [4, 8, 6] |  2 |  (3, 23)
4 |  [8, 5, 4, 8] |  3 |  (4, 23)
5 |  [4, 4, 8, 8, 9] |  3 |  (5, 23)
6 |  [7, 5, 5, 5, 4, 9] |  4 |  (6, 23)
7 |  [9, 7, 8, 6, 9, 6, 7] |  5 |  (7, 23)
8 |  [4, 5, 7, 5, 6, 4, 6, 4] |  4 |  (8, 23)
9 |  [7, 6, 8, 4, 8, 4, 9, 6, 4] |  6 |  (9, 23)
10 |  [5, 8, 6, 7, 10, 9, 4, 10, 7, 4] |  7 |  (10, 23)
42 |  3 |  [4, 8, 6] |  2 |  (3, 42)
4 |  [7, 7, 10, 4] |  3 |  (4, 42)
5 |  [8, 5, 4, 7, 10] |  4 |  (5, 42)
6 |  [9, 9, 9, 9, 7, 4] |  5 |  (6, 42)
7 |  [9, 7, 7, 6, 5, 10, 9] |  5 |  (7,42)
8 |  [8, 6, 9, 7, 7, 7, 5, 4] |  5 |  (8, 42)
9 |  [7, 10, 4, 10, 9, 5, 8, 5, 9] |  7 |  (9, 42)
10 |  [8, 6, 4, 10, 7, 10, 8, 9, 9, 5] |  7 |  (10, 42)
123 |  3 |  [4, 8, 8] |  2 |  (3, 123)
4 |  [4, 10, 5, 5] |  3 |  (4, 123)
5 |  [5, 6, 5, 6, 9] |  3 |  (5, 123)
6 |  [7, 10, 7, 5, 9, 9] |  5 |  (6, 123)
7 |  [10, 10, 4, 7, 5, 5, 5] |  5 |  (7, 123)
8 |  [9, 9, 5, 6, 9, 5, 8, 7] |  6 |  (8, 123)
9 |  [10, 9, 5, 9, 9, 5, 7, 9, 5] |  7 |  (9, 123)
10 |  [5, 5, 4, 7, 4, 8, 6, 5, 6, 4] |  5 |  (10, 123)
90 |  3 |  [8, 6, 4] |  2 |  (3, 90)
4 |  [8, 5, 7, 6] |  3 |  (4, 90)
5 |  [6, 7, 8, 7, 4] |  3 |  (5, 90)
6 |  [7, 8, 9, 9, 10, 6] |  5 |  (6, 90)
7 |  [6, 4, 4, 4, 8, 9, 6] |  4 |  (7, 90)
8 |  [7, 10, 8, 8, 8, 5, 5, 8] |  6 |  (8, 90)
9 |  [9, 6, 4, 10, 10, 5, 4, 4, 6] |  6 |  (9, 90)
10 |  [9, 6, 8, 7, 8, 10, 9, 6, 9, 10] |  8 |  (10, 90)
510 |  3 |  [5, 8, 6] |  2 |  (3, 510)
4 |  [7, 9, 5, 5] |  3 |  (4, 510)
5 |  [6, 10, 4, 9, 4] |  3 |  (5, 510)
6 |  [5, 5, 9, 10, 8, 6] |  4 |  (6, 510)
7 |  [9, 7, 9, 4, 10, 10, 8] |  6 |  (7, 510)
8 |  [9, 10, 8, 9, 4, 4, 9, 5] |  6 |  (8, 510)
9 |  [5, 9, 10, 9, 7, 8, 4, 10, 6] |  7 |  (9, 510)
10 |  [10, 5, 9, 5, 8, 9, 7, 4, 6, 9] |  7 |  (10, 510)

Table 1: The selected set of bin packing instances serves as a testbed to
evaluate all models under consideration. Each instance in the set is
characterized by specific parameters, which are organized into columns for
easy reference and analysis. The first column denotes the seed of the random
generator utilized in generating the instances. The second column indicates
the number of items that need to be placed within the bins. The third column
comprises an array representing the weight of each individual item. The fourth
column provides the $L_{1}$ lower bound for the given instance[39]. Lastly,
the last column assigns a unique label to facilitate identification and
reference to each specific instance.

##### Solving methods

The quantum solver employed in this study is the _D-Wave Advantage 4.1_ ,
featuring a total of $5640$ physical qubits. For the purpose of comparison,
two classical solvers are also utilized. The first is simulated annealing,
which is considered the classical counterpart of quantum annealing from a
methodological perspective. Both quantum annealing and simulated annealing are
available in the D-Wave Python library, enabling a thorough assessment of the
correctness of the QAL-BP formulation without encountering any errors inherent
to real quantum hardware. In addition, we also solve the ILP formulation from
Eq. (1)-(3) via the Gurobi optimizer [40], as a representative for state-of-
the-art solution that relies on the branch-and-bound technique to efficiently
find the optimal solution.

As previously mentioned, the rationale behind choosing quantum annealing over
alternative quantum approaches, such as QAOA, is to enable a direct
performance comparison between the most potent and dependable current quantum
technology and the state-of-the-art classical optimizer. This comparison aims
to demonstrate the comparable capability of current quantum computation in
solving the BPP in relation to the best available classical solution. It is
important to recognize that quantum hardware is still in its nascent stage,
making such a comparative analysis critical in assessing the advancements and
potential of quantum computing in tackling combinatorial optimization problems
like BPP.

##### Metrics

In order to assess the performance of QAL-BP, we employ three distinct
metrics: _Time-To-Solution_ (TTS), _solution quality_ , and _probability_. TTS
represents the time required by the solver to produce the final solution,
measured in microseconds ($\mu$s). For quantum annealing (QA), the
_qpu_sampling_time_ 111docs.dwavesys.com/docs/latest/c_qpu_timing.html is
considered, for an independent assessment of possible connection delays with
the cloud and possible recalibrations. Differently from QA, simulated
annealing (SA) and Gurobi run locally. Thus the TTS of SA is the time to run
the annealing function, while for Gurobi, TTS is calculated as the time
required to obtain the solution for a given problem instance.

The _solution quality_ is defined as the number of bins utilized in a given
solution. The lower the number, the better the solution. The _probability_
measures the odds, for the solution with minimum energy, to be a feasible one.
In particular, this probability is computed, for each class of instances, as
the ratio between the number of instances for which the minimum represents a
feasible solution and the total number of instances in that class. This metric
is only meaningful for simulated and quantum annealing and allows assessing
the performance of current quantum technology in returning the solution for
the QAL-BP formulation.

##### Models parameters

According to the analytical penalty estimation previously described, the
multipliers are determined as follows: $\delta=0.15;\ \lambda=0.1389;\
\rho=0.0278;\ \theta=2;\ \gamma=1$.

### Results

Figure 2 illustrates the experimental results concerning TTS. The performance
of Simulated Annealing (SA) deteriorates rapidly with the number of items.

Figure 2: TTS comparison between Quantum Annealing (QA), Simulated Annealing
(SA) and Gurobi.

Comparing QA with Gurobi, the latter consistently outperforms the former.
However, it is important to note that the asymptotic behavior of these two
methods appears to be equivalent. Due to the limitations in running larger
problem instances, a conclusive assessment of performance cannot be made.
Nevertheless, the quantum annealing’s ability to scale with a runtime linear
in the input size has been observed in other contexts, such as coalition
formation in multi-agent systems [4, 3]. On the other hand, Gurobi, relying on
an exact branch-and-bound approach, is guaranteed to have exponential worst-
case complexity, since the bin packing problem is NP-hard.

Considering these aspects, the results from QA suggest the potential for
outperforming the state-of-the-art classical solver Gurobi once more reliable
quantum technology becomes available. This highlights the promise of quantum
annealing in achieving superior performance in solving combinatorial
optimization problems compared to classical approaches when leveraging
advancements in quantum hardware.

In terms of solution quality, QA demonstrates proficiency in finding the
correct solution only for small problem instances, while encountering
challenges as the problem size increases. Figures 3 and 4 respectively present
the number of used bins in a given solution (referred to as solution quality)
and the probability of that solution being feasible. To gain a comprehensive
understanding of QA’s performance, these two plots must be jointly analyzed.

Figure 3: Comparison of the number of bins of the best solution found by
Gurobi, Simulated and Quantum annealing solvers.

In Fig. 3, Gurobi solutions serve as a benchmark, representing the correct
minimum number of bins. The results obtained by SA confirm the correctness of
the QAL-BP formulation in identifying only feasible solutions for any problem
instance and, in most cases, the best solution of QAL-BP corresponds to the
global optimum of the optimization problem. The only exceptions are observed
in instances $(8,23)$, $(10,123)$, and $(9,510)$, where the solutions returned
by SA, although feasible, are sub-optimal.

With regard to the results obtained by QA, it is observed that for relatively
larger instances (instances with more than $7$ items), QA tends to select a
higher number of bins compared to Gurobi and SA. Specifically, in the
instances $(10,23)$, $(8,123)$, $(10,23)$, and $(9,90)$, QA utilizes a higher
number of bins compared to SA and Gurobi, leading to sub-optimal solutions.
Additionally, in $(7,23)$, $(9,23)$, and $(10,42)$, the solutions returned by
QA are not feasible, as the number of bins is lower than the number of bins
selected by Gurobi. This implies that in these cases, the bins are overfilled
beyond their capacity, resulting in infeasible solutions. Again, this
limitation is due attributed to the current constraints of quantum hardware
when dealing with large problem instances. A similar situation occurs in the
case of $(9,510)$, where QA achieves a different solution compared to SA, but
this solution is inconsistent with the QAL-BP formulation.

To further validate the impediments of current QA hardware in correctly
finding the solutions for the QAL-BP formulation, Fig. 4 presents the
probability of obtaining a feasible solution as the average of different
problem instances of the same size. For instances with more than $7$ items,
the probability of obtaining a feasible solution significantly decreases for
QA, while it remains consistently at $1$ for SA.

In summary, the experimental results demonstrate that the QAL-BP formulation
allows for the accurate analytical estimation of penalty parameters, leading
to almost all cases resulting in the global optimum. However, as the problem
size increases, QA’s performance diminishes due to limitations in current
quantum hardware. To fully leverage the potential of QA in effectively
tackling larger BPP instances, further advancements in quantum technology and
solver optimization are imperative.

Figure 4: Probability of being a feasible solution for the minimum obtained
using annealing-based approaches (simulated and quantum).

## Conclusion

This paper introduces QAL-BP, a novel quantum formulation based on the
augmented Lagrangian method for efficiently solving the Bin Packing Problem
(BPP) using quantum annealing. QAL-BP offers an analytical estimation of
penalty terms in the model for a specific class of problem instances,
eliminating the need for recursive approximation methods to empirically
estimate Lagrangian multipliers. This enhancement amplifies the
generalizability of our approach to diverse input instances and improves
efficiency by reducing the number of QUBO variables compared to alternative
quantum formulations.

We demonstrated the effectiveness of our approach by solving larger problem
instances than any previous QUBO formulation for the BPP. Additionally, we
present the first experimental comparison of classical and quantum solutions
for the BPP, validating that QAL-BP is an analytically correct QUBO
formulation obviating the need for empirical estimation of penalty terms.
Nevertheless, while our implementation on a quantum annealer does not
outperform the state-of-the-art classical solver Gurobi, its TTS exhibits
efficient scaling as the problem size increases, considering the current
limitations of available quantum technology.

However, several limitations and challenges remain. Firstly, the
generalizability of our model to generic BPP instances or other combinatorial
optimization problems requires further investigation. Secondly, the limited
number of qubits on current quantum annealers poses a significant challenge,
restricting the size of problem instances that can be effectively solved.
Consequently, testing our model on larger instances and evaluating scalability
across a wider range of inputs is currently unattainable. Furthermore, noise
and errors in the quantum annealer significantly impact the quality of
provided solutions, particularly evident when dealing with larger problems, as
demonstrated in experimental results compared to simulated annealing.

To address these challenges, future research will explore advanced quantum
hardware with improved qubit accuracy and a greater qubit count. Another
promising avenue involves investigating hybrid quantum annealing approaches
that leverage classical and quantum methods in tandem, facilitating the
solution of larger problem sizes beyond the capabilities of current QPUs.
These endeavors are critical to further harnessing the potential of quantum
computing in combinatorial optimization problems and propelling the field
forward.

## References

  * [1] Garey, M. R. & Johnson, D. S.  “ strong ” np-completeness results: Motivation, examples, and implications.  _J. ACM_ 25, 499–508, DOI: 10.1145/322077.322090 (1978).
  * [2] Delorme, M., Iori, M. & Martello, S.  Bin packing and cutting stock problems: Mathematical models and exact algorithms.  _European Journal of Operational Research_ 255, 1–20, DOI: https://doi.org/10.1016/j.ejor.2016.04.030 (2016).
  * [3] Venkatesh, S. M., Macaluso, A. & Klusch, M.  Gcs-q: Quantum graph coalition structure generation.  In _International Conference on Computational Science_ , 138–152 (Springer, 2023).
  * [4] Venkatesh, S. M., Macaluso, A. & Klusch, M.  Bilp-q: quantum coalition structure generation.  In _Proceedings of the 19th ACM International Conference on Computing Frontiers_ , 189–192 (2022).
  * [5] Venkatesh, S. M., Macaluso, A. & Klusch, M.  Quacs: Variational quantum algorithm for coalition structure generation in induced subgraph games.  _arXiv preprint arXiv:2304.07218_ (2023).
  * [6] Macaluso, A., Klusch, M., Lodi, S. & Sartori, C.  Maqa: a quantum framework for supervised learning.  _Quantum Information Processing_ 22, 159 (2023).
  * [7] Macaluso, A., Orazi, F., Klusch, M., Lodi, S. & Sartori, C.  A variational algorithm for quantum single layer perceptron.  In _International Conference on Machine Learning, Optimization, and Data Science_ , 341–356 (Springer, 2022).
  * [8] Farhi, E., Goldstone, J. & Gutmann, S.  A quantum approximate optimization algorithm.  _arXiv preprint arXiv:1411.4028_ (2014).
  * [9] Hestenes, M. R.  Multiplier and gradient methods.  _Journal of Optimization Theory and Applications_ 4, 303–320, DOI: 10.1007/BF00927673 (1969).
  * [10] Martello, S. & Toth, P.  _Knapsack Problems: Algorithms and Computer Implementations_ (John Wiley and sons, 1990).
  * [11] Eisemann, K.  The trim problem.  _Management Science_ 3, 279–284 (1957).
  * [12] Gilmore, R. & Gomory, R.  A linear programming approach to the cutting stock problem i.  _Oper Res_ 9, DOI: 10.1287/opre.9.6.849 (1961).
  * [13] Kampke, T.  Simulated annealing: Use of new tool in bin packing.  _Ann. Oper. Res._ 16, 327–332, DOI: 10.1007/BF02283751 (1988).
  * [14] Loh, K.-H., Golden, B. & Wasil, E.  Solving the one-dimensional bin packing problem with a weight annealing heuristic.  _Computers & Operations Research_ 35, 2283–2291, DOI: https://doi.org/10.1016/j.cor.2006.10.021 (2008).  Part Special Issue: Includes selected papers presented at the ECCO’04 European Conference on combinatorial Optimization.
  * [15] Scholl, A., Klein, R. & Jürgens, C.  Bison: a fast hybrid procedure for exactly solving the one-dimensional bin packing problem.  _Computers and Operations Research_ 24, 627–645 (1997).
  * [16] Vahrenkamp, R.  Random search in the one-dimensional cutting stock problem.  _European Journal of Operational Research_ 95, 191–200, DOI: https://doi.org/10.1016/0377-2217(95)00198-0 (1996).
  * [17] Burke, E., Hyde, M. & Kendall, G.  Evolving bin packing heuristics with genetic programming.  860–869, DOI: 10.1007/11844297_87 (2006).
  * [18] Falkenauer, E.  A hybrid grouping genetic algorithm for bin packing.  _Journal of Heuristics_ 2, DOI: 10.1007/BF00226291 (1996).
  * [19] Falkenauer, E. & Delchambre, A.  A genetic algorithm for bin packing and line balancing.  In _Proceedings 1992 IEEE International Conference on Robotics and Automation_ , 1186–1192 vol.2, DOI: 10.1109/ROBOT.1992.220088 (1992).
  * [20] Quiroz, M. _et al._ A grouping genetic algorithm with controlled gene transmission for the bin packing problem.  _Computers & Operations Research_ 55, 52–64, DOI: 10.1016/j.cor.2014.10.010 (2014).
  * [21] Sim, K. & Hart, E.  Generating single and multiple cooperative heuristics for the one dimensional bin packing problem using a single node genetic programming island model.  DOI: 10.1145/2463372.2463555 (2013).
  * [22] Bai, R., Blazewicz, J., Burke, E., Kendall, G. & Mccollum, B.  B.: A simulated annealing hyper-heuristic methodology for flexible decision support.  _4OR_ 10, DOI: 10.1007/s10288-011-0182-8 (2012).
  * [23] Lopez-Camacho, E., Terashima-Marin, H. & Ross, P.  A hyper-heuristic for solving one and two-dimensional bin packing problems.  _Proceedings of the 13th annual conference companion on Genetic and evolutionary computation_ (2011).
  * [24] Sosa, A., Terashima-Marín, H., Ortiz-Bayliss, J. C. & Conant-Pablos, S.  Grammar-based selection hyper-heuristics for solving irregular bin packing problems:.  111–112, DOI: 10.1145/2908961.2908970 (2016).
  * [25] Sim, K., Hart, E. & Paechter, B.  A hyper-heuristic classifier for one dimensional bin packing problems: Improving classification accuracy by attribute evolution.  In Coello, C. A. C. _et al._ (eds.) _Parallel Problem Solving from Nature - PPSN XII_ , 348–357 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).
  * [26] Gomez-Meneses, P. & Randall, M.  A hybrid extremal optimisation approach for the bin packing problem.  vol. 5865, DOI: 10.1007/978-3-642-10427-5_24 (2009).
  * [27] Eilon, S. & Christofides, N.  The loading problem.  _Management Science_ 17, 259–268, DOI: 10.1287/mnsc.17.5.259 (1971).  https://doi.org/10.1287/mnsc.17.5.259.
  * [28] Gupta, J. N. D. & Ho, J. C.  A new heuristic algorithm for the one-dimensional bin-packing problem.  _Production Planning & Control_ 10, 598–603, DOI: 10.1080/095372899232894 (1999).  https://doi.org/10.1080/095372899232894.
  * [29] Lewis, R.  A general-purpose hill-climbing method for order independent minimum grouping problems: A case study in graph colouring and bin packing.  _Computers & Operations Research_ 36, 2295–2310, DOI: 10.1016/j.cor.2008.09.004 (2009).
  * [30] M.R.Rao.  On the cutting stock problem.  _Journal of the Computer Society of India_ (1976).
  * [31] Dyckhoff, H.  A new linear programming approach to the cutting stock problem.  _Operations Research_ 29, 1092–1104, DOI: 10.1287/opre.29.6.1092 (1981).  https://doi.org/10.1287/opre.29.6.1092.
  * [32] Stadtler, H.  A comparison of two optimization procedures for 1- and 1 1/2-dimensional cutting stock problems.  _OR Spektrum_ 10, 97–111, DOI: 10.1007/BF01720208 (1988).
  * [33] Lodewijks, B.  Mapping np-hard and np-complete optimisation problems to quadratic unconstrained binary optimisation problems, DOI: 10.48550/ARXIV.1911.08043 (2019).
  * [34] Montanez-Barrera, A., Willsch, D., Maldonado-Romo, A. & Michielsen, K.  Unbalanced penalization: A new approach to encode inequality constraints of combinatorial problems for quantum optimization algorithms, DOI: 10.48550/ARXIV.2211.13914 (2022).
  * [35] Garcia-de Andoin, M., Oregi, I., Villar-Rodriguez, E., Osaba, E. & Sanz, M.  Comparative benchmark of a quantum algorithm for the bin packing problem, DOI: 10.48550/ARXIV.2207.07460 (2022).
  * [36] Kuramata, M., Katsuki, R. & Nakata, K.  Larger sparse quadratic assignment problem optimization using quantum annealing and a bit-flip heuristic algorithm.  In _2021 IEEE 8th International Conference on Industrial Engineering and Applications (ICIEA)_ , DOI: 10.1109/iciea52957.2021.9436749 (IEEE, 2021).
  * [37] Lucas, A.  Ising formulations of many NP problems.  _Frontiers in Physics_ 2, DOI: 10.3389/fphy.2014.00005 (2014).
  * [38] Kuramata, M., Katsuki, R. & Nakata, K.  Larger sparse quadratic assignment problem optimization using quantum annealing and a bit-flip heuristic algorithm.  556–565, DOI: 10.1109/ICIEA52957.2021.9436749 (2021).
  * [39] Martello, S. & Toth, P.  Lower bounds and reduction procedures for the bin packing problem.  _Discrete Applied Mathematics_ 28, 59–70, DOI: https://doi.org/10.1016/0166-218X(90)90094-S (1990).
  * [40] Gurobi optimization.  Website.  Accessed: 22/07/2023.

## Code availability

The code is available on GitHub at the link: https://github.com/Lorenz92/QAL-
BP

## Author contributions statement

L.C. developed the theoretical formalism, performed the analytic calculations
and numerical experiments. Both A.M. and M.L. supervised the project. L.C. and
A.M. wrote the manuscript in consultation with M.L. All authors provided
critical feedback and helped shape the research, analysis, and manuscript.

## Competing interests

The authors declare no competing interests.


longtable Ginwidth=,height=,keepaspectratio

# Revealing Patterns of Symptomatology in Parkinson’s Disease: A Latent Space
Analysis with 3D Convolutional Autoencoders

E. Delgado de las Heras

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada

F.J. Martinez-Murcia

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada
<EMAIL_ADDRESS>

I.A. Illán

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada

C. Jiménez-Mesa

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada

D. Castillo-Barnes

Dpt. of Communications Engineering, University of Malaga.
Malaga

J. Ramírez

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada

J.M. Górriz

Dpt. Signal Theory, Networking and Communications, University of Granada.
Granada
0000-0001-8146-7056

###### Abstract

This work proposes the use of 3D convolutional variational autoencoders
(CVAEs) to trace the changes and symptomatology produced by neurodegeneration
in Parkinson’s disease (PD). In this work, we present a novel approach to
detect and quantify changes in dopamine transporter (DaT) concentration and
its spatial patterns using 3D CVAEs on Ioflupane (FPCIT) imaging. Our approach
leverages the power of deep learning to learn a low-dimensional representation
of the brain imaging data, which then is linked to different symptom
categories using regression algorithms. We demonstrate the effectiveness of
our approach on a dataset of PD patients and healthy controls, and show that
general symptomatology (UPDRS) is linked to a d-dimensional decomposition via
the CVAE with R2>0.25. Our work shows the potential of representation learning
not only in early diagnosis but in understanding neurodegeneration processes
and symptomatology.

##  1 Introduction

Parkinson’s disease (PD) is a progressive neurodegenerative disorder that
affects more than 6.2 million people worldwide [1]. PD is characterized by a
loss of dopamine-producing neurons in the brain, causing tremors, stiffness,
and cognitive decline among other symptoms. FPCIT (ioflupane) SPECT is the
most widely extended neuroimaging technique for the diagnosis of PD. FPCIT
binds to the presynaptic dopamine transporters (DaTs), allowing to visualize
and quantify the DaT concentration at the striata, that characterizes PD [2].
This modality may allow to tackle one of the biggest challenges for PD: the
early detection and monitoring of disease progression [3].

Recent advances in deep learning, and in particular, convolutional neural
networks (CNNs) have shown promising results in medical imaging applications
such as segmentation, registration, and classification [4]. Autoencoders
(AEs), a type of self-supervised neural network, can learn a low-dimensional
representation of high-dimensional data, making them well-suited for image
compression, denoising, and anomaly detection [5].

In this paper, we propose a novel approach to detect and quantify subtle
changes in DaT concentration and distribution in the brain using 3D
convolutional variational AEs (CVAEs). Our approach leverages the power of
deep learning to learn a low-dimensional representation of brain imaging data,
which enables us to longitudinally compare images and identify patterns of
change that are indicative of neurodegeneration. We demonstrate the
effectiveness of our approach on a dataset of PD patients and healthy
controls, and show that latent spaces can capture the variability of
individual symptom categories, as well as the overall disease stage.

##  2 Materials and Methods

###  2.1 Dataset

Data for this study were obtained from the Parkinson’s Progression Markers
Initiative (PPMI) database. For the most recent information on the study,
please visit www.ppmi-info.org. We used the standardized “Original Cohort BL
to Year 5” that includes individuals initially diagnosed as either controls
(CTL) or PD affected subjects with varying levels of severity. After removing
subjects with no FPCIT scans, the CTL group comprises 101 males and 53
females, while the PD group comprises 284 males and 159 females, followed for
up to 5 years, with a total 1399 sessions of available data for studying the
progression of imaging biomarkers and their link to PD-specific progression
indicators. Symptomatology is assessed via the MDS-UPDRS scale [6], either in
its aggregated form (UPDRS -total-) or its 4 parts: 1) non-motor aspects of
daily living; 2) motor aspects of daily living (tremor, walking, etc); 3)
motor examination (rigidity, posture, gait etc) and 4) motor complications
(dyskinesia, fluctuation items).

FPCIT scans labeled as “Reconstructed” were used, having all of them a
consistent orientation and similar sizes (for a scan of shape $91\times
109\times 91$). However, no spatial normalization is performed, as it was
demonstrated to have small impact when using CNN architectures [4]. Intensity
was normalized by substracting the background average intensity, and then
dividing all voxel intensities by a non-specific reference, defined as the
average voxel intensity of the cerebellum and occipital lobe. Finally, to
favor convergence of the CVAE training, the upper-end values of the intensity
distribution were compressed using a sigmoid function.

###  2.2 3D Convolutional Variational Autoencoders

An autoencoder is a neural network architecture that consists of an encoder
network and a decoder network. In this work we used a 3D convolutional
encoder-decoder architecture, capable of dealing with volumetric 3D images.
The 3D encoder can capture spatial features at multiple scales, similarly to
well-known 3D convolutional NNs [4]. Its output is a lower-dimensional
representation which encodes the most salient features of the input data.
Under the manifold hypothesis [5], this representation is often considered as
coordinates in a “latent space” underlying the dataset. The decoder network
performs the inverse operation.

Variational autoencoders consider that the latent space is indeed the
parameter space of $D$ Gaussian distributions, from which the input of the
decoder is sampled [7], allowing for non-sparse sampling and a quantification
of uncertainty in the point estimates. The $\beta$-VAE loss function [8] was
used (Equation 1):

$\hypertarget{eq-
vae}{}{\mathcal{L}_{VAE}=\mathcal{L}_{recon}+\beta\mathcal{L}_{KLD}}$ (1)

This loss function allows controlling for the proportion of the reconstruction
error and the Kullback-Leibler Divergence (KLD) (Equation 3) using the
parameter $\beta$. Reconstruction error was estimated assuming that the voxel
distribution among patients was guassian, which means that we could use the
volumetric Mean Squared Error (MSE) (Equation 2) between the original voxel
intensities $x_{j}$ and the output voxel of the decoder $\hat{x}_{j}$ for all
$j$ voxels in the images. Each subject’s loss was added to conform the batch
loss (reduction “sum”), which proved to improve training convergence.

$\hypertarget{eq-
recon}{}{\mathcal{L}_{recon}=\sum_{j}(x_{j}^{2}-\hat{x}_{j}^{2})}$ (2)
$\hypertarget{eq-
kld}{}{\mathcal{L}_{KLD}=-\frac{1}{2}\sum_{i=0}^{D}\left(1+\log\sigma_{i}^{2}-\mu_{i}^{2}-\sigma_{i}^{2}\right)}$
(3)

The specific architecture and parameters of the 3D-CVAE is shown at Table
LABEL:tbl-encoder-params and Table LABEL:tbl-decoder-params. Different latent
dimensionality (3, 8 and 20) was tested for 400 epochs with Adam optimizer
(lr=1e-3).

Table 1: Encoder architecture parameters Type # in # out KS Stride Padding Conv3D (1) 1 32 3 2 1 Conv3D (2) 32 64 3 2 1 Conv3D (3) 64 128 3 2 1 Conv3D (4) 128 256 3 2 1 Linear 73728 512 - - - 2 $\times$ Linear ($\mu$ and $\sigma$) 512 $D$ - - - Table 2: Decoder architecture parameters |  |  |  |  |
---|---|---|---|---|---
Type |  # in |  # out |  KS |  Stride |  Padding
Linear |  $D$ |  73728 |  – |  – |  –
ConvTranspose3d (1) |  256 |  128 |  3 |  2 |  1
ConvTranspose3d (2) |  128 |  64 |  3 |  2 |  1
ConvTranspose3d (3) |  64 |  32 |  3 |  2 |  1
ConvTranspose3d (4) |  32 |  1 |  3 |  2 |  1

###  2.3 Regression

The purpose of this work is to predict 4 categories of symptomatology using
the latent distribution of the FPCIT image dataset. To perform such task, we
explore the predictive power of two widely-used regression algorithms:
Decision Trees (DT) and XGBoost. DTs are relatively simple to understand and
have been widely used in medicine, because the resulting trees are directly
interpretable [9]. XGboost is a scalable tree boosting system that uses
approximate tree learning, and is widely extended in data science, achieving
high performance in many challenges [10].

###  2.4 Experimental Setup

The machine learning pipeline is defined by first training the dataset and
using the $\mu$ parameters of the latent distribution as features.
Additionally, K-Means features (KMF), which capture the distance to a K-means
cluster center, were generated to capture the salient characteristics of the
data. The number of clusters used to generate KMFs is established as the
maximum value between 8 and the logarithm of the dimensionality times 8, as in
[11]. Then, the two regression algorithms are trained with or without these
KMFs to predict the different UPDRS categories. 10-Fold cross-validation is
used to compute the performance values, including the Mean Absolute Error
(MAE), Root Mean Squared Error (RMSE) and coefficient of determination $R^{2}$
between the real and predicted values for UPDRS.

To more fully understand the output of the pipeline, we applied SHapley
Additive exPlanations (SHAP) [12] to the outputs of each system. SHAP yields
metrics such as feature importance, and the contribution of each feature to
the output of the algorithm, and thus can be used to track the most relevant
characteristics on the data.

##  3 Results and Discussion

###  3.1 Regression Results

The performance of the proposed regression algorithms on the 3, 8 and
20-dimensional latent space of the CVAE is shown at Figure 1. We obtained R2>0
for almost all target symptomatology categories and d-dimensional latent
spaces. We can observe that, while low latent dimensionality (3, 8) is -in
general- better for individual symptom categories, the details of the overall
symptomatology score (UPDRS – total) is better accounted for by the
20-dimensional (20-D) CVAE, achieving a MAE = 12.21 and a R2=0.26. This speaks
of the complexity of the composite, but also of its relationship to the FPCIT
distribution patterns. Similarly, UPDRS 3, more dependent on measured motor
symptoms, benefits from the 20-D space.

Figure 1: R2 performance for the models trained with the five symptomatology
scales and 3, 8 and 20 latent-dimension.

A detailed view of the best results obtained for each combination of target
scale and dimensionality of the latent space D, including their MAE, RMSE and
$R^{2}$ values is shown at Table LABEL:tbl-performance.

Table 3: Classification and regression results. DT stands for Decission Trees, XGB for XG-Boost and KMF for K-Means features target | d | model | MAE | RMSE | R2
---|---|---|---|---|---
| 3 | DT | 3.53 | 4.5 | 0.15
updrs1 | 8 | DT (KMF) | 3.58 | 4.55 | 0.13
| 20 | DT (KMF) | 3.69 | 4.71 | 0.07
| 3 | XGB (KMF) | 3.96 | 5.23 | 0.11
updrs2 | 8 | XGB (KMF) | 3.68 | 4.86 | 0.23
| 20 | XGB (KMF) | 3.74 | 5.08 | 0.16
| 3 | XGB (KMF) | 10.78 | 13.11 | 0.06
updrs3 | 8 | XGB (KMF) | 10.21 | 12.93 | 0.09
| 20 | XGB (KMF) | 10.03 | 12.59 | 0.14
| 3 | DT | 0.7 | 1.44 | 0
updrs4 | 8 | XGB | 0.62 | 1.39 | 0.07
| 20 | DT (KMF) | 0.71 | 1.43 | 0.02
| 3 | XGB (KMF) | 13.15 | 16.42 | 0.17
updrs (total) | 8 | XGB | 12.72 | 15.75 | 0.24
| 20 | XGB | 12.21 | 15.51 | 0.26

###  3.2 Intepretability of the Space

The regression results show that there exist a predictive power in the
manifold representation of the FPCIT dataset. The interpretation of the
effective patterns captured by latent variables can be, in consequence, of key
importance for clinical validation. While DTs (best for UPDRS 1) are visual by
nature, ensemble methods such as XGBoost are more difficult to interpret.
Here, SHAP can pave the way to interpret the contribution of each of the
latent features to the output variables.

In this regard, the prediction of total symptomatology (UPDRS - total) is of
special interest. In both the 8- and 20-dimensional spaces XGBoost achieves
the largest predictive ability in terms of $R^{2}$. Consequently, the most
important features -measured by SHAP- can capture the patterns that lead to
this high performance. The SHAP analysis of the XGBoost regressor with the
20-dimensional space reveals the behavior shown at Figure 2.

Figure 2: SHAP dependence plot of the XGBoost regressor with the
20-dimensional space

SHAP dependence plots show how the SHAP value of importance depends on a given
feature, ordered from highest to lowest importance. In Figure 2 we see that
the top-3 features that contribute to the output of the algorithm are
variables 2, 10 and 0. Color shows the values of a second feature that may
have an interaction effect with the feature we are plotting. Remember that
these are the mean parameters of a Gaussian distribution that encodes the
pattern for each subject. In this regard, we see an almost linear dependency
between importance to the algorithm output and values of the variables, which
indicates that there exist a clear relationship between them. To visualize the
patterns captured by features 2 and 10, we generated brain images sampling
these two features with the decoder and setting 0-mean to all other features.
The result is shown at Figure 3.

Figure 3: Latent manifold for variables 10 and 2 of the 20-dimensional latent
space of the 3D CVAE.

We observe that the composition of the two variables account for two relevant
characteristics of FPCIT: the overall intensity of the striata, the separation
between them and the ratio between uptake at the anterior and posterior parts
of the striata. General intensity and separation patterns seem to be encoded
by our variable 10, whereas the striata anterior-posterior length is generally
encoded by variable 2. In this regard, the overall intensity (uptake) at the
striata is indeed the objective of FPCIT as a biomarker, and therefore
expected. As for the anterior-posterior striata patterns, many works that used
FDOPA, a presynaptic dopaminergic marker, reported a distinct
anterior–posterior gradient of uptake as PD progresses [13]. This could be
then a relevant marker for assessing the progression of the disease. Since
SHAP importance is higher for this variable 2, we can assume that its
contribution to the modelling of the progression is higher, perhaps showing a
more smooth transition as the symptoms progress, in contrast to overall drug
uptake, more related to late neurodegeneration.

##  4 Conclusions

We hypothesize that the non-linear self-supervised latent space of a 3D
Convolutional Variational Autoencoder (CVAE) is linked to the symptoms shown
by affected subjects. The latent features of a trained CVAEs was related to
different aspects of the MDS-UPDRS scale with R2>0.20, proving a link between
FPCIT spatial patterns and Parkinson’s Disease symptomatology. The most
relevant variables for the predictive algorithms reveal that the CVAE captures
patterns related to overall intensity, striata separation and differences
between anterior-posterior parts of the striata, and that this last feature is
more relevant for assessing the progression of symptoms than the former. This
comprehensive approach may help to better understand the complex relationship
between clinical disease and imaging biomarkers, beyond the obvious
relationship between the average binding-ratio and the disease.

## Acknowledgements

This work was partly supported by the MINECO/FEDER under the
RTI2018-098913-B-I00 projects, and in by the Consejería de Economía,
Innovación, Ciencia y Empleo (Junta de Andalucía) and FEDER under the
P20-00525 project and PPJIA2021-17. Work by F.J.M.M. was supported by the
MICINN “Ramón y Cajal” RYC2021-030875-I. Work by C.J.M. is supported by
Ministerio de Universidades under the FPU18/04902 grant.

## References

re [1] V.L. Feigin, A.A. Abajobir, K.H. Abate, et al., “Global, regional, and
national burden of neurological disorders during 1990–2015: A systematic
analysis for the global burden of disease study 2015,”. _The Lancet
Neurology_. vol. 16, no. 11, pp. 877–897, 2017.

pre [2] F. Segovia, J.M. Górriz, J. Ramírez, F.J. Martínez-Murcia, and D.
Salas-Gonzalez, “Preprocessing of 18F-DMFP-PET Data Based on Hidden Markov
Random Fields and the Gaussian Distribution,”. _Frontiers in aging
neuroscience_. vol. 9, p. 326, 2017.

pre [3] E. Tolosa, A. Garrido, S.W. Scholz, and W. Poewe, “Challenges in the
diagnosis of parkinson’s disease,”. _The Lancet Neurology_. vol. 20, no. 5,
pp. 385–397, 2021.

pre [4] F.J. Martinez-Murcia, J.M. Górriz, J. Ramírez, and A. Ortiz,
“Convolutional neural networks for neuroimaging in parkinson’s disease: Is
preprocessing needed?”. _International journal of neural systems_. vol. 28,
no. 10, p. 1850035, 2018.

pre [5] F.J. Martinez-Murcia, A. Ortiz, J.-M. Gorriz, J. Ramirez, and D.
Castillo-Barnes, “Studying the Manifold Structure of Alzheimer’s Disease: A
Deep Learning Approach Using Convolutional Autoencoders,”. _IEEE Journal of
Biomedical and Health Informatics_. vol. 24, no. 1, pp. 17–26, 2020.

pre [6] C.G. Goetz, B.C. Tilley, S.R. Shaftman, et al., “Movement disorder
society-sponsored revision of the unified parkinson’s disease rating scale
(MDS-UPDRS): Scale presentation and clinimetric testing results,”. _Movement
disorders: official journal of the Movement Disorder Society_. vol. 23, no.
15, pp. 2129–2170, 2008.

pre [7] D.P. Kingma and M. Welling, “Auto-encoding variational bayes,”. _arXiv
preprint arXiv:1312.6114_. 2013.

pre [8] I. Higgins, L. Matthey, A. Pal, et al., “Beta-vae: Learning basic
visual concepts with a constrained variational framework,”. In: _International
conference on learning representations_ (2017).

pre [9] J.E. Arco, A. Ortiz, J. Ramírez, F.J. Martínez-Murcia, Y.-D. Zhang,
and J.M. Górriz, “Uncertainty-driven ensembles of multi-scale deep
architectures for image classification,”. _Information Fusion_. vol. 89, pp.
53–65, 2023.

pre [10] T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,”.
In: _Proceedings of the 22nd acm sigkdd international conference on knowledge
discovery and data mining_. pp. 785–794 (2016).

pre [11] C. Boutsidis, A. Zouzias, and P. Drineas, “Random projections for
$k$-means clustering,”. _Advances in neural information processing systems_.
vol. 23, 2010.

pre [12] S.M. Lundberg and S.-I. Lee, “A unified approach to interpreting
model predictions,”. In: I. Guyon, U.V. Luxburg, S. Bengio, et al., Eds.
_Advances in neural information processing systems_. _Curran Associates, Inc._
(2017).

pre [13] V. Dhawan, Y. Ma, V. Pillai, et al., “Comparative analysis of
striatal FDOPA uptake in parkinson’s disease: Ratio method versus graphical
approach,”. _Journal of Nuclear Medicine_. vol. 43, no. 10, pp. 1324–1330,
2002.

p