<figure><img src="image_1.png" /></figure>

HAL
open science

Benchmarking 3D generative autoencoders for
pseudo-healthy reconstruction of brain
18F-fluorodeoxyglucose positron emission tomography

Ravi Hassanaly, Maëlys Solal, Olivier Colliot, Ninon Burgos

► To cite this version:

Ravi Hassanaly, Maëlys Solal, Olivier Colliot, Ninon Burgos. Benchmarking 3D generative autoencoders for pseudo-healthy reconstruction of brain 18F-fluorodeoxyglucose positron emission tomography. Journal of Medical Imaging, 2025, 12 (5), pp.054005. 10.1117/1.JMI.12.5.054005. hal-04445378v2

HAL Id: hal-04445378

https://inria.hal.science/hal-04445378v2

Submitted on 16 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Benchmarking 3D generative autoencoders for pseudo-healthy reconstruction of brain
$^{18}$F-fluorodeoxyglucose positron emission tomography

Ravi Hassanaly<sup>a</sup>, Maëlys Solal<sup>a</sup>, Olivier Colliot<sup>a</sup>, Ninon Burgos<sup>a, *</sup>,
for the Alzheimer's Disease Neuroimaging Initiative<sup>†</sup>

<sup>a</sup>Sorbonne Université, Institut du Cerveau - Paris Brain Institute - ICM, CNRS, Inria, Inserm, AP-HP,
Hôpital de la Pitié-Salpêtrière, F-75013, Paris, France

**Abstract.** **Purpose** Many deep generative models have been proposed to reconstruct pseudo-healthy images for anomaly detection. Among these models, the variational autoencoder (VAE) has emerged as both simple and efficient. While significant progress has been made in refining the VAE within the field of computer vision, these advancements have not been extensively applied to medical imaging applications.

**Approach** We present a benchmark that assesses the ability of multiple VAEs to reconstruct pseudo-healthy neuroimages for anomaly detection in the context of dementia. We first propose a rigorous methodology to define the optimal architecture of the vanilla VAE and select through random searches the best hyper-parameters of the VAE variants. Relying on a simulation-based evaluation framework, we thoroughly assess the ability of 20 VAE models to reconstruct pseudo-healthy images for the detection of dementia-related anomalies in 3D brain $^{18}$F-fluorodeoxyglucose (FDG) positron emission tomography (PET) and compare their performance.

**Results** This benchmark demonstrated that the majority of the VAE models tested were able to reconstruct images of good quality and generate healthy looking images from simulated images presenting anomalies.

**Conclusions** Even if no model clearly outperformed all the others, the benchmark allowed identifying a few models that perform slightly better than the vanilla VAE. It further showed that many VAE-based models can generalize to the detection of anomalies of various intensities, shapes and locations in 3D brain FDG PET.

**Keywords:** Variational autoencoder, Unsupervised anomaly detection, Deep generative models, PET, Alzheimer's disease.

*Ninon Burgos, ninon.burgos@cnrs.fr

# 1 Introduction

The synergy between innovations in imaging technologies, the growing volume of medical data, and sophisticated machine learning algorithms have given rise to algorithms capable of performing complex tasks such as anomaly detection for computer-aided diagnosis.¹ A first approach consists in using supervised learning algorithms to segment lesions. Whilst this approach performs well, it requires human annotations that are costly in time and money, and leads to very specialized models. Another strategy consists in using unsupervised algorithms to learn the distribution of healthy data and to detect anomalies as deviations from this learned prior knowledge.² This is particularly useful as it does not require annotated data, and should generalize to any type of anomalies, without having seen them before.

A prevailing approach in unsupervised anomaly detection (UAD) consists in using generative models to reconstruct healthy looking images, or pseudo-healthy images.¹⁻³ The pseudo-healthy image can then be compared with the real image to detect anomalies, and possibly further

<footnote>†Data used in preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf.</footnote>
segment lesions. The underlying assumption is that a model trained on images of subjects diagnosed as healthy should generate images that do not contain pathology-specific features and look like a healthy image, even when a pathological image is provided as input.

Pseudo-healthy reconstruction approaches that have been developed for medical imaging often rely on generative models such as variational autoencoders (VAEs),⁴ generative adversarial networks (GANs)⁵ and more recently diffusion models.⁶ Although diffusion models have shown remarkable performance for image generation, they do not easily scale to 3D images,⁷ mainly due to memory issues. On the GAN side, after the foundational work of Schlegl et al., AnoGAN⁸ and f-AnoGAN,⁹ only a few works have been published. They either use cycle GANs¹⁰ or combine GANs with autoencoders.¹¹,¹² On the other hand, even though VAEs' image generation quality is lower, they are easy to train, they scale well to high-dimensional data, provide good interpretation capacity thanks to their regularized latent space, and are able to handle small datasets. Many new VAE extensions have shown their efficacy in the computer vision literature,¹³⁻²⁹ but only a handful have been applied to medical imaging.³,³⁰⁻³³

In 2021, Baur et al.³ compared VAE-based approaches to the best GANs for unsupervised anomaly segmentation in brain structural magnetic resonance imaging (MRI). It was conducted on models that had already been employed for UAD in the medical imaging context, such as Context VAE,³⁴ Constrained AAE,³⁰ or AnoVAEGAN.³⁵ They showed that the vanilla VAE used for density-based restoration³⁶ outperforms other models, including GAN approaches, at the cost of a longer inference time. They compared the performance of their models using segmentation metrics such as the dice similarity coefficient (DSC) and the area under the precision-recall curve (AUPRC) computed between the residual (i.e. the difference between the input and the reconstructed image) and the ground truth anomaly mask provided in the datasets they used. This study focused on the segmentation of glioblastoma and multiple sclerosis lesions, which consist of sharp and intense anomalies that are segmented in 2D slices extracted from MRI volumes.

In this work, our objective is to apply UAD methods to identify metabolic changes associated with Alzheimer's disease and other dementias³⁷ that are visible in brain ¹⁸F-fluorodeoxyglucose (FDG) positron emission tomography (PET) images. ¹⁸F-FDG PET images are 3D images that highlight the concentration of administered FDG, a tracer used to localize hypometabolism in the case of neurodegeneration.³⁸ This application is particularly interesting as deep learning methods for UAD have rarely been applied for the diagnosis of dementia,³¹,³⁹,⁴⁰ whereas this approach could enable early diagnosis since changes visible in neuroimaging can occur years before the onset of initial symptoms.⁴¹ Metabolic abnormalities can be difficult to detect as they are diffuse and sometimes subtle (limited difference in intensity between normal tissues and areas with hypometabolism),⁴² contrary to glioblastoma or white matter hyper-intensities usually studied on structural MRI.³,¹⁰,³⁰,³⁴

Following the work of Baur et al.,³ we introduce a methodology to compare generative models for pseudo-healthy reconstruction and apply it to benchmark 20 VAE-based models, focusing on the pseudo-healthy reconstruction of 3D FDG PET images for anomaly detection in the context of dementia. We compare many VAE-based models that have not been applied to medical image analysis yet, thanks to the software package of Chadebec et al.⁴³ In contrast to computer vision works, where datasets typically contain several tens of thousands images, it will be interesting to examine the performance of such models when trained on a relatively small dataset, comprising only a few hundred images, which is typical in medical imaging. Our contributions are threefold:

1. first, we propose a rigorous methodology, and provide the associated software tool, that we used to both define the optimal architecture of the vanilla VAE and select the best hyper-parameters of the VAE variants in the context of neuroimaging. It consists of a
first random-search on the architectural parameters of the VAE that fixes the encoder and decoder structures for all the models, followed by a random hyper-parameter search to find optimal ones for each model;

2. then, we put in application the evaluation framework introduced by Hassanaly et al.⁴⁴ to thoroughly assess the ability of 20 VAE models to reconstruct pseudo-healthy images for the detection of dementia-related anomalies in 3D brain FDG PET and compare their performance;

3. finally, we conclude on the best performing models, providing a state-of-the-art on the use of 3D convolutional VAEs in such context.

A preliminary version of this work was published as a conference paper.⁴⁵ The present article extends the previous work mainly with: (i) the addition of three new VAE-based models; (ii) the optimization of the best encoder-decoder architecture for the vanilla VAE through an extensive random search; (iii) the optimization of the hyper-parameters for all the models based on literature results and extensive random search; (iv) the use of full resolution 3D brain FDG PET instead of downsampled images; (v) and the thorough evaluation of the different models using the evaluation protocol described in Hassanaly et al.⁴⁴

# 2 Pseudo-healthy image reconstruction with variational autoencoders

## 2.1 Variational autoencoder framework

The VAE framework⁴ assumes that a latent variable z is involved in the generation process of the input data x: $p_{\theta}(x) = \int_{z} p(z)p_{\theta}(x|z)dz$ where z ~ $p_{\theta}(z)$ is the prior distribution on the latent space and $p_{\theta}(x|z)$ is the generative model (or the decoder). To compute the appropriate z for each data input x of the dataset, the posterior distribution $p_{\theta}(z|x)$ needs to be modeled. Since it is untractable, the posterior distribution is approximated using variational inference by introducing another model $q_{\phi}(z|x)$ such that $q_{\phi}(z|x) \approx p_{\theta}(z|x)$. $q_{\phi}(z|x)$ is the inference model (or encoder). Both the decoder and encoder are parametric models whose parameters $\theta$ and $\phi$ are given by a neural network.

The objective is to maximize the likelihood of $p_{\theta}(x)$, which is equivalent to maximizing the evidence lower bound, which defines the loss function $\mathcal{L}_{\theta,\phi}$

$$ \log(p_{\theta}(\mathbf{x})) \geq \mathcal{L}_{\theta,\phi}(\mathbf{x}) = \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})} \left[ \log(p_{\theta}(\mathbf{x}|\mathbf{z})) \right] - D_{\text{KL}}(q_{\phi}(\mathbf{z}|\mathbf{x})||p_{\theta}(\mathbf{z})) \quad (1) $$

with $D_{\text{KL}}$ the Kullback-Leibler divergence.

The VAE is particularly suited for pseudo-healthy reconstruction. Let's consider D a set of medical images collected following a similar protocol. D contains healthy and pathological images and is the union of two complementary subsets $D_h$ and $D_p$. For instance, $D_h$ could be a set of healthy FDG PET images $x \in D$ whose distribution is $p(x)$. The goal of pseudo-healthy image reconstruction is to reconstruct an FDG PET image of healthy appearance given an input $x \in D$. During the training process, an approximation of the posterior distribution $q_{\phi}(z|x)$ is learned for $x \in D_h$ as the model is trained using only healthy subjects. In other words, the healthy image true distribution $p(x)$ is approximated with the learned parametric distribution $p_{\theta}(x)$ such that $p_{\theta}(x) \approx p(x)$. During reconstruction, this approximate posterior is used to estimate the latent variable z for $x \in D$ (it can be from $D_h$ or $D_p$), i.e., the images (of healthy subjects or patients) are projected into that “healthy images” learned subspace. Then, the decoder can generate healthy images from z.
## 2.2 Extensions to the variational autoencoder framework

Several contributions have been proposed to improve the VAE framework.43 These contributions can be divided into four categories that correspond to different objectives.

The aim of the first category of approaches is to improve the prior distribution $p(z)$ by using a variational mixture of posteriors as prior (VAMP),27 by using a specific geometry in the latent space such as hyperspherical VAE (SVAE),17 by learning the prior on a discrete latent space with vector quantized-VAE (VQVAE),28 or by substituting the prior with a density estimation method using regularization with a gradient penalty (RAE-GP) or an $\ell^2$ penalty on the decoder (RAE-$\ell^2$).18

Other methods aim to better estimate the lower bound by using importance weighting (IWAE),13 by using linear normalizing flows (VAE LinNF),24 inverse autoregressive flows (VAE-IAF)21 or Markov chain Monte Carlo using Hamiltonian importance sampling (HVAE)15 to better estimate the posterior.

Approaches in the third category encourage disentanglement of the features in the latent space by adding a weight to balance the terms of the loss in Eq. 1 ($\beta$-VAE),19 subsequently improved with a better reconstruction capacity by progressively increasing the KL-divergence term (Disentangled $\beta$-VAE),14 by decomposing the loss to show a total correlation term ($\beta$-TC VAE),16 or by encouraging the distribution of the latent variable $q(z)$ to be factorial (FactorVAE).20

Finally, other methods change the distance computed between the distributions by adding the mutual information between x and z as regularization (InfoVAE),29 using another divergence term in the loss such as the maximum mean discrepancy in the Wasserstein autoencoder (WAE)26 or a discriminator to differentiate a prior's sample from a posterior's sample in the adversarial autoencoder (Adv. AE),23 or by changing the reconstruction metric for another similarity metric such as the multi-scale structural similarity (MS-SSIM VAE),25 or for the prediction of a discriminator on the output of the VAE (VAEGAN).22

All of these models, described in more detail in Appendix D, perform well in computer vision, as shown by Chadebec et al.43 who compared 19 of them on classic computer vision datasets (MNIST, CIFAR10 and CELEBA) on five tasks: image reconstruction, image generation, classification, clustering and interpolation. However, they have not been compared in the context of medical imaging.

## 2.3 Evaluation, comparison and selection of the models

When evaluating unsupervised anomaly detection approaches, two aspects are usually assessed: their ability to reconstruct images of high quality and their ability to detect anomalies. The first aspect is assessed when reconstructing images of healthy subjects. Commonly used metrics are the mean-squared error (MSE), the peak signal-to-noise ratio (PSNR) and the structural similarity (SSIM).46 These paired metrics are computed between the input and reconstructed images.47 To assess the second aspect, most studies3, 10, 12, 30, 34, 48–55 rely on ground truth anomaly masks, making such evaluation similar to that of supervised segmentation methods using metrics such as the DSC or the AUPRC.3 To go further, Xia et al.10 defined two new metrics: an identity metric to measure whether the image reconstructed corresponds to the same subject as the image given in input, and a healthiness score to measure whether the reconstructed image appears healthy. However, these metrics use the a priori information of the anomaly mask, which limits their application. In situations when the ground truth anomaly mask is unavailable, the pseudo-healthy reconstructions can be evaluated using a classifier trained to differentiate pathological and healthy images,42 or using an anomaly score derived from the reconstruction error (supposed to be low for healthy images and high for pathological images).31, 54
In our context, the detection of hypometabolism in brain FDG PET images, no ground truth anomaly masks are usually available. We thus previously proposed a simulation framework<sup>44,56</sup> that consists in simulating the effect of a disease on images of healthy subjects by reducing the PET uptake within areas of the brain associated with different dementias<sup>42</sup> defined using a mask *M*. After locally reducing the intensity of the image by a certain percentage within *M*, a Gaussian smoothing is applied to generate a realistic result and diffuse anomalies. This approach effectively replicates realistic regional hypometabolism and provides pairs of diseased images with the original healthy scan that is used as the target ground truth for the pseudo-healthy reconstruction. Hypometabolism can be simulated in different parts of the brain, which allows assessing the models' ability to generalize to anomalies caused by different dementia subtypes. Furthermore, for a given dementia subtype, different severity degrees can be simulated by adjusting the intensity reduction. This allows investigating the sensitivity of the UAD approaches to both subtle and severe anomalies. All the details about the masks used and the simulation pipeline are in our previous works.<sup>44,56</sup>

In Hassanaly et al.,<sup>44</sup> we also proposed a metric, which relies on the simulation framework, to evaluate whether a model is able to reconstruct images that are looking healthy. This healthiness score *H* is defined as

$$H = \frac{\mu_M}{\mu_{\bar{M}}}, \qquad (2)$$

with $\mu_M$ the average uptake within the mask *M* used to simulate the anomaly and $\mu_{\bar{M}}$ the average uptake of voxels in the brain excluding the mask *M*. The healthiness score compares the average uptake in the region in which we simulate the disease and the other regions of the brain. It is supposed to be around 1 for images of healthy subjects, lower than 1 for simulated images, and expected to be around 1 again for the pseudo-healthy reconstructions.

As we consider that to accurately detect anomalies a model should reconstruct pseudo-healthy images of high quality, we use the pairwise performance measures as a first step in our evaluation. We especially rely on the SSIM, rather than the MSE or PSNR, as it is a perceptual metric that appears more informative than a pixel-wise difference, and because it is a different metric than the optimization criterion, which is MSE for all the models except for the MS-SSIM VAE.<sup>25</sup> In particular, the SSIM is used as selection criterion when searching for the best hyper-parameters' configurations and selecting the best trained models, and is combined with the MSE when searching for the best encoder-decoder architecture. We use the simulation framework with the healthiness score in a second step to push further the evaluation of the trained models being compared.

## 3 Materials

In this study, we use FDG PET scan from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database.<sup>57,58</sup> The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial MRI, PET, other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment and early AD.

We selected FDG PET images co-registered, averaged and uniformized to a resolution of 8 mm full width at half maximum to reduce the variability due to the use of different scanners. We only used PET images for which a T1-weighted (T1w) MR image was available at the same session for preprocessing purposes. The images were then processed using Clinica's<sup>59</sup> pet-linear pipeline: they were registered using a rigid transformation to the corresponding T1w MRI of the same session, and then affinely registered to the MNI ICBM 2009c Nonlinear
<table><caption>**Table 1** Summary of participant demographics at baseline for the different training/validation splits and test sets considered. Note that split *s* corresponds to using fold *s* from the 6-fold cross-validation as validation set and the other folds as training set.</caption><thead><tr><th></th><th>Set</th><th># subjects (%F)</th><th># images</th><th>Avg age (± SD)</th><th>Age range</th></tr></thead><tbody><tr><td rowspan="2">split 0</td><td>train</td><td>265 (52.8%)</td><td>536</td><td>74.9 ± 6.1</td><td>55.8, 95.0</td></tr><tr><td>validation</td><td>53 (41.5%)</td><td>53</td><td>72.8 ± 5.5</td><td>62.3, 85.3</td></tr><tr><td rowspan="2">split 1</td><td>train</td><td>265 (53.5%)</td><td>533</td><td>74.7 ± 6.1</td><td>55.8, 95.0</td></tr><tr><td>validation</td><td>53 (37.7%)</td><td>53</td><td>73.4 ± 5.6</td><td>59.9, 88.9</td></tr><tr><td rowspan="2">split 2</td><td>train</td><td>265 (49.8%)</td><td>537</td><td>74.6 ± 5.9</td><td>55.8, 95.0</td></tr><tr><td>validation</td><td>53 (56.6%)</td><td>53</td><td>74.4 ± 6.6</td><td>63.8, 93.6</td></tr><tr><td rowspan="2">split 3</td><td>train</td><td>265 (48.6%)</td><td>538</td><td>74.7 ± 6.0</td><td>55.8, 95.0</td></tr><tr><td>validation</td><td>53 (62.2%)</td><td>53</td><td>73.9 ± 6.2</td><td>61.2, 86.2</td></tr><tr><td rowspan="2">split 4</td><td>train</td><td>265 (49.1%)</td><td>511</td><td>75.0 ± 5.9</td><td>59.7, 95.0</td></tr><tr><td>validation</td><td>53 (60.4%)</td><td>53</td><td>72.4 ± 6.5</td><td>55.8, 84.7</td></tr><tr><td rowspan="2">split 5</td><td>train</td><td>265 (51.6%)</td><td>510</td><td>74.6 ± 5.9</td><td>55.8, 93.6</td></tr><tr><td>validation</td><td>53 (47.2%)</td><td>53</td><td>73.1 ± 6.6</td><td>59.7, 92.8</td></tr><tr><td rowspan="2">Test</td><td>CN test</td><td>60 (63.3%)</td><td>60</td><td>73.5 ± 6.6</td><td>59.8, 85.8</td></tr><tr><td>AD test</td><td>353 (41.1%)</td><td>353</td><td>75.3 ± 7.6</td><td>55.1, 90.3</td></tr></tbody></table>

Symmetric template⁶⁰,⁶¹ using the transformation computed with the t1-linear pipeline. They were then normalized in intensity using the average PET uptake in a region comprising cerebellum and pons, and cropped. In the end, the dimension of the PET scan is 169 × 208 × 179 with 1 mm isotropic voxels. To filter out potential PET images not well registered to the MNI template, we used the quality control pipeline of ClinicaDL.⁶²

A total of 3511 FDG PET scans from 1600 participants were downloaded from the ADNI database. This includes 554 cognitively normal (CN) subjects (1010 images) that we selected since UAD models are trained only on images from healthy subjects. We know that physiological changes can appear several years before the first clinical symptoms, so to ensure that images really correspond to a healthy brain, we kept only scans from subjects that are CN for at least three years after the session considered. All the details about data selection are in Appendix A. Our final dataset comprises 739 images from 378 CN subjects.

We split our dataset of 378 CN subjects into training/validation and test sets at the subject's level to avoid any form of data leakage,⁶³ stratifying by sex and age to reduce biases. Only baseline sessions were kept in the test set to avoid biased results. The test set, comprising 60 CN subjects (60 images), is used to assess whether the healthy images are reconstructed as healthy. We denote it as "Test CN". We then performed a six-fold cross-validation on the training/validation data to estimate the variance due to data splitting.⁶⁴ 53 subjects (53 images) belong to the validation sets to monitor the training and 265 subjects (between 510 and 538 images depending on the fold) are used to train our models. The training split and test set statistics are summarized in Table 1.

We also use the 60 images from the CN test set to build new test sets by using our simulation method.⁴⁴ We simulate AD with nine intensity levels from 5% to 70% hypometabolism and simulate six dementia subtypes at 30% hypometabolism: AD, behavioral variant frontotemporal dementia (bvFTD), logopenic variant primary progressive aphasia (lvPPA), semantic variant PPA (svPPA), non-fluent variant PPA (nfvPPA), and posterior cortical atrophy (PCA), which results in a total of 14 simulated test sets. These test sets are denoted using the dementia simulated and
the hypometabolism intensity. For instance, “Test AD 30” corresponds to images simulating AD with a 30% hypometabolism.

In addition to the images of CN subjects, the data downloaded from the ADNI dataset include FDG PET images of AD patients. After applying the same data selection procedure as the one used for the CN subjects (explained in Appendix A), we keep 353 images at baseline from 353 AD patients for testing purposes.

# 4 Model selection: searching for the best architecture and hyper-parameters

We present a methodology to compare 20 AE and VAE-based models. For the comparison to be meaningful, we defined a protocol to find the best architectures and parameters for each model. We decided to use the same encoder-decoder architecture for all the models as it would have been too long to find an optimal architecture for each model, and as we believe it makes the comparison fairer. The architecture was obtained using a random search on the vanilla VAE. After having fixed the encoder and decoder architecture, the latent space dimension and some training parameters, we attempted to find optimal hyper-parameters for each variant through either a random search or a grid search. Once the best parameters for each model were found, we trained them on all the six splits of the cross-validation, and we selected the best. The best trained models were finally evaluated using the simulation framework and metric presented in Section 2.3. The procedure is summarized in Figure 1. The random search and evaluation procedures are implemented in ClinicaDL⁶² while the VAE-based models are implemented in Pythae,⁴³ which are both open-source tools (clinicadl.readthedocs.io, pythae.readthedocs.io).

<figure><img src="image_3.png" /><figcaption>**Fig 1** Diagram summarizing the benchmark steps. We represent steps performed on training sets in purple (random search and training), the selections on the validation sets are represented in blue and the final evaluation on test sets is in green.</figcaption></figure>

All the models were trained for 200 epochs on an HPC with Nvidia Tesla V100 GPUs that have 32 GB of dedicated memory. The choice for the batch size and the learning rate will be discussed further in this section. We used the same environment to train all the models.

## 4.1 Selection of the encoder-decoder architecture

The training parameters and the encoder-decoder architecture were selected with a random search for the vanilla VAE. We trained the models on a random selection of three splits as a trade-off between reducing the variance due to data splitting and the computational time required to train the models. We then selected the models based on the average SSIM and MSE, computed within the full image field of view, on the validation sets.

We defined a modular architecture for the encoder and decoder which is shown in Figure 2. The encoder (shown in Figure 2a) is composed of a number $B_e$ of blocks, each containing a number $S_e$ of sub-blocks. Similarly, the decoder is composed of a number $B_d$ of blocks, each containing a number $S_d$ of sub-blocks (Figure 2b). For a chosen architecture, the sub-blocks can either all be convolutional or all be residual (see Figure 2c). In both cases, the convolution layers
<figure><img src="image_4.png" /><figcaption>(a) Encoder architecture with $B_e$ blocks and $S_e$ sub-blocks. Blocks and sub-blocks in dashed lines are optional.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Decoder architecture with $B_d$ blocks and $S_d$ sub-blocks. Blocks and sub-blocks in dashed lines are optional, and layers in dotted lines are only included if `last Conv = True`.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Sub-block architecture. Left: convolutional sub-blocks. Right: residual sub-blocks. Black sub-blocks: image size is divided by 2 along each dimension by using 4x4x4 convolutions with stride 2 (layers in darker gray), used as the first sub-block of each encoder block. Gray sub-blocks: image size remains the same by using 3x3x3 convolutions with stride 1 (layers in light gray). BN: batch normalization, SiLU: swish activation function</figcaption></figure>

**Fig 2** Encoder-decoder modular architecture. The number of convolution kernels in each sub-block is indicated under the sub-block (e.g. $C_e$).

are followed by a batch normalization and we use a swish activation function as suggested in Vahdat et al.⁶⁵

In the encoder, the number of channels is doubled by the first convolution in each block. At the same time, the size of the image is divided by 2 along each dimension by using a 3D convolution with kernel size (4, 4, 4), stride (2, 2, 2) and padding (1, 1, 1). The following sub-blocks are optional and their convolution operations have kernel size (3, 3, 3), stride (1, 1, 1) and padding (1, 1, 1). In the decoder, the last sub-block of each block is preceded by an upsampling layer, to be roughly symmetrical with the encoder. Convolution operations in the decoder have kernel size (3, 3, 3), stride (1, 1, 1) and padding (1, 1, 1). This architecture was inspired from ResNet models⁶⁶ and VGG models.⁶⁷

The parameters of this modular architecture (summarized in Table 2) are therefore the following: the latent space size, the number of blocks in the encoder $B_e$, the number of blocks in the decoder $B_d$, the number of sub-blocks per encoder block $S_e$, the number of sub-blocks per decoder block $S_d$, the number of channels for the first encoder block $C_e$, the number of channels for the last decoder block $C_d$, and the layer type (convolution or residual). We implement a random search to explore this parameter space, and choose possible values for each parameter
<table><caption>**Table 2** Hyper-parameters included in our encoder-decoder VAE architecture random search</caption><thead><tr><th>Hyper-parameter</th><th>Label</th><th>Search space</th><th>Selected value</th></tr></thead><tbody><tr><td>Number of encoder blocks</td><td>B<sub>e</sub></td><td>{4, 5, 6}</td><td>5</td></tr><tr><td>Number of sub-blocks per encoder block</td><td>S<sub>e</sub></td><td>{1, 2, 3}</td><td>1</td></tr><tr><td>Number of channels for the first encoder sub-block</td><td>C<sub>e</sub></td><td>{16, 32}</td><td>16</td></tr><tr><td>Number of decoder blocks</td><td>B<sub>d</sub></td><td>{4, 5, 6}</td><td>5</td></tr><tr><td>Number of decoder sub-blocks</td><td>S<sub>d</sub></td><td>{1, 2, 3}</td><td>1</td></tr><tr><td>Number of channels for the last decoder sub-block</td><td>C<sub>d</sub></td><td>{16, 32}</td><td>16</td></tr><tr><td>Latent space size</td><td></td><td>{256, 512, 1024}</td><td>256</td></tr><tr><td>Learning rate</td><td></td><td>{10<sup>-3</sup>, 10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>10<sup>-4</sup></td></tr><tr><td>Block type</td><td></td><td>{conv, res}</td><td>conv</td></tr><tr><td>Added convolution in last decoder block</td><td>lastoj conv</td><td>{True, False}</td><td>False</td></tr></tbody></table>

based on previous experiments, intermediate results (as we launched the random search in successive batches) and intuition. We decided to set the batch size of the data loader to 8. Even though this is a constraint for configurations that would require more memory, this choice allows flexibility; in scenarios where certain VAE variants require more memory, we can reduce the batch size while maintaining a reasonable number of images per batch (e.g., 6 or 4). Details of all the parameters tested and their impact are discussed in Appendix B.

After comparing around 200 configurations, the encoder architecture selected is composed of five blocks, each with one sub-block, each containing a convolutional layer, a batch normalization and a swish activation function. These blocks are followed by a flatten and a fully connected layer. The latent space has size 256. The decoder is symmetrical, it is composed of a fully connected layer followed by five blocks, each with one sub-block, each composed of an upsampling layer, a convolutional layer, a batch normalization and a swish activation. This model has 16 channels after the first encoder block and before the last decoder block. This is shown in Figure 3 and detailed in Table 9.

<figure><img src="image_7.png" /><figcaption>**Fig 3** Diagram of the selected VAE architecture</figcaption></figure>
## 4.2 Selection of the models' hyper-parameters

Once we found an encoder-decoder architecture that gave good performance, we used it for the AE and 18 VAE variants presented in Section 2.2. However, all of these variants, except the SVAE,¹⁷ have supplementary hyper-parameters that may have significant impact on the models' performance. We therefore searched for the best configuration of hyper-parameters for each model in the context of 3D brain FDG PET reconstruction by launching either a random search (when we searched for more than one hyper-parameter) or grid search (when there is only one hyper-parameter). Similarly to the architecture search, we train each configuration on three splits. We then selected the best set of hyper-parameters for each VAE-based model using the best average SSIM on the validation set as criterion.

As there were many models to optimize, we limited the number of random searches to $N \times 10$ with $N$ the number of parameters to search. For instance, when there was only one hyper-parameter to tune, we launched a grid search with maximum 10 different values for that parameter, if there were two parameters, we trained a maximum of 20 models and so on. This may not be the fairest decision as it does not allow exploring the same percentage of the parameter space depending on $N$ (as it scales to the power $N$ and not linearly), but it accounts for the fact that a model with more parameters may be more tedious to tune. Moreover, we carefully chose a range of values to test for each hyper-parameter of each model based on the original implementation papers, our prior knowledge of the models, and the work done by Chadebec et al.⁴³ Note that some of the hyper-parameters were excluded from our search when an optimal value was provided in the literature, which allowed reducing the number of configurations trained.

Following the vanilla VAE training, the different configurations were trained by default with a batch size of 8 and a learning rate of $10^{-4}$ on 200 epochs. When some set of hyper-parameters led to memory errors, we gradually reduced the batch size to 6, 4 or 2, and when they lead to errors in the computation of the loss, we reduced the learning rate to $10^{-5}$. In spite of this, combinations leading to errors were removed, further reducing the size of the hyper-parameter space.

The details relating to the different VAE-based models, their hyper-parameters, the random search, the trained configurations and the results are provided in Appendix D. A summary is proposed in Table 3.

We report in Table 4 reconstruction metrics for the 18 VAE variants with the best configuration of hyper-parameters that we tested. Out of all the models, only three did not perform well on the validation set (highlighted in gray): the VAMP²⁷ with an average SSIM of 0.702, the MS-SSIM VAE²⁵ with an average SSIM of 0.472 and the SVAE with a very low average SSIM of 0.151. We found it quite surprising that the MS-SSIM VAE²⁵ performed so poorly in terms of average SSIM, since it optimizes a perceptual metric related to the SSIM, namely the multi-scale SSIM.⁶⁸ These results could be explained by the fact that the MS-SSIM computation in 3D is very costly, meaning that the only kernel size that allowed training in a reasonable amount of time was 2, potentially leading to a poor estimation of the metric, especially since the kernel size suggested in the MS-SSIM original implementation is 11.⁶⁸ In the end, only three MS-SSIM VAE models with three different combinations of parameters were trained successfully, possibly explaining why we did not find a configuration giving acceptable reconstruction. Regarding the SVAE,¹⁷ we observed that no model could be trained with a latent space dimension of 64 or higher. This limitation is due to the computation of the modified Bessel function in the KL-divergence term of the loss, which is highly unstable, particularly for high-dimensional latent spaces.⁶⁹,⁷⁰ This instability not only explains why the model struggles to train with larger latent spaces but also
<table><caption>**Table 3** Summary of the hyper-parameters random search for each VAE variant. The hyper-parameters are detailed in Appendix D.</caption><thead><tr><th>Models</th><th>Hyper-parameters</th><th>Search space</th><th>Selected value</th></tr></thead><tbody><tr><td>Adv. AE<sup>23</sup></td><td>adv. loss scale</td><td>{0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99}</td><td>0.9</td></tr><tr><td rowspan="3">β-TC VAE<sup>16</sup></td><td>β</td><td>{0.001, 0.005, 0.01, 0.05, 0.1, 1, 2, 5, 10}</td><td>2</td></tr><tr><td>α</td><td>{1, 3}</td><td>1</td></tr><tr><td>γ</td><td>{1, 3}</td><td>1</td></tr><tr><td>β-VAE<sup>19</sup></td><td>β</td><td>{0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 2, 5, 10, 100}</td><td>10</td></tr><tr><td rowspan="3">Dis. β-VAE<sup>14</sup></td><td>β</td><td>{0.01, 0.1, 1, 5, 10}</td><td>10</td></tr><tr><td>C</td><td>{5, 25, 50}</td><td>50</td></tr><tr><td>warm-up epochs</td><td>{100, 1000}</td><td>1000</td></tr><tr><td>FactorVAE<sup>20</sup></td><td>γ</td><td>{2, 5, 10, 15, 20, 30, 40, 50, 100, 200}</td><td>40</td></tr><tr><td rowspan="3">HVAE<sup>15</sup></td><td>n<sub>lf</sub></td><td>{1, 2, 10, 15, 20}</td><td>10</td></tr><tr><td>ε<sub>lf</sub></td><td>{0.00001, 0.0001, 0.001, 0.01}</td><td>0.00001</td></tr><tr><td>β<sub>0</sub></td><td>{0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}</td><td>0.8</td></tr><tr><td rowspan="3">InfoVAE<sup>29</sup></td><td>kernel choice</td><td>{rbf, imq}</td><td>rbf</td></tr><tr><td>α</td><td>{0.0, 0.2, 0.4, 0.6, 0.8, 1.0}</td><td>1</td></tr><tr><td>λ</td><td>{0.01, 0.1, 1, 10, 100}</td><td>0.1</td></tr><tr><td rowspan="2">IWAE<sup>13</sup></td><td>kernel bandwidth</td><td>{0.01, 0.1, 0.5, 1, 5, 10, 100}</td><td>0.1</td></tr><tr><td>n samples</td><td>{2, 3, 4, 5, 6, 8, 10, 12, 15, 20}</td><td>6</td></tr><tr><td rowspan="2">MS-SSIM VAE<sup>25</sup></td><td>β</td><td>{0.01, 0.1, 1, 10, 100}</td><td>100</td></tr><tr><td>window size</td><td>{2, 3, 5, 11}</td><td>2</td></tr><tr><td rowspan="2">RAE-ℓ<sup>218</sup></td><td>embedding weight</td><td>{0.00001, 0.0001, 0.001, 0.01, 0.1, 1}</td><td>0.0001</td></tr><tr><td>reg. weight</td><td>{0.00001, 0.0001, 0.001, 0.01, 0.1, 1}</td><td>1</td></tr><tr><td rowspan="2">RAE-GP<sup>18</sup></td><td>embedding weight</td><td>{0.00001, 0.0001, 0.001, 0.01, 0.1, 1}</td><td>0.01</td></tr><tr><td>reg. weight</td><td>{0.00001, 0.0001, 0.001, 0.01, 0.1, 1}</td><td>0.0001</td></tr><tr><td>SVAE<sup>17</sup></td><td>latent space size</td><td>{8, 16, 32}</td><td>16</td></tr><tr><td rowspan="2">VAEGAN<sup>22</sup></td><td>adv. loss scale</td><td>{0.3, 0.5, 0.7, 0.9}</td><td>0.5</td></tr><tr><td>reconstruction layer</td><td>{1, 2, 3, 4}</td><td>1</td></tr><tr><td rowspan="3">VAE-IAF<sup>21</sup></td><td>n made blocks</td><td>{2, 4, 6, 8}</td><td>4</td></tr><tr><td>n hidden in made</td><td>{2, 3, 4, 5}</td><td>4</td></tr><tr><td>hidden size</td><td>{64, 128, 256}</td><td>128</td></tr><tr><td>VAE LinNF<sup>24</sup></td><td>flows</td><td>{10P, 10R, 5P, 5R, 5P5R, 5R5P, 5PR, 5RP, 2PR, 2RP}</td><td>10R</td></tr><tr><td rowspan="2">VAMP<sup>27</sup></td><td>number components</td><td>{10, 20, 30, 40, 50}</td><td>20</td></tr><tr><td>linear scheduling steps</td><td>{0, 20, 40}</td><td>40</td></tr><tr><td rowspan="2">VQVAE<sup>28</sup></td><td>quantization loss factor</td><td>{0.25, 0.5, 0.75, 0.9, 1, 1.5, 2, 4}</td><td>2</td></tr><tr><td>n embeddings</td><td>{128, 256, 512, 1024}</td><td>512</td></tr><tr><td rowspan="3">WAE<sup>26</sup></td><td>kernel choice</td><td>{rbf, imq}</td><td>rbf</td></tr><tr><td>reg. weight</td><td>{0.01, 0.1, 0.5, 1, 5, 10, 100}</td><td>0.1</td></tr><tr><td>kernel bandwidth</td><td>{0.01, 0.1, 0.5, 1, 5, 10, 100}</td><td>5</td></tr></tbody></table>
<table><caption>**Table 4** Reconstruction metrics obtained for the best configuration of each VAE variant on the validation sets (mean ± std computed over the three splits randomly selected). Gray shading indicates models with results substantially worse than the others; these models were subsequently removed from the benchmark.</caption><thead><tr><th>Models</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th><th>PSNR ↑</th></tr></thead><tbody><tr><td>β-TC VAE<sup>16</sup></td><td>0.870 ± 0.002</td><td>1.901 ± 0.123</td><td>27.41 ± 0.27</td></tr><tr><td>β-VAE<sup>19</sup></td><td>0.868 ± 0.003</td><td>1.995 ± 0.067</td><td>27.17 ± 0.14</td></tr><tr><td>Dis. β-VAE<sup>14</sup></td><td>0.874 ± 0.006</td><td>2.004 ± 0.153</td><td>27.18 ± 0.30</td></tr><tr><td>FactorVAE<sup>20</sup></td><td>0.876 ± 0.003</td><td>1.895 ± 0.084</td><td>27.47 ± 0.14</td></tr><tr><td>HVAE<sup>15</sup></td><td>0.873 ± 0.007</td><td>1.862 ± 0.068</td><td>27.48 ± 0.14</td></tr><tr><td>InfoVAE<sup>29</sup></td><td>0.877 ± 0.006</td><td>1.813 ± 0.075</td><td>27.63 ± 0.13</td></tr><tr><td>IWAE<sup>13</sup></td><td>0.865 ± 0.007</td><td>2.087 ± 0.146</td><td>27.02 ± 0.24</td></tr><tr><td>MS-SSIM VAE<sup>25</sup></td><td>0.472 ± 0.034</td><td>70.174 ± 5.660</td><td>11.61 ± 0.36</td></tr><tr><td>RAE-GP<sup>18</sup></td><td>0.880 ± 0.006</td><td>1.715 ± 0.105</td><td>27.84 ± 0.26</td></tr><tr><td>RAE-ℓ<sup>218</sup></td><td>0.884 ± 0.005</td><td>1.815 ± 0.049</td><td>27.61 ± 0.11</td></tr><tr><td>SVAE<sup>17</sup></td><td>0.151 ± 0.001</td><td>632.694 ± 5.106</td><td>1.99 ± 0.03</td></tr><tr><td>VAEGAN<sup>22</sup></td><td>0.860 ± 0.013</td><td>2.241 ± 0.193</td><td>26.64 ± 0.38</td></tr><tr><td>VAE-IAF<sup>21</sup></td><td>0.823 ± 0.005</td><td>2.272 ± 0.057</td><td>26.65 ± 0.08</td></tr><tr><td>VAE LinNF<sup>24</sup></td><td>0.871 ± 0.001</td><td>1.855 ± 0.125</td><td>27.54 ± 0.21</td></tr><tr><td>VAMP<sup>27</sup></td><td>0.702 ± 0.097</td><td>5.581 ± 0.874</td><td>22.73 ± 0.72</td></tr><tr><td>VQVAE<sup>28</sup></td><td>0.881 ± 0.003</td><td>1.805 ± 0.032</td><td>27.62 ± 0.07</td></tr><tr><td>WAE<sup>26</sup></td><td>0.881 ± 0.005</td><td>1.862 ± 0.075</td><td>27.54 ± 0.08</td></tr></tbody></table>

contributes to the overall difficulty in training the model. For the remaining of the study, we decided not to consider the VAMP,<sup>27</sup> MS-SSIM VAE<sup>25</sup> and SVAE.<sup>17</sup>

## 4.3 Selection of the best trained models

Once the best parameters were selected through the random search, all 17 models (AE, VAE and the 15 remaining VAE-based models) were trained on the six splits of the cross-validation. We kept the same training parameters as for the random search: the models were trained on 200 epochs, with a learning rate of 10<sup>-4</sup> and a batch size of 8. There were a few exceptions: the VAE-IAF<sup>21</sup> was trained with a learning rate of 10<sup>-5</sup> to avoid errors during training. The RAE-GP<sup>18</sup> was trained with a batch size of 6, the VAEGAN<sup>22</sup> with a batch size of 4 and the IWAE<sup>13</sup> with a batch size of 2 because of the high memory usage of these models.

We then selected the best fold for each model using the average SSIM on the validation sets. The performance of all 17 models on the six splits are presented in Appendix E, with the best split of each model highlighted in bold. We can notice that the splits 2 and 3 are over-represented among the selected models. This can be explained by the fact that the cross-validation is not stratified, and the distributions of age and sex between training and validation sets for split 2 and 3 are more similar than for the other splits (Table 1).

# 5 Results obtained for the best models on the test sets

Once all the models with a correct reconstruction were trained and the best model was selected (optimal set of parameters among those tested and best split), we could evaluate each model using the procedure defined in Section 2.3. Pseudo-healthy images were reconstructed for each of the
<table><caption>**Table 5** Reconstruction metrics obtained for images from Test CN (mean ± std computed over images from the test set). The top three results for each metric are highlighted in bold.</caption><thead><tr><th>Models</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th><th>PSNR ↑</th></tr></thead><tbody><tr><td>AE</td><td>0.882 ± 0.026</td><td><strong>1.649 ± 0.613</strong></td><td><strong>28.00 ± 1.10</strong></td></tr><tr><td>Adv. AE<sup>23</sup></td><td>0.882 ± 0.028</td><td>1.707 ± 0.610</td><td>27.83 ± 1.06</td></tr><tr><td>β-TC VAE<sup>16</sup></td><td>0.878 ± 0.025</td><td>1.720 ± 0.565</td><td>27.79 ± 1.05</td></tr><tr><td>β-VAE<sup>19</sup></td><td>0.876 ± 0.027</td><td>1.846 ± 0.638</td><td>27.49 ± 1.04</td></tr><tr><td>Dis. β-VAE<sup>14</sup></td><td>0.880 ± 0.023</td><td>1.841 ± 0.634</td><td>27.50 ± 1.06</td></tr><tr><td>FactorVAE<sup>20</sup></td><td>0.879 ± 0.026</td><td>1.651 ± 0.584</td><td>27.98 ± 1.06</td></tr><tr><td>HVAE<sup>15</sup></td><td>0.882 ± 0.024</td><td>1.809 ± 0.635</td><td>27.59 ± 1.09</td></tr><tr><td>InfoVAE<sup>29</sup></td><td><strong>0.883 ± 0.024</strong></td><td>1.704 ± 0.594</td><td>27.84 ± 1.04</td></tr><tr><td>IWAE<sup>13</sup></td><td>0.876 ± 0.026</td><td>1.859 ± 0.564</td><td>27.44 ± 1.03</td></tr><tr><td>RAE-GP<sup>18</sup></td><td><strong>0.887 ± 0.023</strong></td><td><strong>1.605 ± 0.671</strong></td><td><strong>28.14 ± 1.13</strong></td></tr><tr><td>RAE-<i>l</i><sup>218</sup></td><td>0.882 ± 0.024</td><td><strong>1.631 ± 0.531</strong></td><td><strong>28.02 ± 1.02</strong></td></tr><tr><td>VAE<sup>4</sup></td><td>0.873 ± 0.028</td><td>1.736 ± 0.566</td><td>27.75 ± 1.02</td></tr><tr><td>VAEGAN<sup>22</sup></td><td>0.866 ± 0.027</td><td>2.195 ± 0.641</td><td>26.72 ± 1.04</td></tr><tr><td>VAE-IAF<sup>21</sup></td><td>0.837 ± 0.027</td><td>2.099 ± 0.720</td><td>26.92 ± 1.00</td></tr><tr><td>VAE LinNF<sup>24</sup></td><td>0.881 ± 0.023</td><td>1.807 ± 0.610</td><td>27.58 ± 1.05</td></tr><tr><td>VQVAE<sup>28</sup></td><td><strong>0.884 ± 0.026</strong></td><td><strong>1.649 ± 0.593</strong></td><td>27.99 ± 1.10</td></tr><tr><td>WAE<sup>26</sup></td><td><strong>0.883 ± 0.026</strong></td><td>1.651 ± 0.618</td><td>27.99 ± 1.10</td></tr></tbody></table>

15 test sets (the test set with the images of healthy subjects and the 14 test sets with simulated images) in order to measure the performance of the models both qualitatively by visualizing the pseudo-healthy reconstructions, and quantitatively by computing the reconstruction metrics and the healthiness score.

## 5.1 Quantitative evaluation of the pseudo-healthy reconstructions from images of control subjects

We first assessed whether the different models could correctly reconstruct images of healthy subjects from the test set by computing the SSIM, MSE and PSNR between the input and the pseudo-healthy reconstruction. Results are reported in Table 5. We observe that the reconstruction metrics of all but two models are in the same order of magnitude, with an SSIM on average between 0.873 (VAE⁴) and 0.887 (RAE-GP¹⁸), an MSE on average between 1.6 × 10⁻³ for the RAE-GP¹⁸ and 1.859 × 10⁻³ for the IWAE,¹³ and a PSNR on average between 26.7 (VAEGAN²²) and 28.1 (RAE-GP¹⁸). On the other hand, the VAEGAN²² and the VAE-IAF²¹ perform the worst, with respectively an average SSIM of 0.866 and 0.837, and an average MSE of 2.195 × 10⁻³ and 2.099 × 10⁻³, which is even worse than the vanilla VAE and the AE.

## 5.2 Quantitative evaluation of the pseudo-healthy reconstructions from images with simulated dementia

The first evaluation step with simulated data is to compute reconstruction metrics between the pseudo-healthy reconstructions obtained from these simulated data and the ground truth images used to simulate hypometabolic images, which are the targets. These results are reported in Table 6. For all the models, the reconstructions are slightly worse than for images reconstructed from the ground truth itself (Table 5) with an average SSIM between 0.854 (VAEGAN²²) and
<table><caption>**Table 6** Reconstruction metrics obtained between pseudo-healthy reconstructions obtained from the simulated images of Test AD 30 and the ground truth images (mean ± std computed over images from the test set). The top three results for each metric are highlighted in bold.</caption><thead><tr><th>Models</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th><th>PSNR ↑</th></tr></thead><tbody><tr><td>AE</td><td>0.876 ± 0.025</td><td>2.067 ± 0.643</td><td>26.99 ± 1.03</td></tr><tr><td>Adv. AE<sup>23</sup></td><td><b>0.878 ± 0.027</b></td><td><b>2.021 ± 0.620</b></td><td><b>27.07 ± 0.99</b></td></tr><tr><td>β-TC VAE<sup>16</sup></td><td>0.871 ± 0.024</td><td>2.094 ± 0.605</td><td>26.92 ± 1.02</td></tr><tr><td>β-VAE<sup>19</sup></td><td>0.873 ± 0.026</td><td>2.107 ± 0.678</td><td>26.90 ± 1.00</td></tr><tr><td>Dis. β-VAE<sup>14</sup></td><td>0.874 ± 0.023</td><td>2.158 ± 0.667</td><td>26.80 ± 1.02</td></tr><tr><td>FactorVAE<sup>20</sup></td><td>0.873 ± 0.024</td><td>2.157 ± 0.603</td><td>26.78 ± 0.94</td></tr><tr><td>HVAE<sup>15</sup></td><td>0.876 ± 0.023</td><td>2.071 ± 0.646</td><td>26.98 ± 1.04</td></tr><tr><td>InfoVAE<sup>29</sup></td><td>0.878 ± 0.022</td><td><b>2.044 ± 0.595</b></td><td><b>27.02 ± 0.97</b></td></tr><tr><td>IWAE<sup>13</sup></td><td>0.864 ± 0.027</td><td>2.265 ± 0.571</td><td>26.55 ± 0.91</td></tr><tr><td>RAE-GP<sup>18</sup></td><td><b>0.878 ± 0.022</b></td><td>2.118 ± 0.690</td><td>26.88 ± 0.99</td></tr><tr><td>RAE-<i>l</i><sup>218</sup></td><td>0.877 ± 0.023</td><td><b>1.997 ± 0.564</b></td><td><b>27.12 ± 0.99</b></td></tr><tr><td>VAE<sup>4</sup></td><td>0.870 ± 0.027</td><td>2.075 ± 0.589</td><td>26.95 ± 0.95</td></tr><tr><td>VAEGAN<sup>22</sup></td><td>0.854 ± 0.027</td><td>2.650 ± 0.662</td><td>25.88 ± 0.97</td></tr><tr><td>VAE-IAF<sup>21</sup></td><td>0.842 ± 0.025</td><td>2.322 ± 0.735</td><td>26.47 ± 0.97</td></tr><tr><td>VAE LinNF<sup>24</sup></td><td>0.876 ± 0.022</td><td>2.179 ± 0.614</td><td>26.74 ± 0.97</td></tr><tr><td>VQVAE<sup>28</sup></td><td><b>0.878 ± 0.025</b></td><td>2.089 ± 0.596</td><td>26.92 ± 0.97</td></tr><tr><td>WAE<sup>26</sup></td><td>0.877 ± 0.025</td><td>2.087 ± 0.650</td><td>26.95 ± 1.04</td></tr></tbody></table>

0.878 (Adv. AE,²³ RAE-GP¹⁸ and VQVAE²⁸), an average MSE between 1.997 × 10⁻³ for the RAE-<i>l</i><sup>218</sup> and 2.650 × 10⁻³ for the VAEGAN,²² and an average PSNR between 25.88 (VAEGAN²²) and 27.12 (RAE-<i>l</i><sup>218</sup>). The only exception is the VAE-IAF,²¹ for which the SSIM increases from 0.837 on average to 0.842. However, the reconstruction metrics are still quite high, meaning that the reconstructions from simulated hypometabolic images are similar to their target.

## 5.3 Qualitative evaluation of the pseudo-healthy reconstructions

Examples of pseudo-healthy reconstructions obtained from the original image of a control subject and images with simulated dementia are displayed in Figure 4. We first observe that all the models are able to reconstruct the input image of a healthy subject. We can recognize the shape of the brain, the areas with high metabolism (gray matter) and the others with a lower metabolism (white matter, ventricles). The VAE-IAF²¹ reconstruction has an artifact in the precuneus, which appears as a spherical hypermetabolism. This probably explains why the average SSIM is lower for the VAE-IAF²¹ than for other models. We can also see that the VAEGAN²² tends to reconstruct the image with a higher average intensity, as shown by the fact that the difference map is mostly negative (meaning that the reconstruction's voxel values are superior to the input's voxel values).

When reconstructing images with different degrees of simulated AD, we observe that all the models are able to reconstruct images that are visibly healthy by correcting the hypometabolism simulated. On the difference maps, we can recognize the mask used for the simulation as an anomaly, meaning that the model is able to reconstruct pseudo-healthy images. From this qualitative analysis, the models that seem to perform the best in terms of anomaly detection are
the VAE-IAF<sup>21</sup> (excluding the fact that it reconstructs an artifact), the β-VAE,<sup>19</sup> the disentangled β-VAE<sup>14</sup> and the HVAE,<sup>15</sup> at least for images with low (AD 15) and medium (AD 30) severity. It is indeed possible to better distinguish the abnormal area in both hemispheres on the difference maps and the reconstruction errors do not hide the anomaly.

Additional examples of pseudo-healthy reconstructions obtained for different simulated dementias are displayed in Figure 12.

## 5.4 Quantitative evaluation with the healthiness score

After qualitatively analyzing the pseudo-healthy reconstructions, we computed the healthiness score defined in Equation 2 for the different simulated test sets.

Figure 5 displays the distribution of the healthiness score for the ground truth (i.e., the images of healthy subjects), the images simulating AD with 30% hypometabolism (AD 30) and the reconstructions obtained for the different models from the AD 30 images. As expected, the healthiness of the ground truth is between 1.0 and 1.08, and it drops to between 0.83 and 0.90 when simulating AD with a hypometabolism intensity of 30%. Studying the healthiness of the pseudo-healthy reconstruction for each model, we first observe that all the models are able to reconstruct images that are healthier than the simulated input as the healthiness of the reconstructions (around 1) is superior to the healthiness of the simulated images they were reconstructed from (around 0.87). We can observe that three models seem to perform slightly better than the others, namely the β-VAE,<sup>19</sup> the disentangled β-VAE<sup>14</sup> and the VAE-IAF<sup>21</sup> with a healthiness between 0.97 and 1.04 for the first two and 0.96 and 1.03 for the third. On the other hand, the VAEGAN<sup>22</sup> appears to be the model with the worst performance (with a healthiness between 0.93 and 1.0), followed by the FactorVAE<sup>20</sup> and the RAE-GP<sup>18</sup> (which have a healthiness score between 0.94 and 1.01).

These results are consistent with the qualitative analysis, as we observed that the β-VAE,<sup>19</sup> disentangled β-VAE<sup>14</sup> and VAE-IAF<sup>21</sup> seemed to better highlight the simulated anomalies, while the VAEGAN's<sup>22</sup> poor reconstructions tended to hide the anomalies.

We also analyzed the impact of the severity of the simulated disease on the healthiness score. Figure 6 displays the evolution of the healthiness for all the models with increasing severity of simulated AD from 5% to 70%. We notice that all the models reconstruct images that are decreasingly healthy according to this metric when increasing the severity of the simulated disease. As in the previous experiment, the β-VAE<sup>19</sup> and disentangled β-VAE perform the best for high hypometabolism, followed by the VAE-IAF.<sup>21</sup> The VAEGAN<sup>22</sup> and RAE-GP<sup>18</sup> have the worst performance. However, the healthiness of the reconstruction remains above the one of simulated data, which means that all the models can reconstruct pseudo-healthy images.

Figure 7 displays for various dementia subtypes (PCA, bvFTD, lvPPA, svPPA and nfvPPA simulated at 30%) the distribution of the healthiness computed for the ground truth, the simulated images and the images reconstructed with all the models. All the models have very similar performance with a healthiness between 0.95 and 1 when simulating PCA, between 0.9 and 1.1 for bvFTD, between 0.96 and 1.6 for lvPPA, between 0.68 and 0.87 for svPPA, and between 1.0 and 1.1 for nfvPPA. As for AD, the VAEGAN's<sup>22</sup> performance is slightly lower than that of the other models, and the β-VAE<sup>19</sup> and disentangled β-VAE<sup>14</sup> seem to perform slightly better than the average. We notice that the healthiness of the ground truth, derived from the images of CN subjects, depend on the simulated dementia, given the different masks used for computation. For example, in the case of svPPA, the ground truth's healthiness tends to be lower than that of AD (falling between 0.67 and 0.92). This difference comes from the mask's location in the temporal pole for svPPA, where FDG uptake is naturally lower even in healthy images, in contrast to other regions.<sup>71</sup>
<figure><img src="image_8.png" /><figcaption>**Fig 4** Examples of reconstructions obtained with the different VAE variants from the original image of a cognitively normal subject (images of the first column, Test CN) and from the same subject with AD simulated at different intensity degrees (AD 15, AD 30, AD 50 and AD 70). The first row shows the input image in odd columns and the mask of the simulated disease in even columns when the input is a simulated image. All the other rows are the pseudo-healthy reconstructions of the models in odd columns and the difference between the pseudo-healthy reconstruction and the input in even columns. Coronal and sagittal views are displayed in Figures 10 and 11, respectively.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig 5** Ridgeline plot showing the distribution of the healthiness score for images from Test AD 30. The first row corresponds to the healthiness of the ground truth, the second row to the healthiness of the images simulating AD with 30% hypometabolism used as input, and the remaining rows to the healthiness of the pseudo-healthy reconstructions obtained with the VAE models.</figcaption></figure>

To push further the comparison of the models, we jointly analyzed their performance in terms of reconstruction accuracy and healthiness. Figure 8 displays a joint density plot of the SSIM and healthiness score computed for pseudo-healthy reconstructions obtained from images simulating AD at 30% of hypometabolism. We carefully selected four models that we compare to the VAE: the VAEGAN²² that performs worse than most models, both in terms of reconstruction accuracy and healthiness, the RAE-GP¹⁸ that has a good reconstruction but a low healthiness performance, the VAE-IAF²¹ that has the worst reconstruction accuracy but a good healthiness, and finally the β-VAE¹⁹ that has both good reconstruction and healthiness performance. This analysis confirms that, among the selected models, the β-VAE is the one that performs the best, and the VAEGAN is the one performing the worst.

**5.5 Qualitative analysis of the pseudo-healthy reconstructions obtained from real AD patients**

Even though no ground truth is available, it is important to analyze the behavior of the VAE models on data from real patients, here with AD. Figure 9 displays examples of pseudo-healthy reconstructions obtained from the image of an AD patient. This patient presents a typical hypometabolism in the parietal and temporal lobes, which is detected by all the models. However, we also observe for all the models what appears as hypermetabolism in the frontal lobe, which is
<figure><img src="image_10.png" /><figcaption>**Fig 6** healthiness score depending on the severity of the anomalies simulated. The healthiness of the ground truth, which is constant, is displayed as reference. The healthiness of the images simulating AD rapidly drops, with the hypometabolism increasing from 5% to 70%. The other curves correspond to the healthiness of the reconstructions obtained with the VAE models. Each dot represents the mean value of the healthiness, and the error bar represents the standard deviation.</figcaption></figure>

not typical of AD and probably results from reconstruction inaccuracies as this tendency was
also visible for the CN subject displayed in Figure 4, better seen in Figure 11.

# 6 Discussion

This benchmark assessed the ability of 20 VAE models to reconstruct pseudo-healthy 3D
brain FDG PET images for anomaly detection. We first searched for the best encoder-decoder
architecture for the vanilla VAE. We then optimized the hyper-parameters of all the VAE-based
models. After discarding the models with low reconstruction performance, we trained the 17
remaining ones on all the splits of the cross-validation to select the best split for each model.
Finally, we compared the trained models using conventional reconstruction metrics, as well as
the simulation framework<sup>44,56</sup> paired with the healthiness score<sup>44</sup> we previously proposed. We
found that most models tested generate images of good quality, which are healthy looking even
when given as input images with simulated anomalies. This study allowed the identification of
slightly better performing models (which can be useful for future works), whilst demonstrating
that many VAE-based models perform good anomaly detection for various anomaly intensities,
shapes and locations in 3D brain FDG PET.

## 6.1 Selection of the best configuration for each model

We performed an extensive random search to define the optimal encoder-decoder architecture for
the vanilla VAE. 200 models were trained for a total of approximately 5000 GPU hours. Since
<figure><img src="image_11.png" /><figcaption>**Fig 7** Distribution of the healthiness score depending on the dementia simulated at 30% hypometabolism: PCA, bvFTD, IvPPA, svPPA and nfVPPA. Each box plot displays the median, the lower and upper quartiles and the minimum and maximum (excluding outliers) of the healthiness. The first box (top row) shows the healthiness of the ground truth, the second one the healthiness of the simulated images used as input and the remaining ones the healthiness of the pseudo-healthy reconstructions obtained with the VAE models.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig 8** Joint density plot of the healthiness score (x-axis) and SSIM (y-axis) computed for pseudo-healthy reconstructions obtained from images from Test AD 30. A good model should appear on the top right part of the graph (high SSIM and healthiness close to 1).</figcaption></figure>

our objective was to find a small model that allows fitting 3D high resolution images in the GPU memory, the architecture we obtained is similar to what we could have implemented following examples and guidelines found in the literature: the encoder and decoder are symmetric, they are composed of five blocks, each containing only one 3D convolution layer, a batch normalization and a swish activation.⁶⁵ This architecture is for instance very similar to the one we tested in our
<figure><img src="image_13.png" /><figcaption>**Fig 9** Example of reconstructions obtained from the different VAE variants on an AD patient (on axial, coronal and sagittal slices). The first row shows the input image in odd columns. The rows below are the pseudo-healthy reconstruction of the models in odd columns, and the difference between the pseudo-healthy reconstruction and the input in even columns.</figcaption></figure>
previous work.⁴⁴ Having a small encoder and decoder proves especially advantageous in this benchmark for models with heightened memory-requirements, like the VAEGAN²² (due to its extra discriminator network), the IWAE¹³ (since it uses several samples from the latent space), and the VAE-IAF²¹ (since it has extra layers for the auto-regressive flows in the latent space). This architecture was used for all the VAE-based models. Whilst optimizing the architecture for the vanilla VAE may give this model an advantage, it was not conceivable for us, given our computational resources, to optimize the encoder-decoder architecture separately for all the models.

To optimize the hyper-parameters of each VAE variant, 324 models were trained for a total of approximately 18,000 hours of GPU use. At this stage we removed three models from the study, as they led to poor reconstructions in comparison with the others: the SVAE,¹⁷ the MS-SSIM VAE²⁵ and the VAMP.²⁷ For each remaining model, it was possible to find a set of hyper-parameters that led to good reconstruction performance.

After training the models with the selected hyper-parameters on the six splits of the cross-validation, we selected the best split for each of them. We observed that splits 2 and 3 gave the best results for 13 models out of 17 (Table 28). This might be explained by the fact that the cross-validation was not stratified, and so the validation set may not be representative of the training population for some of the splits (Table 1). This may have biased the selection of the hyper-parameters since some models were not trained on splits 2 and 3 when randomly selecting three folds out of six, which may result in underestimated performance for these configurations. However, it would have been too long to train all the configurations on all the splits; and we appraise that we still found a satisfying combination of parameters with respect to the reconstruction metrics, even though it may not be the best one.

All the selection steps were based on the validation sets, potentially leading to over-fitting on these validation sets. Performing a 6-fold cross-validation and randomly selecting the splits reduced this risk.

In this context of extensive hyper-parameter optimization and potential biases in split selection, another important observation concerns the impact of the *β* value in the *β*-VAE models. While the *β*-VAE emerged as the best-performing model in terms of both reconstruction and healthiness, the choice of *β* = 10 as the optimal value is somewhat surprising. One might expect that smaller values of *β* (e.g., *β* < 1) would favor reconstruction performance by reducing the influence of the KL divergence term. However, as shown in Table 12 in Appendix D.3, *β* = 10 yielded the best validation performance, appearing to be an exception rather than a consistent trend. This suggests that model performance is highly sensitive to the choice of *β*, and potentially to interactions with other hyper-parameters or training splits. A more thorough investigation of this hyper-parameter, using a broader range of values and possibly visualizing the results with joint plots (as in Figure 8), could help disentangle these effects and better understand how *β* influences the balance between reconstruction quality and latent space regularization.

## 6.2 Comparison of the VAE variants

To evaluate the different VAE variants, we applied the evaluation procedure presented in Hassanaly et al.⁴⁴ This evaluation consists in two main steps: first measuring the reconstruction performance of the models using reconstruction metrics for images of healthy subjects, and then using simulated data in order to evaluate the ability of the models to reconstruct pseudo-healthy images (i.e. whether the reconstructions appear healthy).

In terms of reconstruction metrics, for the images of healthy subjects (Table 5), all the trained models led to similar performance, with the nine best models having an average SSIM above 0.88, seven models having an average SSIM between 0.86 and 0.88, and only one having an
average SSIM below 0.86. All the models were able to reconstruct realistic brain images, as shown in Figure 4. Many models performed better than the vanilla VAE on average according to both the MSE and the SSIM, but not substantially. The reconstruction metrics computed between the reconstructions obtained from simulated images and the ground truth (Table 6) show that the reconstructed images are quite similar to their healthy target (i.e., the original images used to simulate hypometabolic images), indicating that the reconstruction capacity of the models is not affected when using images with anomalies as input.

In terms of healthiness, on images simulating AD, it appears that the $\beta$-VAE,¹⁹ the disentangled $\beta$-VAE¹⁴ and the VAE-IAF²¹ performed better than the other VAEs, whereas the VAEGAN²² and the RAE-GP¹⁸ gave the worst results (Figure 5). Interestingly, the healthiness distribution of the ground truth images is multi-modal, and so is logically the distribution of the healthiness computed for the simulated images. This is also the case of the healthiness distributions of the reconstruction for most of the models, especially for the FactorVAE²⁰ and the VQVAE²⁸ for which we can properly recognize the shape of the distribution. However, the reconstructions of the two best performing models, the $\beta$-VAE¹⁹ and the disentangled $\beta$-VAE,¹⁴ have a uni-modal healthiness distribution, potentially explaining why they are not the best performing models in terms of reconstruction metrics. This may be explained by the fact that, in both cases, we set $\beta = 10 (\gg 1)$, giving more weight to the KL term than the reconstruction term in the loss.

As illustrated in Figure 7, most of the models were able to reconstruct images of healthy appearance also for dementia subtypes different from AD. As for AD, the best performing models on average were the $\beta$-VAE¹⁹ and the disentangled $\beta$-VAE.¹⁴ The VAE LinNF²⁴ and the HVAE¹⁵ also seem to have a higher healthiness than the other models. On the other hand, the VAEGAN,²² the RAE-GP¹⁸ and the FactorVAE²⁰ were the models with the lowest performance. We could further see that the healthiness score was not optimal for all the dementia subtypes. For instance, for PCA and svPPA, the healthiness of simulated images was not substantially different from that of the ground truth. Nevertheless, we observed that the healthiness of the reconstruction was close to that of the ground truth and higher than that of simulated data, which is sufficient to assess the healthiness of reconstructed images.

In general, we observed that all the models were able to reconstruct images with a healthiness substantially above the healthiness of simulated images, regardless of the kind of simulated anomalies, and almost equal to the healthiness of the ground truth, indicating that the reconstructions are indeed healthy looking.

To push further the model comparison, we jointly analyzed reconstruction and healthiness score (Figure 8). The RAE-GP¹⁸ was the model with the best reconstruction, but was ranked among the worst in terms of healthiness. Even though the reconstructions look healthy when compared to the simulated input, it means that the RAE-GP¹⁸ did not learn the healthy image distribution as well as other models, but rather learned to reconstruct the input as is. On the contrary, the VAE-IAF²¹ was the model with the worst reconstruction, but was among the best in terms of healthiness. This can be explained by the presence of a reconstruction artifact that impacts the reconstruction score. The $\beta$-VAE¹⁹ was the best model in terms of healthiness and was average in terms of reconstruction, and the VAEGAN²² under-performed both in terms of reconstruction and healthiness.

A surprising result highlighted by the benchmark is that the simple AE performed on a par with more complex models, especially according to the healthiness for simulated data. Although it was expected that this model would be able to reconstruct images of healthy subjects, there was no certainty that the AE would be able to reconstruct healthy looking images from simulated images, especially when simulating severe hypometabolism (50% and more). It would be interesting to assess the performance of this model when given real images from AD patients.
In our preliminary benchmark study,⁴⁵ we compared a subset of the models that we present here on down-sampled 3D brain FDG PET. Another major difference with the present work is that we had trained the models with default hyper-parameters' values. We observe that some of the models that performed poorly in this previous study, such as the VAEGAN²² and the VAE LinNF,²⁴ performed much better after searching for optimal hyper-parameters, whereas the VAMP²⁷ and the MS-SSIM VAE²⁵ still perform poorly even after hyper-parameter tuning. Even though not surprising, this highlights the benefit and need of optimizing each model, even though this step does not guarantee reaching good performance.

## 6.3 Limitations and perspectives

Despite the absence of ground truth masks for the anomalies we aim to detect, this benchmark demonstrated the effectiveness of the simulation-based framework we previously proposed for evaluating anomaly detection methods based on pseudo-healthy reconstruction. Specifically, the framework enables model evaluation by leveraging pairs of abnormal and healthy images from the same subjects–spanning various dementia subtypes and severity levels–to assess the quality of the reconstructed pseudo-healthy images. The evaluation framework also introduced the healthiness score that automatically quantifies whether a reconstruction is pseudo-healthy. This framework serves as a valuable first evaluation step that does not require the involvement of a clinician; an important advantage, given that it would be impractical to ask clinicians to assess reconstructions from 20 different models.

However, this also highlights the main limitation of our work: the absence of evaluation of the models' capacity to detect real anomalies using these pseudo-healthy reconstructions. Without ground truth annotations or expert evaluation, it is impossible to compute metrics such as Dice score or average precision, as done in the literature.³ Therefore, our models are only compared on their ability to reconstruct pseudo-healthy images, which is insufficient to assess whether they enable the detection of real anomalies. As discussed in Section 6.4, reconstruction-based methods alone are insufficient. A solution has been proposed by Hassanaly et al.;⁴⁴ however, it is applicable only to group-level studies, thereby eliminating the capacity for subject-specific anomaly detection evaluation. A comprehensive evaluation would ultimately require performing a reader study, where the results obtained from real images with real anomalies would be reviewed by clinicians.

Beyond this, another important limitation of the current evaluation is that it focuses solely on the quality of the reconstructions and their degree of healthiness, without directly assessing how well each model learned the healthy distribution. An interesting direction for future work would be to compare the latent distributions of the trained models to assess whether the posterior learned by the different models is the same for images from healthy and diseased subjects. This could be achieved using the simulation framework⁴⁴ and comparing the latent representations of both the original and simulated images. Such an analysis could provide insights into the performance differences observed across the various VAE models and potentially suggest avenues for improvement.

The models were compared on a single modality, FDG PET. It would be further interesting to test these models on structural MRI, which have different characteristics, such as sharp structures. This would also allow us to compare the performance of these VAE variants with other approaches in the literature, as many have been developed to detect lesions in structural MRI. Similarly, it would be interesting to include other VAE models that performed well in computer vision such as Hierarchical VAEs,⁶⁵ or compare VAEs to other generative models such as GANs and diffusion models. Finally, we could combine multiple modalities and use multi-modal models.⁷²
## 6.4 Perspectives on the use of VAEs to reconstruct pseudo-healthy images for anomaly detection

From a broader perspective, this discussion addresses the limitations of using VAE pseudo-healthy reconstructions for anomaly detection. We demonstrated that VAE-based models are effective for generating pseudo-healthy images. These models are easy to implement and train, and they produce subject-specific pseudo-healthy reconstructions by learning the distribution of healthy data.

However, detecting anomalies from pseudo-healthy reconstructions alone is not straightforward and often requires several additional steps. The first step is building anomaly maps. This can be done by computing the residual, i.e. the voxel-wise difference between the patient's image and the model's pseudo-healthy reconstruction, or z-scores for greater accuracy.⁴²,⁷¹ Then, to build an anomaly score, the anomaly maps must be further processed, which may involve determining appropriate thresholds or defining suitable scores.¹⁰,⁷³ This means that relying solely on reconstruction-based metrics does not allow for a full assessment of a method's anomaly detection capabilities.⁷⁴ This issue is particularly crucial for VAE-based reconstruction, as VAEs are known to produce blurry images.

Another limitation of VAEs is their tendency to overgeneralize.⁷⁵,⁷⁶ VAEs are powerful generative models that can generalize well, sometimes reconstructing both normal and anomalous data with similarly low errors. This makes it difficult to distinguish anomalies based solely on reconstruction, as anomalies do not necessarily produce high reconstruction errors.⁷⁵ This is particularly challenging for anomaly detection in brain FDG PET, where individual metabolism variations can closely resemble subtle hypometabolism caused by dementia. For instance, if the training dataset includes numerous healthy subjects with hypometabolism due to aging, the model would learn to reconstruct age-related hypometabolism, potentially leading to accurate reconstruction of dementia-related hypometabolism as well. This means that this approach may not be appropriate when the anomalies to be detected closely resemble the inherent diversity of the dataset.

## 6.5 Reproducibility

In order to make this study as reproducible as possible,⁷⁷ we tried to follow the guidelines of the MICCAI reproducibility checklist*:

* the publicly available dataset and final cohort we work with is mainly described in Section 3 with details of the preprocessing and data selection steps presented in Appendix A. We provide a summary of participant demographics for the train, validation and test splits in Table 1;

* the architecture choices for the VAE and the impact of those choices are detailed in Section 4.1 and Appendix B;

* the VAE variants are described in Appendix D with the range of hyper-parameters considered for each of them;

* the training protocol and the method to tune and select hyper-parameters are described in Section 4.2 and Appendix D;

* we also provide a clear definition of the specific evaluation metrics and statistics used to report results in Section 2.3.

Moreover, most of the code that we used is available in ClinicaDL,⁶² an open-source software that is developed to enable reproducible deep learning studies in neuroimaging. Pipelines are available to perform the following steps:

<footnote>*https://miccai2021.org/files/downloads/MICCAI2021-Reproducibility-Checklist.pdf</footnote>
* selecting subjects from a neuroimaging dataset,

* rigorously separating data into independent training and testing sets,

* rigorously splitting the training set using a cross-validation,

* launching random searches to optimize architecture and hyper-parameters,

* easily training VAE-based models on neuroimages,

* constructing new test sets by generating simulated data using the method described in Hassanaly et al.,⁴⁴

* reconstructing pseudo-healthy images from trained models for the test sets and computing the reconstruction metrics used in evaluation.

All the VAE-based models are implemented in Pythae,⁴³ an open-source Python library that aims at unifying the implementation of VAE-based models, and facilitating benchmarks. Moreover, all the preprocessing pipelines are available in Clinica,⁵⁹ an open-source software for reproducible processing of neuroimaging datasets. In the present work, Clinica has been used to:

* curate and organize the ADNI dataset following a community standard, namely the brain imaging data structure (BIDS),⁷⁸

* perform linear registration and intensity normalization of the FDG PET scans (pet-linear pipeline).

Finally, all the code for random searches, model training and evaluation is available in the following repository: github.com/ravih18/UAD_VAE_benchmark. This repository includes dependencies and software versions used.

# 7 Conclusion

In summary, we presented a benchmark of twenty VAE-based models for the unsupervised detection of dementia related anomalies in 3D brain FDG PET. The aim was to introduce the use of recent VAE variants with medical imaging data of high dimension and compare their performance. We proposed a random search method to find the optimal architecture for the vanilla VAE, as well as a random search method to tune the hyper-parameters of the implemented models.

We observed that 17 of the 20 models had a good reconstruction quality. Using our previously proposed evaluation framework,⁴⁴ we showed that the 17 models were able to reconstruct pseudo-healthy images when fed with simulated abnormal images. By simulating AD with varying intensity and dementia other than AD, we also showed that these models were able to generalize to anomalies of different shapes, localizations and intensities. If no model clearly outperformed the others, the β-VAE¹⁹ and disentangled β-VAE¹⁴ slightly outperformed the other models, while remaining easy to tune and not being noticeably computationally costly.

Even if it is recognized that VAEs generate blurry images, all these experiments showed that most of the models were able to reconstruct good quality pseudo-healthy 3D FDG PET. The VAE variants showed similar performance and did not systematically outperform the vanilla VAE (or even the simple AE).

Finally, we can conclude that most VAEs are well suited for pseudo-healthy reconstruction of brain FDG PET images.
## A Data selection

The ADNI database comprises at least 3511 FDG PET scans from 1600 participants. This includes 1010 images of 554 cognitively normal (CN) subjects. To ensure that images really correspond to a healthy brain, as we know that physiological changes can appear several years before the first clinical symptoms, we kept only scans from subjects that are CN for at least three years after the session considered. We discarded 78 subjects (129 images) for whom diagnosis progresses to AD, 72 subjects (72 images) for whom there is a unique session (which is not enough to assess the reliability of the CN label) and 21 subjects (49 images) for whom there are multiple conversions or regressions. We finally kept 383 stable CN subjects (760 images).

The ADNI database also includes 560 AD patients (791 images). We removed 2 patients (2 images) with unstable AD diagnosis, 3 patients (3 images) of regressive AD, 189 patients (189 images) for whom there is a unique session, and 4 subjects that were already in the training set. We kept the 362 baseline sessions of the remaining AD patients for testing purposes and discarded all the other images.

In addition, we ran quality control and discarded in total 30 images:⁶² 18 images from CN subjects and 9 images from AD patients after t1-linear quality control, and 3 images from CN subjects after pet-linear quality control. In the end, we have 378 CN subjects (739 images) and 353 AD patients (353 images).

## B Details of the encoder-decoder architecture selection procedure

In previous works,⁴⁴,⁵⁶ we set the latent space size to 128 as a trade-off between performance and resources but observed that a larger latent space would lead to better reconstructions. We therefore decided to try sizes from the set {256, 512, 1024}. For choosing the number of blocks for the encoder and decoder, $B_e$ and $B_d$, we initially tried the integer range from 3 to 7. This parameter influences the size of the last feature map before the fully connected layer, and therefore the number of parameters in that layer. We noticed that having an encoder with 3 blocks leads to a very large number of parameters in the fully connected layer (around 750,000 if $C_e = 16$ or double if we double $C_e$), whereas an encoder with 7 blocks would lead to very small feature maps (1 × 1 × 1). We therefore reduced the range, and chose values between 4 and 6. We kept the number of sub-blocks in the encoder and decoder $S_e$ and $S_d$ relatively small to restrain the number of parameters of our model whilst still testing deep architectures. We chose the number of channels $C_e$ and $C_d$ based on previous experiments and decided to set it to either 16 or 32. We also added the possibility to add a convolution layer in our last decoder block (shown by dotted lines in Figure 2). We also included the learning rate and the optimizer as parameters in our random search. We first performed experiments where the learning rate was chosen from {10⁻⁵, 10⁻⁴, 10⁻³}, but setting it to 10⁻³ led to errors in the computation of the loss, so we pursued our search with only 10⁻⁵ and 10⁻⁴ as options. The optimizer could be either Adam or Adamax, following the suggestions from.⁶⁵ The parameters included in our random search are summarized in Table 2.

After attempting to train 200 models, a pattern emerged and we could select and test an additional architecture following our intuition. These results are summarized in Table 7 and Table 8. Certain parameters, such as the learning rate, the latent space size, and the number of channels $C_e$ and $C_d$ were easy to choose as a clear relation with the reconstruction metrics could be observed, allowing us to choose to set the learning rate to 10⁻⁴, the latent space size to 256 and $D_e = C_d = 16$. We particularly struggled to train models with residual sub-blocks due to memory constraints, and those that were able to train did not give very good results, so we focused our efforts on models with convolutional sub-blocks. The optimizer did not seem to
make a large difference, so we chose Adam. For the remaining parameters, we observed that there was no need for a very deep architecture (or large number of sub-blocks within blocks), whereas a large number of blocks was beneficial in terms of memory as it induced a smaller number of parameters before the fully connected layer. After analyzing these results and noticing these patterns, we decided to test an extra model which we designed to be symmetrical (for the sake of simplicity) and as light as possible in terms of memory (128 MB instead of 286 MB), as we knew that some of the models that we would train later with this architecture are much more memory greedy. We selected random splits by drawing 3 cards from a deck of 6 cards to train our model and found that this model performed similarly to the best performing models from our random search (equivalent SSIM and best MSE). We therefore decided to select this architecture as it was simpler (because symmetrical) and smaller in terms of memory and number of parameters.
<table><caption>**Table 7** Results of the random search on the VAE architecture: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected). The selected configuration is highlighted in gray.</caption><thead><tr><th>Id</th><th>block type</th><th>C<sub>e</sub></th><th>B<sub>e</sub></th><th>S<sub>e</sub></th><th>latent size</th><th>B<sub>d</sub></th><th>S<sub>d</sub></th><th>C<sub>d</sub></th><th>last conv</th><th>learning rate</th><th>SSIM ↑</th><th>MSE ↓</th></tr></thead><tbody><tr><td>1</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>256</td><td>5</td><td>2</td><td>16</td><td>True</td><td>0.0001</td><td>0.866 ± 0.006</td><td>2.158 ± 0.043</td></tr><tr><td>2</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>256</td><td>4</td><td>3</td><td>32</td><td>False</td><td>0.0001</td><td>0.864 ± 0.008</td><td>2.221 ± 0.111</td></tr><tr><td>3</td><td>conv</td><td>16</td><td>5</td><td>3</td><td>1024</td><td>4</td><td>3</td><td>16</td><td>False</td><td>0.0001</td><td>0.861 ± 0.017</td><td>2.202 ± 0.082</td></tr><tr><td>4</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>256</td><td>5</td><td>1</td><td>32</td><td>False</td><td>0.00001</td><td>0.861 ± 0.004</td><td>2.081 ± 0.021</td></tr><tr><td>5</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>512</td><td>5</td><td>2</td><td>32</td><td>False</td><td>0.0001</td><td>0.857 ± 0.001</td><td>2.153 ± 0.056</td></tr><tr><td>6</td><td>conv</td><td>16</td><td>5</td><td>1</td><td>256</td><td>5</td><td>1</td><td>16</td><td>False</td><td>0.0001</td><td>0.856 ± 0.004</td><td>1.919 ± 0.088</td></tr><tr><td>7</td><td>conv</td><td>32</td><td>5</td><td>3</td><td>512</td><td>5</td><td>1</td><td>16</td><td>True</td><td>0.0001</td><td>0.852 ± 0.002</td><td>2.138 ± 0.101</td></tr><tr><td>8</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>512</td><td>4</td><td>1</td><td>32</td><td>True</td><td>0.0001</td><td>0.851 ± 0.006</td><td>2.303 ± 0.187</td></tr><tr><td>9</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>1024</td><td>4</td><td>1</td><td>16</td><td>True</td><td>0.0001</td><td>0.850 ± 0.007</td><td>2.572 ± 0.145</td></tr><tr><td>10</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>1024</td><td>5</td><td>1</td><td>16</td><td>True</td><td>0.0001</td><td>0.840 ± 0.024</td><td>2.014 ± 0.083</td></tr></tbody></table>

<table><caption>**Table 8** Results of the random search on the VAE architecture: ranking according to the MSE of the 10 best configurations (mean ± std over the three folds randomly selected). The selected configuration is highlighted in gray. The Id corresponds to the rank of the model when sorting them according to the SSIM (Table 7).</caption><thead><tr><th>Id</th><th>block type</th><th>C<sub>e</sub></th><th>B<sub>e</sub></th><th>S<sub>e</sub></th><th>latent size</th><th>B<sub>d</sub></th><th>S<sub>d</sub></th><th>C<sub>d</sub></th><th>last conv</th><th>learning rate</th><th>SSIM ↑</th><th>MSE ↓</th></tr></thead><tbody><tr><td>6</td><td>conv</td><td>16</td><td>5</td><td>1</td><td>256</td><td>5</td><td>1</td><td>16</td><td>False</td><td>0.0001</td><td>0.856 ± 0.004</td><td>1.919 ± 0.088</td></tr><tr><td>10</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>1024</td><td>5</td><td>1</td><td>16</td><td>True</td><td>0.0001</td><td>0.840 ± 0.024</td><td>2.014 ± 0.083</td></tr><tr><td>14</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>1024</td><td>4</td><td>3</td><td>16</td><td>False</td><td>0.00001</td><td>0.826 ± 0.013</td><td>2.062 ± 0.065</td></tr><tr><td>4</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>256</td><td>5</td><td>1</td><td>32</td><td>False</td><td>0.00001</td><td>0.861 ± 0.004</td><td>2.081 ± 0.021</td></tr><tr><td>11</td><td>conv</td><td>16</td><td>5</td><td>1</td><td>256</td><td>4</td><td>1</td><td>32</td><td>True</td><td>0.0001</td><td>0.800 ± 0.049</td><td>2.081 ± 0.107</td></tr><tr><td>13</td><td>conv</td><td>32</td><td>6</td><td>1</td><td>256</td><td>6</td><td>2</td><td>32</td><td>False</td><td>0.00001</td><td>0.832 ± 0.009</td><td>2.091 ± 0.113</td></tr><tr><td>12</td><td>conv</td><td>16</td><td>5</td><td>3</td><td>512</td><td>4</td><td>2</td><td>16</td><td>True</td><td>0.0001</td><td>0.818 ± 0.032</td><td>2.111 ± 0.069</td></tr><tr><td>7</td><td>conv</td><td>32</td><td>5</td><td>3</td><td>512</td><td>5</td><td>1</td><td>16</td><td>True</td><td>0.0001</td><td>0.853 ± 0.002</td><td>2.128 ± 0.101</td></tr><tr><td>5</td><td>conv</td><td>32</td><td>5</td><td>1</td><td>512</td><td>5</td><td>2</td><td>32</td><td>False</td><td>0.0001</td><td>0.857 ± 0.001</td><td>2.153 ± 0.056</td></tr><tr><td>1</td><td>conv</td><td>16</td><td>5</td><td>2</td><td>256</td><td>5</td><td>2</td><td>16</td><td>True</td><td>0.0001</td><td>0.866 ± 0.006</td><td>2.158 ± 0.043</td></tr></tbody></table>
C Detailed description of the final encoder-decoder architecture

<table><caption>**Table 9** Summary of the layer parameters in the final VAE architecture</caption><thead><tr><th>Layer</th><th>Input shape</th><th>Output shape</th><th>Kernel</th><th>Stride</th><th>Padding</th><th>BN</th><th>Activation</th></tr></thead><tbody><tr><td rowspan="7">Encoder</td><td>Conv3D</td><td>(1, 169, 208, 179)</td><td>(16, 84, 104, 89)</td><td>(4, 4, 4)</td><td>(2, 2, 2)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Conv3D</td><td>(16, 84, 104, 89)</td><td>(32, 42, 52, 44)</td><td>(4, 4, 4)</td><td>(2, 2, 2)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Conv3D</td><td>(32, 42, 52, 44)</td><td>(64, 21, 26, 22)</td><td>(4, 4, 4)</td><td>(2, 2, 2)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Conv3D</td><td>(64, 21, 26, 22)</td><td>(128, 10, 13, 11)</td><td>(4, 4, 4)</td><td>(2, 2, 2)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Conv3D</td><td>(128, 10, 13, 11)</td><td>(256, 5, 6, 5)</td><td>(4, 4, 4)</td><td>(2, 2, 2)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Flatten</td><td>(256, 5, 6, 5)</td><td>(38400)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>FC × 2</td><td>(38400)</td><td>(256) × 2</td><td>-</td><td>-</td><td>-</td><td>False</td><td>SiLU</td></tr><tr><td rowspan="13">Decoder</td><td>FC</td><td>(256)</td><td>(38400)</td><td>-</td><td>-</td><td>-</td><td>False</td><td>ReLU</td></tr><tr><td>Unflatten</td><td>(38400)</td><td>(256, 5, 6, 5)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Upsample</td><td>(256, 5, 6, 5)</td><td>(256, 10, 13, 11)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Conv3D</td><td>(256, 10, 13, 11)</td><td>(128, 10, 13, 11)</td><td>(3, 3, 3)</td><td>(1, 1, 1)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Upsample</td><td>(128, 10, 13, 11)</td><td>(128, 21, 26, 22)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Conv3D</td><td>(128, 21, 26, 22)</td><td>(64, 21, 26, 22)</td><td>(3, 3, 3)</td><td>(1, 1, 1)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Upsample</td><td>(64, 21, 26, 22)</td><td>(64, 42, 52, 44)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Conv3D</td><td>(64, 42, 52, 44)</td><td>(32, 42, 52, 44)</td><td>(3, 3, 3)</td><td>(1, 1, 1)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Upsample</td><td>(32, 42, 52, 44)</td><td>(32, 84, 104, 89)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Conv3D</td><td>(32, 84, 104, 89)</td><td>(16, 84, 104, 89)</td><td>(3, 3, 3)</td><td>(1, 1, 1)</td><td>(1, 1, 1)</td><td>True</td><td>SiLU</td></tr><tr><td>Upsample</td><td>(16, 84, 104, 89)</td><td>(16, 169, 208, 179)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Conv3D</td><td>(16, 169, 208, 179)</td><td>(1, 169, 208, 179)</td><td>(3, 3, 3)</td><td>(1, 1, 1)</td><td>(1, 1, 1)</td><td>False</td><td>Sigmoid</td></tr></tbody></table>

BN: batch normalization, FC: fully connected, SiLU: swish activation function
# D Description of the VAE variants and of their hyper-parameter selection procedure

This section describes all the VAE variants. Hyper-parameters were chosen following implementations and recommendations from the original papers and the benchmark previously done by Chadebec et al.⁴³ The results of the random searches are reported for each of the models. For each model, we selected the configuration with the best average SSIM on the validation folds.

## D.1 Adversarial Autoencoder

The adversarial autoencoder²³ is a probabilistic autoencoder model that uses the GAN framework to perform variational inference in the latent space. It uses a discriminator network to differentiate a prior's sample from a posterior's sample as a form of regularization. Its objective and training are quite similar to that of a VAE

$$ \mathcal{L}_{\text{Adv. AE}} = \mathbb{E}_{z \sim q_{\phi}(z|x)} [\log p_{\theta}(x|z)] + \alpha \mathcal{L}_{\text{GAN}}, $$

where

$$ \mathcal{L}_{\text{GAN}} = \mathbb{E}_{\tilde{z} \sim p_z(z)} [\log(1 - D(\tilde{z}))] + \mathbb{E}_{x \sim p_{\theta}} [\mathbb{E}_{z \sim q_{\phi}(z|x)} [\log D(z)]]. $$

We set the discriminator to be the same as in Chadebec et al.,⁴³ that is, a multilayer perceptron with a single hidden layer with 256 units and ReLU activation. We performed a grid search of 10 configurations for

$$ \alpha \in \{0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99\}. $$

The results are reported in Table 10.

<table><caption>**Table 10** Results of the random search on the hyper-parameters of the Adv. AE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>adversarial loss scale</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>0.9</td><td>0.873 ± 0.005</td><td>1.770 ± 0.083</td></tr><tr><td>0.01</td><td>0.872 ± 0.003</td><td>1.771 ± 0.144</td></tr><tr><td>0.1</td><td>0.869 ± 0.005</td><td>1.846 ± 0.115</td></tr><tr><td>0.5</td><td>0.869 ± 0.015</td><td>1.811 ± 0.094</td></tr><tr><td>0.75</td><td>0.869 ± 0.006</td><td>1.784 ± 0.142</td></tr><tr><td>0.25</td><td>0.866 ± 0.002</td><td>1.841 ± 0.155</td></tr><tr><td>0.99</td><td>0.865 ± 0.014</td><td>1.863 ± 0.094</td></tr><tr><td>0.05</td><td>0.863 ± 0.009</td><td>1.779 ± 0.103</td></tr><tr><td>0.001</td><td>0.863 ± 0.007</td><td>1.860 ± 0.075</td></tr><tr><td>0.95</td><td>0.856 ± 0.001</td><td>1.814 ± 0.118</td></tr></tbody></table>

## D.2 β-TC VAE

The β-TC VAE, or Total Correlation VAE,¹⁶ is an extension of the β-VAE,¹⁹ which aims at further isolating sources of disentanglement by rewriting the ELBO in the following way:

$$ \mathcal{L}_{\text{β-TC VAE}} = \mathbb{E}_{z \sim q_{\phi}(z|x)} [\log p_{\theta}(x|z)] - \mathcal{L}_{\text{reg}}, $$
where

$$ \mathcal{L}_{\text{reg}} = \alpha \mathcal{D}_{\text{KL}} [q_{\phi}(z, x) || q_{\phi}(z) p_{\theta}(x)] + \beta \mathcal{D}_{\text{KL}} \left[ q_{\phi}(z) || \prod_{j} q_{\phi}(z_j) \right] + \gamma \sum_{j} \mathcal{D}_{\text{KL}} [q_{\phi}(z_j) || p_{z}(z_j)] . $$

The regularization term is therefore the sum of the mutual information between x and z, the total correlation, which models the dependence between dimensions of the latent vector, and the dimension-wise KL divergence, which prevents each dimension of the latent variable from diverging too far from its prior.

Following the authors' suggestion, we set $\alpha = \gamma = 1$ for most of the models and performed a grid search for parameter $\beta \in \{0.001, 0.005, 0.01, 0.05, 0.1, 1, 2, 5, 10\}$. We also tried the configurations $(\beta, \alpha, \gamma) = (1, 1, 3)$ and $(\beta, \alpha, \gamma) = (1, 3, 1)$, which made a total of 12 configurations. The results are reported in Table 11.

<table><caption>**Table 11** Results of the random search on the hyper-parameters of the β-TC VAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>β</th><th>α</th><th>γ</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>2</td><td>1</td><td>1</td><td>0.870 ± 0.002</td><td>1.901 ± 0.123</td></tr><tr><td>0.05</td><td>1</td><td>1</td><td>0.866 ± 0.002</td><td>1.953 ± 0.082</td></tr><tr><td>1</td><td>1</td><td>3</td><td>0.866 ± 0.003</td><td>1.993 ± 0.077</td></tr><tr><td>5</td><td>1</td><td>1</td><td>0.864 ± 0.004</td><td>1.923 ± 0.072</td></tr><tr><td>0.005</td><td>1</td><td>1</td><td>0.864 ± 0.009</td><td>1.871 ± 0.113</td></tr><tr><td>0.001</td><td>1</td><td>1</td><td>0.863 ± 0.005</td><td>1.903 ± 0.138</td></tr><tr><td>1</td><td>3</td><td>1</td><td>0.862 ± 0.010</td><td>1.810 ± 0.034</td></tr><tr><td>10</td><td>1</td><td>1</td><td>0.862 ± 0.008</td><td>1.969 ± 0.095</td></tr><tr><td>0.01</td><td>1</td><td>1</td><td>0.860 ± 0.010</td><td>1.917 ± 0.096</td></tr><tr><td>0.1</td><td>1</td><td>1</td><td>0.855 ± 0.010</td><td>1.864 ± 0.100</td></tr></tbody></table>

## D.3 β-VAE

The β-VAE¹⁹ was introduced to encourage the disentanglement of features in the latent space by adding a weight $\beta$ in front of the KL term to adjust the balance between reconstruction and regularization. The objective is:

$$ \mathcal{L}_{\text{β-VAE}} = \mathbb{E}_{z \sim q_{\phi}(z|x)} [\log p_{\theta}(x|z)] - \beta \mathcal{D}_{\text{KL}} [q_{\phi}(z|x) || p_{z}(z)] , $$

where setting $\beta > 1$ leads to stronger disentanglement whereas using a smaller $\beta$ can favor better reconstruction abilities.

We performed a grid search of 10 configurations for $\beta \in \{0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 2, 5, 10, 100\}$. The results are reported in Table 12.
<table><caption>**Table 12** Results of the random search on the hyper-parameters of the β-VAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>β</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>10</td><td>0.868 ± 0.003</td><td>1.995 ± 0.067</td></tr><tr><td>0.005</td><td>0.868 ± 0.006</td><td>1.785 ± 0.142</td></tr><tr><td>0.01</td><td>0.867 ± 0.006</td><td>1.755 ± 0.122</td></tr><tr><td>0.001</td><td>0.866 ± 0.005</td><td>1.825 ± 0.072</td></tr><tr><td>0.05</td><td>0.866 ± 0.008</td><td>1.859 ± 0.090</td></tr><tr><td>2</td><td>0.863 ± 0.009</td><td>1.9 ± 0.072</td></tr><tr><td>0.1</td><td>0.863 ± 0.011</td><td>1.894 ± 0.103</td></tr><tr><td>0.5</td><td>0.858 ± 0.011</td><td>1.835 ± 0.117</td></tr><tr><td>5</td><td>0.856 ± 0.011</td><td>1.969 ± 0.095</td></tr><tr><td>100</td><td>0.816 ± 0.008</td><td>3.716 ± 0.292</td></tr></tbody></table>

## D.4 Disentangled β-VAE

The disentangled β-VAE¹⁴ introduces a way to progressively increase the latent encoding capacity to improve the reconstruction accuracy in comparison with the β-VAE.¹⁹ The objective becomes

$$ \mathcal{L}_{\text{disentangled β-VAE}} = \mathcal{L}_{\text{rec}} - \beta |\mathcal{D}_{KL}(q_{\phi}(z|x)||p(z)) - C|, $$

with C the value of the KL divergence term we would like to approach.

We performed a random search on the three parameters: $\beta \in \{10^{-2}, 10^{-1}, 1, 5, 10\}$, $C \in \{5, 25, 50\}$ and the number of epochs (warm-up epochs) during which the KL divergence in the ELBO will increase from 0 to C, which can be 100 or 1000. We trained a total of 20 configurations (out of 60 possible combinations), and the results of the random search are given in the Table 13.

<table><caption>**Table 13** Results of the random search on the hyper-parameters of the Dis. β-VAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>β</th><th>C</th><th>warm-up epoch</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>10</td><td>50</td><td>1000</td><td>0.874 ± 0.006</td><td>2.004 ± 0.153</td></tr><tr><td>0.1</td><td>25</td><td>1000</td><td>0.873 ± 0.002</td><td>1.821 ± 0.056</td></tr><tr><td>1</td><td>50</td><td>100</td><td>0.873 ± 0.007</td><td>1.852 ± 0.092</td></tr><tr><td>0.01</td><td>5</td><td>1000</td><td>0.871 ± 0.004</td><td>1.755 ± 0.055</td></tr><tr><td>0.1</td><td>50</td><td>100</td><td>0.871 ± 0.009</td><td>1.869 ± 0.073</td></tr><tr><td>0.01</td><td>25</td><td>1000</td><td>0.870 ± 0.003</td><td>1.753 ± 0.074</td></tr><tr><td>10</td><td>5</td><td>1000</td><td>0.870 ± 0.003</td><td>2.053 ± 0.110</td></tr><tr><td>0.1</td><td>5</td><td>1000</td><td>0.869 ± 0.014</td><td>1.815 ± 0.036</td></tr><tr><td>1</td><td>5</td><td>100</td><td>0.869 ± 0.008</td><td>1.879 ± 0.064</td></tr><tr><td>5</td><td>25</td><td>1000</td><td>0.867 ± 0.002</td><td>2.009 ± 0.068</td></tr></tbody></table>

## D.5 Factor VAE

Kim et al²⁰ proposed a new metric for disentanglement that encourages the latent representation to be factorial, and independent across each dimension of the latent space. The loss function is
the following:

$$ \mathcal{L}_{\text{FactorVAE}} = \mathcal{L}_{\text{VAE}} - \gamma \mathcal{D}_{KL}(q_{\phi}(z) || \bar{q}_{\phi}(z))) , $$

with $\bar{q}_{\phi}(z) := \prod_{j=1}^{d} q_{\phi}(z_j)$ for a model with a latent space of dimension $d$.

We performed a grid search of 10 configurations to find the optimal $\gamma \in \{2, 5, 10, 15, 20, 30, 40, 50, 100, 200\}$ The results are reported in Table 14.

<table><caption>**Table 14** Results of the random search on the hyper-parameters of the FactorVAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>γ</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>40</td><td>0.876 ± 0.003</td><td>1.895 ± 0.084</td></tr><tr><td>100</td><td>0.875 ± 0.007</td><td>1.827 ± 0.092</td></tr><tr><td>15</td><td>0.874 ± 0.004</td><td>1.872 ± 0.048</td></tr><tr><td>20</td><td>0.869 ± 0.010</td><td>1.875 ± 0.090</td></tr><tr><td>50</td><td>0.866 ± 0.011</td><td>1.850 ± 0.070</td></tr><tr><td>200</td><td>0.864 ± 0.008</td><td>1.820 ± 0.032</td></tr><tr><td>10</td><td>0.864 ± 0.011</td><td>1.859 ± 0.086</td></tr><tr><td>30</td><td>0.864 ± 0.020</td><td>1.805 ± 0.096</td></tr><tr><td>5</td><td>0.862 ± 0.019</td><td>1.890 ± 0.075</td></tr><tr><td>2</td><td>0.852 ± 0.016</td><td>1.901 ± 0.081</td></tr></tbody></table>

## D.6 Hamiltonian VAE

Caterini et al.¹⁵ introduced a new method to obtain a low variance unbiased estimation of the ELBO using Markov chain Monte Carlo with Hamiltonian importance sampling⁷⁹ and by proposing a method to select optimal reverse kernels, building the Hamiltonian VAE with the following loss:

$$ \mathcal{L}_{\text{HVAE}} = \mathbb{E}_{z_0 \sim q_{\theta, \phi}^0(\cdot, \cdot)} \left[ \log p_{\theta}(x, z_K) - \frac{1}{2} \rho_K^T \rho_K - \log q_{\theta, \phi}^0(z_0) \right] + \frac{l}{2} $$

where $(z_0, \rho_0) = \mathcal{H}_{\theta,\phi}(z_0, \gamma_0/\sqrt{\beta_0})$, $\mathcal{H}$ is the Hamiltonian importance sampling,⁷⁹ $\beta_0$ is the inverse temperature and $\gamma_0 \sim N(\cdot|0, I)$.

There are three hyper-parameters that we randomly searched for: the number of step in the leapfrog $n_{lf} \in \{1, 2, 10, 15, 20\}$, the leapfrog step size $\epsilon_{lf} \in \{10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}\}$ and $\beta_0 \in \{0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\}$ the tempering factor in the Hamiltonian Monte Carlo Sampler. We trained 20 configurations out of 220 possible combinations. The results are reported in Table 15. Note that some configurations were really long to train, sometimes exceeding the time limit of the HPC used to train the models.
<table><caption>**Table 15** Results of the random search on the hyper-parameters of the HVAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>n<sub>lf</sub></th><th>ε<sub>lf</sub></th><th>β<sub>0</sub></th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>10</td><td>0.00001</td><td>0.8</td><td>0.873 ± 0.007</td><td>1.862 ± 0.068</td></tr><tr><td>2</td><td>0.00001</td><td>0.7</td><td>0.870 ± 0.002</td><td>1.905 ± 0.079</td></tr><tr><td>2</td><td>0.001</td><td>0.2</td><td>0.870 ± 0.004</td><td>1.847 ± 0.082</td></tr><tr><td>1</td><td>0.001</td><td>0.5</td><td>0.869 ± 0.008</td><td>1.853 ± 0.101</td></tr><tr><td>15</td><td>0.00001</td><td>0.2</td><td>0.868 ± 0.009</td><td>1.854 ± 0.075</td></tr><tr><td>1</td><td>0.001</td><td>0.7</td><td>0.865 ± 0.004</td><td>1.890 ± 0.102</td></tr><tr><td>2</td><td>0.01</td><td>1</td><td>0.865 ± 0.005</td><td>1.911 ± 0.066</td></tr><tr><td>15</td><td>0.001</td><td>0.4</td><td>0.865 ± 0.008</td><td>1.805 ± 0.024</td></tr><tr><td>15</td><td>0.001</td><td>0.1</td><td>0.864 ± 0.009</td><td>1.908 ± 0.084</td></tr><tr><td>10</td><td>0.0001</td><td>0.9</td><td>0.863 ± 0.003</td><td>1.882 ± 0.104</td></tr></tbody></table>

## D.7 Info VAE MMD

To improve both the generative model and the amortized inference distribution, Zhao et al.²⁹ proposed to add the mutual information between *z* and *x* in the objective function of the VAE. To be optimized, the loss is rewritten as follows:

$$ \mathcal{L}_{\text{InfoVAE}} = \mathbb{E}_{p(\mathcal{D}(x))} \mathbb{E}_{q_{\phi}(z|x)} [\log p_{\theta}(x|z)] - (1-\alpha) \mathbb{E}_{p(\mathcal{D}(x))} \mathcal{D}_{KL}(q_{\phi}(z|x) || p(z)) - (\alpha+\lambda-1) \mathcal{D}(q_{\phi}(z) || p(z)) $$

with $\mathcal{D}$ the maximum mean discrepancy (MMD).

We performed a random search of the following parameters: $\alpha \in \{0.0, 0.2, 0.4, 0.6, 0.8, 1.0\}$, $\lambda \in \{0.01, 0.1, 1, 10, 100\}$, the choice of the kernel for the MMD $\in \{\text{rbf}, \text{imq}\}$ (rbf: radial basis function, imq: inverse multi-quadratic) and the kernel bandwidth $\in \{0.01, 0.1, 0.5, 1, 5, 10, 100\}$. We trained 30 configurations out of 420 possible combinations. The results are reported in Table 16.

<table><caption>**Table 16** Results of the random search on the hyper-parameters of the InfoVAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>kernel choice</th><th>α</th><th>λ</th><th>kernel bandwidth</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>rbf</td><td>1</td><td>0.1</td><td>0.1</td><td>0.877 ± 0.006</td><td>1.813 ± 0.075</td></tr><tr><td>rbf</td><td>0.4</td><td>1</td><td>0.5</td><td>0.875 ± 0.006</td><td>1.804 ± 0.052</td></tr><tr><td>rbf</td><td>1</td><td>0.1</td><td>0.5</td><td>0.874 ± 0.003</td><td>1.770 ± 0.077</td></tr><tr><td>rbf</td><td>0.00001</td><td>100</td><td>1</td><td>0.873 ± 0.002</td><td>1.852 ± 0.095</td></tr><tr><td>rbf</td><td>0.00001</td><td>10</td><td>0.5</td><td>0.873 ± 0.004</td><td>1.846 ± 0.094</td></tr><tr><td>imq</td><td>0.4</td><td>100</td><td>0.1</td><td>0.872 ± 0.008</td><td>1.866 ± 0.045</td></tr><tr><td>imq</td><td>1</td><td>10</td><td>1</td><td>0.872 ± 0.006</td><td>1.830 ± 0.068</td></tr><tr><td>rbf</td><td>0.6</td><td>0.01</td><td>5</td><td>0.871 ± 0.007</td><td>1.832 ± 0.088</td></tr><tr><td>rbf</td><td>1</td><td>100</td><td>0.01</td><td>0.870 ± 0.004</td><td>1.768 ± 0.079</td></tr><tr><td>imq</td><td>0.2</td><td>0.01</td><td>0.01</td><td>0.870 ± 0.005</td><td>1.830 ± 0.107</td></tr></tbody></table>

## D.8 IWAE

Instead of relying on a single sample for estimating the posterior, the IWAE¹³ utilizes importance weights during the sampling process in the latent space on multiple samples (Monte Carlo
estimator), assigning higher weights to more probable samples. This provides a new ELBO that becomes tighter when the number of samples increases. The loss is the following:

$$ \mathcal{L}_{\text{IWAE}} = \mathbb{E}_{z_1, \dots, z_k \sim q_{\phi}(z|x)} \left[ \log \frac{1}{k} \sum_{i=1}^{k} \frac{p_{\theta}(x, z_i)}{q_{\phi}(z_i|x)} \right] $$

with $k \in \{2, 3, 4, 5, 6, 8, 10, 12, 15, 20\}$ the number of samples to use in the Monte Carlo estimator.

When $k$ grows, the IWAE becomes very memory greedy and time-consuming during training, especially with 3D images. We had to reduce the batch size to 2, and, in spite of this, the model would crash because of memory when setting $k > 6$. The results are reported in Table 17.

<table><caption>**Table 17** Results of the random search on the hyper-parameters of the IWAE: ranking according to the SSIM of the 3 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>number of samples</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>6</td><td>0.865 ± 0.007</td><td>2.087 ± 0.146</td></tr><tr><td>3</td><td>0.861 ± 0.002</td><td>2.048 ± 0.064</td></tr><tr><td>4</td><td>0.854 ± 0.013</td><td>2.178 ± 0.201</td></tr></tbody></table>

## D.9 MS-SSIM VAE

Snell et al.²⁵ proposed an extension of the VAE, called the expected loss VAE, where the pixel-wise reconstruction loss can be replaced by any deterministic reconstruction loss. For this, the probabilistic decoder $p_{\theta}$ is replaced by a deterministic equivalent $f_{\theta}$ so that the reconstruction $\hat{x}$ of $x$ given $z \sim q_{\phi}(z|x)$ is given by $\hat{x} = f_{\theta}(x)$ and the reconstruction loss is given by $\Delta(x, \hat{x})$. The objective becomes

$$ \mathcal{L}_{\text{EL-VAE}} = \mathbb{E}_{q_{\phi}(z|x)} [\Delta(x, \hat{x})] - \beta \mathcal{D}_{\text{KL}} [q_{\phi}(z|x) || p_z(z)] . $$

Following the authors' suggestion, we use the MS-SSIM, or multi-scale structural similarity, as our reconstruction loss.

We performed a random search on $\beta$ and the window size used in the computations of the MS-SSIM, where $\beta$ is sampled from $\{0.01, 0.1, 1, 10, 100\}$ and the window size is sampled from $\{2, 3, 5, 11\}$. We trained 10 configurations out of 20 possible combinations. The results are reported in Table 18.

The training time was too long for configurations with a window size different from 2, explaining why Table 18 contains only five configuration with a window size of 2.
<table><caption>**Table 18** Results of the random search on the hyper-parameters of the MS-SSIM VAE: ranking according to the SSIM of the 5 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>β</th><th>window size</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>100</td><td>2</td><td>0.472 ± 0.034</td><td>70.174 ± 5.660</td></tr><tr><td>1</td><td>2</td><td>0.453 ± 0.050</td><td>75.443 ± 11.098</td></tr><tr><td>1</td><td>2</td><td>0.448 ± 0.018</td><td>74.110 ± 1.158</td></tr><tr><td>1</td><td>2</td><td>0.445 ± 0.050</td><td>76.354 ± 8.802</td></tr><tr><td>0.01</td><td>2</td><td>0.393 ± 0.039</td><td>84.579 ± 7.183</td></tr></tbody></table>

## D.10 Regularized auto-encoder

Ghosh et al.¹⁸ claimed that the probabilistic sampling in VAE is equivalent to a noise injection to the decoder, acting as a stochastic regularization of the latent space. The authors proposed a new approach that consists in replacing the random noise injection by a deterministic regularization in the decoder. The training objective becomes

$$ \mathcal{L}_{\text{RAE}} = \|x - \hat{x}\|_{2}^{2} + \beta \cdot \mathcal{L}_{Z}^{\text{RAE}} + \lambda \cdot \mathcal{L}_{\text{REG}}, $$

with $\mathcal{L}_{\text{REG}}$ the regularization term for the decoder and $\mathcal{L}_{Z}^{\text{RAE}} = 1/2\|z\|_{2}^{2}$ a constraint on the latent space. The authors suggested two different regularization terms for the decoder:

* the first option is a $\mathcal{L}_{2}$ norm on the weights of the decoder $\mathcal{L}_{\text{REG}} = \|\theta\|_{2}^{2}$, giving the RAE-$\ell^{2}$ model;

* another choice is to apply a gradient penalty on the discriminator $\mathcal{L}_{\text{REG}} = \|\nabla D_{\theta}(E_{\phi}(x))\|_{2}^{2}$, giving the RAE-GP model.

We performed a random search on both $\lambda$ and $\beta$, that are both sampled from $\{10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}, 1\}$. We trained 20 configurations for both the RAE-$\ell^{2}$ and the RAE-GP out of 36 possible combinations for each model. The results are respectively reported in Tables 19 and 20.

<table><caption>**Table 19** Results of the random search on the hyper-parameters of the RAE-$\ell^{2}$: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>embedding weight</th><th>reg weight</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>0.0001</td><td>1</td><td>0.884 ± 0.005</td><td>1.815 ± 0.049</td></tr><tr><td>0.0001</td><td>0.001</td><td>0.883 ± 0.002</td><td>1.765 ± 0.070</td></tr><tr><td>0.0001</td><td>1</td><td>0.879 ± 0.008</td><td>1.848 ± 0.059</td></tr><tr><td>0.0001</td><td>0.01</td><td>0.879 ± 0.009</td><td>1.857 ± 0.064</td></tr><tr><td>0.00001</td><td>0.01</td><td>0.879 ± 0.007</td><td>1.868 ± 0.055</td></tr><tr><td>0.1</td><td>0.001</td><td>0.878 ± 0.007</td><td>1.814 ± 0.052</td></tr><tr><td>0.00001</td><td>0.01</td><td>0.878 ± 0.006</td><td>1.785 ± 0.076</td></tr><tr><td>0.1</td><td>0.0001</td><td>0.878 ± 0.007</td><td>1.853 ± 0.077</td></tr><tr><td>0.00001</td><td>0.1</td><td>0.878 ± 0.007</td><td>1.783 ± 0.107</td></tr><tr><td>0.00001</td><td>0.01</td><td>0.877 ± 0.005</td><td>1.831 ± 0.121</td></tr></tbody></table>
<table><caption>**Table 20** Results of the random search on the hyper-parameters of the RAE-GP: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>embedding weight</th><th>reg weight</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>0.01</td><td>0.0001</td><td>0.880 ± 0.006</td><td>1.715 ± 0.105</td></tr><tr><td>0.0001</td><td>0.0001</td><td>0.877 ± 0.008</td><td>1.744 ± 0.056</td></tr><tr><td>0.1</td><td>0.00001</td><td>0.877 ± 0.009</td><td>1.820 ± 0.093</td></tr><tr><td>1</td><td>0.001</td><td>0.867 ± 0.003</td><td>1.756 ± 0.063</td></tr><tr><td>0.1</td><td>0.01</td><td>0.861 ± 0.011</td><td>1.828 ± 0.031</td></tr><tr><td>0.1</td><td>0.1</td><td>0.845 ± 0.012</td><td>1.750 ± 0.107</td></tr><tr><td>0.0001</td><td>0.1</td><td>0.842 ± 0.010</td><td>1.769 ± 0.079</td></tr><tr><td>0.1</td><td>0.1</td><td>0.839 ± 0.008</td><td>1.799 ± 0.101</td></tr><tr><td>0.00001</td><td>1</td><td>0.825 ± 0.013</td><td>1.906 ± 0.135</td></tr><tr><td>0.1</td><td>1</td><td>0.808 ± 0.004</td><td>1.924 ± 0.145</td></tr></tbody></table>

## D.11 Hyperspherical VAE

The hyperspherical VAE¹⁷ uses a von Mises-Fisher (vMF) distribution as prior, leading to a hyperspherical latent space. This model has the advantage of not having additional hyper-parameters compared to a standard VAE but only works with a small latent space as large values lead to errors when computing the modified Bessel function involved in the probability density function of the vMF distribution. Therefore, we performed a grid search on three different smaller latent space sizes: 8, 16, 32. The results are reported in Table 21.

<table><caption>**Table 21** Results of the random search on the hyper-parameters of the SVAE: ranking according to the SSIM of the 3 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>latent space size</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>16</td><td>0.151 ± 0.001</td><td>632.694 ± 5.106</td></tr><tr><td>32</td><td>0.150 ± 0.002</td><td>640.791 ± 4.328</td></tr><tr><td>8</td><td>0.083 ± 0.028</td><td>189.998 ± 68.394</td></tr></tbody></table>

## D.12 VAEGAN

In the VAEGAN,²² a discriminator is trained on the output of a VAE to enhance the VAE's reconstruction abilities. The idea is to use the learned feature representations from intermediate layers of the GAN discriminator as a basis for the VAE reconstruction objective, assuming that the discriminator can capture high-level structures relevant to the data distribution. Overall, this allows replacing voxel-wise similarity between input and output by feature-wise similarity. For $z \sim p_z(z)$ and $\hat{x} \sim D_\theta(z)$, the objective is given by

$$ \mathcal{L}_{\text{VAEGAN}} = \mathbb{E}_{z \sim q_{\phi}(z|x)} [\log \mathcal{N}(D_l(x)|D_l(\hat{x}), \mathbf{I})] - \mathcal{D}_{\text{KL}} [q_{\phi}(z|x) || p_z(z)] - \log \left( \frac{D(x)}{1 - D(D_{\theta}(z))} \right), $$

where $D$ denotes the discriminator, $D_l$ the hidden representation of the $l$-th layer of the discriminator, and $D_\theta$ the decoder. We also added a hyper-parameter $\alpha$ to the decoder's loss, such that high values of $\alpha$ encourage better reconstruction with respect to the features learned at the layer $l$ of the discriminator.
We set the discriminator to be a neural network with 4 convolutions and 2 fully connected layers, with batch normalization and ReLU activation. We set the margin to 0.4 and the equilibrium to 0.68 as in the original paper. We performed a random search for $\alpha \in \{0.3, 0.5, 0.7, 0.9\}$ and $l \in \{1, 2, 3, 4\}$. This model is particularly long to train and, due to memory constraints, we reduced the batch size to 4 instead of 8 for these models. We trained 10 different configurations out of 16 possible combinations. The results are reported in Table 22. We note that there is a strong correlation between the chosen reconstruction layer in the decoder and the quality of the reconstruction.

<table><caption>**Table 22** Results of the random search on the hyper-parameters of the VAEGAN: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>adversarial loss scale</th><th>reconstruction layer</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>0.5</td><td>1</td><td>0.860 ± 0.013</td><td>2.241 ± 0.193</td></tr><tr><td>0.1</td><td>1</td><td>0.851 ± 0.015</td><td>2.671 ± 0.090</td></tr><tr><td>0.1</td><td>1</td><td>0.849 ± 0.006</td><td>2.547 ± 0.356</td></tr><tr><td>0.3</td><td>2</td><td>0.799 ± 0.050</td><td>3.705 ± 0.559</td></tr><tr><td>0.3</td><td>2</td><td>0.780 ± 0.101</td><td>3.968 ± 0.642</td></tr><tr><td>0.9</td><td>3</td><td>0.714 ± 0.148</td><td>9.727 ± 4.423</td></tr><tr><td>0.7</td><td>3</td><td>0.692 ± 0.101</td><td>9.832 ± 1.185</td></tr><tr><td>0.8</td><td>3</td><td>0.572 ± 0.061</td><td>10.336 ± 5.712</td></tr><tr><td>0.9</td><td>4</td><td>0.560 ± 0.113</td><td>24.463 ± 6.516</td></tr></tbody></table>

## D.13 VAE with inverse auto-regressive flows

The VAE with inverse auto-regressive flows²¹ incorporates a series of inverse auto-regressive flows in the encoder, enhancing the flexibility of the learned posterior distribution, and scaling well to high-dimensional latent spaces. We use masked autoencoder for distribution estimation (MADE)⁸⁰ as normalizing flow, as suggested by Kingma et al.²¹ and implemented in Pythae.⁴³

We performed a random search on the following parameters: the number of MADE blocks $\in \{2, 3, 4, 5, 6, 8\}$, the number of hidden layers in the MADE blocks $\in \{2, 3, 4, 5\}$ and the size of the hidden layers $\in \{64, 128, 256\}$. We trained 30 different configurations out of 72 possible combinations. However, we noticed that the performance was really poor when the number of MADE blocks was odd, reducing the possible values for this parameter to even numbers. The results are reported in Table 23.
<table><caption>**Table 23** Results of the random search on the hyper-parameters of the VAE-IAF: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>n MADE blocks</th><th>n hidden in MADE</th><th>hidden size</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>4</td><td>4</td><td>128</td><td>0.823 ± 0.005</td><td>2.272 ± 0.057</td></tr><tr><td>4</td><td>5</td><td>256</td><td>0.823 ± 0.008</td><td>2.248 ± 0.063</td></tr><tr><td>2</td><td>5</td><td>256</td><td>0.823 ± 0.007</td><td>2.220 ± 0.071</td></tr><tr><td>8</td><td>4</td><td>128</td><td>0.822 ± 0.005</td><td>2.282 ± 0.092</td></tr><tr><td>6</td><td>5</td><td>128</td><td>0.820 ± 0.004</td><td>2.259 ± 0.148</td></tr><tr><td>6</td><td>5</td><td>64</td><td>0.820 ± 0.003</td><td>2.403 ± 0.107</td></tr><tr><td>4</td><td>3</td><td>128</td><td>0.819 ± 0.008</td><td>2.260 ± 0.055</td></tr><tr><td>4</td><td>2</td><td>64</td><td>0.818 ± 0.015</td><td>2.331 ± 0.133</td></tr><tr><td>4</td><td>5</td><td>64</td><td>0.817 ± 0.006</td><td>2.359 ± 0.114</td></tr><tr><td>8</td><td>5</td><td>64</td><td>0.816 ± 0.005</td><td>2.546 ± 0.028</td></tr></tbody></table>

## D.14 VAE with linear normalizing flows

The VAE with linear normalizing flows²⁴ enables a better approximation of the posterior distribution $q_{\phi}(z|x)$ using a series of linear normalizing flows, which are invertible transformations. To get the latent vector $z_K$, $z_0$ is sampled from $q_{\phi}(z|x)$ and passes through $K$ linear normalizing flows $f_k$ such that $z_K = f_K \circ \dots \circ f_2 \circ f_1(z_0)$. These flows enable the model to capture complex distributions in the latent space. The authors suggest using a succession of linear flows, and more precisely planar or radial flows, because it is computationally efficient to compute their Jacobian matrix, as needed to compute the loss

$$ \mathcal{L}_{\text{VAElinNF}} = \mathcal{L}_{\text{rec}} + \log q_{\phi}(z_0) - \log q(z_K) - \sum_{k=1}^{K} \log | \det \frac{\partial f_k}{\partial z} | . $$

We tried 10 different configurations of flows ∈ {10P, 10R, 5P, 5R, 5P5R, 5R5P, 5PR, 5RP, 2PR, 2RP}, with P designating a planar flow and R a radial flow. The results are reported in Table 24. We note that configurations with radial flows clearly outperform configurations with planar flows.

<table><caption>**Table 24** Results of the random search on the hyper-parameters of the VAE LinNF: ranking according to the SSIM of the nine best configurations (mean ± std over the three folds randomly selected). P designate a planar flow, R designate a radial flow.</caption><thead><tr><th>flows</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>10R</td><td>0.871 ± 0.001</td><td>1.855 ± 0.125</td></tr><tr><td>5R</td><td>0.860 ± 0.008</td><td>1.897 ± 0.066</td></tr><tr><td>2PR</td><td>0.827 ± 0.005</td><td>2.262 ± 0.135</td></tr><tr><td>5PR</td><td>0.737 ± 0.058</td><td>3.218 ± 0.205</td></tr><tr><td>5P5R</td><td>0.720 ± 0.095</td><td>3.148 ± 0.763</td></tr><tr><td>5RP</td><td>0.716 ± 0.112</td><td>3.829 ± 0.942</td></tr><tr><td>5R5P</td><td>0.708 ± 0.080</td><td>4.652 ± 1.518</td></tr><tr><td>5P</td><td>0.679 ± 0.098</td><td>4.897 ± 1.878</td></tr><tr><td>10P</td><td>0.570 ± 0.064</td><td>8.082 ± 4.095</td></tr></tbody></table>
## D.15 VAE with VampPrior

The VAE with a “Variational Mixture of Posteriors” prior, or VampPrior,²⁷ aims to replace the simple normal prior with a mixture of distributions (e.g. mixture of Gaussians), allowing capturing more complex data distributions. We optimize the following loss:

$$ \mathcal{L}_{\text{VAMP}} = \mathcal{L}_{\text{rec}} - (\log p_{\lambda}(z) - \log q_{\phi}(z|x)) $$

with $p_{\lambda}(z) = \frac{1}{K} \sum_{k=1}^{K} q_{\phi}(z|u_k)$, K being the number of components, and $u_k$ being the “pseudo-input” learned through back-propagation.

We performed a random search on the number of components $K \in \{10, 20, 30, 40, 50\}$ and the number of linear scheduling steps $\in \{0, 20, 40\}$. The results are reported in Table 25.

<table><caption>**Table 25** Results of the random search on the hyper-parameters of the VAMP: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>number components</th><th>linear scheduling steps</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>20</td><td>40</td><td>0.702 ± 0.097</td><td>5.581 ± 0.874</td></tr><tr><td>10</td><td>20</td><td>0.686 ± 0.019</td><td>7.231 ± 5.478</td></tr><tr><td>10</td><td>20</td><td>0.678 ± 0.025</td><td>5.841 ± 0.902</td></tr><tr><td>20</td><td>20</td><td>0.633 ± 0.007</td><td>3.640 ± 0.025</td></tr><tr><td>20</td><td>0</td><td>0.631 ± 0.003</td><td>3.569 ± 0.108</td></tr><tr><td>20</td><td>20</td><td>0.628 ± 0.002</td><td>3.586 ± 0.120</td></tr><tr><td>30</td><td>20</td><td>0.625 ± 0.005</td><td>3.892 ± 0.150</td></tr><tr><td>30</td><td>0</td><td>0.622 ± 0.001</td><td>3.965 ± 0.091</td></tr><tr><td>40</td><td>40</td><td>0.621 ± 0.005</td><td>4.151 ± 0.226</td></tr><tr><td>40</td><td>0</td><td>0.620 ± 0.004</td><td>4.074 ± 0.169</td></tr></tbody></table>

## D.16 Vector-quantized VAE

Van Den Oord et al.²⁸ suggested using discrete (rather than continuous) latent representations and having a learned (rather than static) prior. The latent space is structured as an $\mathbb{R}^{K \times D}$ vector space. We denote $\mathcal{E} = \{e_1, e_2, ..., e_K\}$ where $e_i \in \mathbb{R}^D$ for $i \in \{1, 2, ..., K\}$. We say that K is the size of the latent embedding space, and D is the dimension of the embedding vectors.

For an embedding size d, the input x is passed through the encoder to obtain the output $z_e(x) \in \mathbb{R}^{d \times D}$, which is then passed through the discretization bottleneck to map it to an element of $z_q(x) \in \mathcal{E}^d$ such that $(z_q(x))_j = e_k$ where $k = \arg \min_l ||z_e(x) - e_l||_2$ for $j \in \{1, 2, ..., d\}$. As the argmin operation lacks differentiability, learning of the embeddings and regularization of the latent space is managed by integrating the stopgradient operator sg into the training objective:

$$ \mathcal{L}_{\text{VQVAE}}(x) = \log p(x|z_q(x)) + ||\text{sg}[z_e(x)] - e||_2^2 + \beta ||z_e(x) - \text{sg}[e]||_2^2 . $$

As suggested by the authors, we replaced the term $||\text{sg}[z_e(x)] - e||_2^2$ in the loss by the exponential moving average (EMA) update with a decay of 0.99. We then considered two hyper-parameters in our random search: the size of the latent embedding space $K \in \{128, 256, 512, 1024\}$ and the regularization weight $\beta \in \{0.25, 0.5, 0.75, 0.9, 1, 1.5, 2, 4\}$. The results are reported in Table 26.
<table><caption>**Table 26** Results of the random search on the hyper-parameters of the VQVAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>commitment loss factor</th><th>quantization loss factor</th><th>num embeddings</th><th>use EMA</th><th>decay</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>0.25</td><td>2</td><td>512</td><td>True</td><td>0.99</td><td>0.881 ± 0.003</td><td>1.805 ± 0.032</td></tr><tr><td>0.25</td><td>0.25</td><td>1024</td><td>True</td><td>0.99</td><td>0.880 ± 0.009</td><td>1.866 ± 0.064</td></tr><tr><td>0.25</td><td>0.5</td><td>256</td><td>True</td><td>0.99</td><td>0.879 ± 0.005</td><td>1.797 ± 0.037</td></tr><tr><td>0.25</td><td>0.25</td><td>512</td><td>True</td><td>0.99</td><td>0.877 ± 0.007</td><td>1.836 ± 0.093</td></tr><tr><td>0.25</td><td>4</td><td>1024</td><td>True</td><td>0.99</td><td>0.876 ± 0.005</td><td>1.854 ± 0.065</td></tr><tr><td>0.25</td><td>0.75</td><td>256</td><td>True</td><td>0.99</td><td>0.874 ± 0.011</td><td>1.896 ± 0.065</td></tr><tr><td>0.25</td><td>0.9</td><td>512</td><td>True</td><td>0.99</td><td>0.870 ± 0.004</td><td>1.927 ± 0.023</td></tr><tr><td>0.25</td><td>1.5</td><td>256</td><td>True</td><td>0.99</td><td>0.870 ± 0.011</td><td>1.856 ± 0.100</td></tr><tr><td>0.25</td><td>4</td><td>1024</td><td>True</td><td>0.99</td><td>0.870 ± 0.011</td><td>1.854 ± 0.056</td></tr><tr><td>0.25</td><td>1.5</td><td>1024</td><td>True</td><td>0.99</td><td>0.868 ± 0.014</td><td>1.827 ± 0.084</td></tr></tbody></table>

## D.17 Wasserstein auto-encoder

The Wasserstein auto-encoder²⁶ introduces the use of a penalized form of the Wasserstein distance to measure the dissimilarity between the model's generated distribution and the true data distribution. This leads to more stable training, mitigating mode collapse and improving the model's ability to generate diverse and realistic samples.

$$ \mathcal{L}_{\text{WAE}} = \mathcal{L}_{\text{rec}} + \lambda \mathcal{D}_Z(p_z(z), q_{\phi}(z)) , $$

with $\mathcal{D}_Z$ an arbitrary divergence. Different divergences $\mathcal{D}_Z$ are suggested by the authors, we here use the maximum mean discrepancy (MMD).

We performed a random search on three parameters: the kernel choice ∈ {rbf, imq}, the regularization weight $\lambda$ ∈ {0.01, 0.1, 0.5, 1, 5, 10, 100} and the kernel bandwidth ∈ {0.01, 0.1, 0.5, 1, 5, 10, 100}. The results are reported in Table 27.

<table><caption>**Table 27** Results of the random search on the hyper-parameters of the WAE: ranking according to the SSIM of the 10 best configurations (mean ± std over the three folds randomly selected).</caption><thead><tr><th>kernel choice</th><th>regularization weight</th><th>kernel bandwidth</th><th>SSIM ↑</th><th>MSE (×10<sup>-3</sup>) ↓</th></tr></thead><tbody><tr><td>rbf</td><td>0.1</td><td>5</td><td>0.881 ± 0.005</td><td>1.862 ± 0.075</td></tr><tr><td>rbf</td><td>0.5</td><td>0.1</td><td>0.880 ± 0.002</td><td>1.835 ± 0.096</td></tr><tr><td>rbf</td><td>0.5</td><td>0.5</td><td>0.879 ± 0.004</td><td>1.798 ± 0.035</td></tr><tr><td>rbf</td><td>0.01</td><td>0.1</td><td>0.879 ± 0.006</td><td>1.866 ± 0.061</td></tr><tr><td>imq</td><td>5</td><td>1</td><td>0.878 ± 0.003</td><td>1.838 ± 0.073</td></tr><tr><td>rbf</td><td>10</td><td>5</td><td>0.877 ± 0.005</td><td>1.838 ± 0.090</td></tr><tr><td>imq</td><td>1</td><td>1</td><td>0.876 ± 0.006</td><td>1.835 ± 0.097</td></tr><tr><td>imq</td><td>1</td><td>0.01</td><td>0.876 ± 0.007</td><td>1.882 ± 0.067</td></tr><tr><td>imq</td><td>5</td><td>100</td><td>0.874 ± 0.002</td><td>1.894 ± 0.070</td></tr><tr><td>imq</td><td>100</td><td>100</td><td>0.874 ± 0.016</td><td>1.865 ± 0.069</td></tr></tbody></table>
E Model training

<table><caption>**Table 28** SSIM obtained on each validation set of the 6-fold cross-validation for the different trained models (mean ± std over the images from the validation set). The best split of each VAE variant is highlighted in bold.</caption><thead><tr><th>Models</th><th>Split 0</th><th>Split 1</th><th>Split 2</th><th>Split 3</th><th>Split 4</th><th>Split 5</th></tr></thead><tbody><tr><td>AE</td><td>0.865 ± 0.032</td><td>0.860 ± 0.026</td><td><strong>0.876 ± 0.036</strong></td><td>0.876 ± 0.024</td><td>0.875 ± 0.030</td><td>0.859 ± 0.040</td></tr><tr><td>Adv. AE<sup>23</sup></td><td>0.866 ± 0.034</td><td>0.852 ± 0.032</td><td>0.868 ± 0.034</td><td><strong>0.881 ± 0.024</strong></td><td>0.871 ± 0.029</td><td>0.862 ± 0.037</td></tr><tr><td>β-TC VAE<sup>16</sup></td><td>0.858 ± 0.036</td><td>0.869 ± 0.024</td><td>0.869 ± 0.030</td><td><strong>0.876 ± 0.024</strong></td><td>0.860 ± 0.035</td><td>0.871 ± 0.037</td></tr><tr><td>β-VAE<sup>19</sup></td><td>0.870 ± 0.029</td><td>0.869 ± 0.024</td><td><strong>0.874 ± 0.030</strong></td><td>0.872 ± 0.027</td><td>0.865 ± 0.032</td><td>0.864 ± 0.036</td></tr><tr><td>Dis. β-VAE<sup>14</sup></td><td>0.870 ± 0.029</td><td>0.868 ± 0.023</td><td><strong>0.879 ± 0.027</strong></td><td>0.868 ± 0.026</td><td>0.865 ± 0.033</td><td>0.865 ± 0.036</td></tr><tr><td>FactorVAE<sup>20</sup></td><td>0.863 ± 0.033</td><td>0.849 ± 0.024</td><td><strong>0.874 ± 0.033</strong></td><td>0.872 ± 0.028</td><td>0.873 ± 0.027</td><td>0.861 ± 0.040</td></tr><tr><td>HVAE<sup>15</sup></td><td>0.840 ± 0.040</td><td>0.869 ± 0.025</td><td>0.864 ± 0.031</td><td>0.867 ± 0.027</td><td><strong>0.878 ± 0.027</strong></td><td>0.858 ± 0.038</td></tr><tr><td>InfoVAE<sup>29</sup></td><td>0.870 ± 0.030</td><td>0.871 ± 0.025</td><td>0.874 ± 0.031</td><td>0.873 ± 0.024</td><td><strong>0.876 ± 0.027</strong></td><td>0.870 ± 0.036</td></tr><tr><td>IWAE<sup>13</sup></td><td>0.861 ± 0.036</td><td>0.860 ± 0.031</td><td>0.867 ± 0.033</td><td>0.868 ± 0.028</td><td><strong>0.875 ± 0.026</strong></td><td>0.861 ± 0.046</td></tr><tr><td>RAE-GP<sup>18</sup></td><td>0.878 ± 0.030</td><td>0.872 ± 0.028</td><td><strong>0.884 ± 0.029</strong></td><td>0.884 ± 0.029</td><td>0.880 ± 0.026</td><td>0.873 ± 0.033</td></tr><tr><td>RAE-ℓ<sup>218</sup></td><td>0.857 ± 0.037</td><td>0.873 ± 0.025</td><td>0.850 ± 0.040</td><td><strong>0.882 ± 0.023</strong></td><td>0.868 ± 0.031</td><td>0.863 ± 0.038</td></tr><tr><td>VAE<sup>4</sup></td><td>0.866 ± 0.030</td><td>0.868 ± 0.024</td><td><strong>0.869 ± 0.030</strong></td><td>0.865 ± 0.029</td><td>0.851 ± 0.041</td><td>0.868 ± 0.037</td></tr><tr><td>VAEGAN<sup>22</sup></td><td>0.804 ± 0.044</td><td>0.863 ± 0.025</td><td>0.846 ± 0.038</td><td><strong>0.866 ± 0.026</strong></td><td>0.855 ± 0.038</td><td>0.858 ± 0.035</td></tr><tr><td>VAE-IAF<sup>21</sup></td><td>0.827 ± 0.037</td><td>0.818 ± 0.032</td><td><strong>0.829 ± 0.034</strong></td><td>0.828 ± 0.027</td><td>0.828 ± 0.034</td><td>0.823 ± 0.037</td></tr><tr><td>VAE LinNF<sup>24</sup></td><td>0.860 ± 0.033</td><td><strong>0.870 ± 0.024</strong></td><td>0.865 ± 0.035</td><td>0.852 ± 0.032</td><td>0.858 ± 0.039</td><td>0.870 ± 0.036</td></tr><tr><td>VQVAE<sup>28</sup></td><td>0.857 ± 0.036</td><td>0.852 ± 0.027</td><td>0.869 ± 0.031</td><td><strong>0.883 ± 0.025</strong></td><td>0.868 ± 0.031</td><td>0.864 ± 0.037</td></tr><tr><td>WAE<sup>26</sup></td><td>0.870 ± 0.032</td><td>0.871 ± 0.024</td><td>0.871 ± 0.033</td><td><strong>0.882 ± 0.024</strong></td><td>0.869 ± 0.034</td><td>0.862 ± 0.038</td></tr></tbody></table>

F Reconstructions
<figure><img src="image_14.png" /><figcaption>**Fig 10** Examples of reconstructions (coronal slices) obtained with the different VAE variants from the original image of a cognitively normal subject (images of the first column, Test CN) and from the same subject with AD simulated at different intensity degrees (AD 15, AD 30, AD 50 and AD 70). The first row shows the input image in odd columns and the mask of the simulated disease in even columns when the input is a simulated image. All the other rows are the pseudo-healthy reconstructions of the models in odd columns and the difference between the pseudo-healthy reconstruction and the input in even columns.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig 11** Examples of reconstructions (sagittal slices) obtained with the different VAE variants from the original image of a cognitively normal subject (images of the first column, Test CN) and from the same subject with AD simulated at different intensity degrees (AD 15, AD 30, AD 50 and AD 70). The first row shows the input image in odd columns and the mask of the simulated disease in even columns when the input is a simulated image. All the other rows are the pseudo-healthy reconstructions of the models in odd columns and the difference between the pseudo-healthy reconstruction and the input in even columns.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig 12** Examples of reconstructions (axial slices) obtained with the different VAE variants from the same subject with different dementia subtypes simulated at 30% intensity degree (bvFTD 30, IvPPA 30, svPPA 30, nfVPPA 30 and PCA 30). The first row shows the input image in odd columns and the mask of the simulated disease in even columns. All the other rows are the pseudo-healthy reconstructions of the models in odd columns and the difference between the pseudo-healthy reconstruction and the input in even columns.</figcaption></figure>
## Disclosures

Competing financial interests related to the present article: none to disclose for all authors.

Competing financial interests unrelated to the present article: Olivier Colliot reports a relationship with Therapanacea that includes consulting or advisory. Olivier Colliot has patent ##PCT/IB2016/0526993 issued to Schiratti J-B, Allassonniere S, Colliot O, Durrleman S.

## Code, Data, and Materials Availability

Most of the code that we used is available in ClinicaDL.⁶² The VAE-based models are implemented in Pythae.⁴³ The preprocessing pipelines are available in Clinica.⁵⁹ The code for random searches, model training and evaluation is available in the following repository: github.com/ravih18/UAD_VAE_benchmark. This repository includes dependencies and software versions used.

Data used in preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu).

## Acknowledgments

The research leading to these results has received funding from the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program (references ANR-19-P3IA-0001 - PRAIRIE - and reference ANR-10IAIHU-06 - IHU ICM) and the "France 2030" program (reference ANR-23-IACL-0008, PRAIRIE-PSAI), and from the European Union's Horizon Europe Framework Program (grant number 101136607, project CLARA). This research was also funded by the Agence Nationale de la Recherche under the project "ANR-23-CE45-0005-01".

This work was granted access to the HPC resources of IDRIS under the allocation AD011011648 made by GENCI (Grand Equipement National de Calcul Intensif).

Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer's Association; Alzheimer's Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.
References

1 T. Fernando, H. Gammulle, S. Denman, et al., “Deep learning for medical anomaly detection – a survey,” *ACM Computing Surveys* **54**(7) (2021).

2 X. Chen and E. Konukoglu, *Unsupervised abnormality detection in medical images with deep generative methods*, 303–324. Elsevier (2022).

3 C. Baur, S. Denner, B. Wiestler, et al., “Autoencoders for unsupervised anomaly segmentation in brain MR images: A comparative study,” *Medical Image Analysis* **69**, 101952 (2021).

4 D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” in *International Conference on Learning Representations (ICLR)*, (2014).

5 I. Goodfellow, J. Pouget-Abadie, M. Mirza, et al., “Generative adversarial nets,” in *Advances in Neural Information Processing Systems*, **27** (2014).

6 J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” *Advances in Neural Information Processing Systems* **33**, 6840–6851 (2020).

7 M. S. Graham, W. H. L. Pinaya, P. Wright, et al., “Unsupervised 3d out-of-distribution detection with latent diffusion models,” in *International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI)*, 446–456, Springer (2023).

8 T. Schlegl, P. Seeböck, S. M. Waldstein, et al., “Unsupervised anomaly detection with generative adversarial networks to guide marker discovery,” in *Information Processing in Medical Imaging, LNCS* (2017).

9 T. Schlegl, P. Seeböck, S. M. Waldstein, et al., “f-AnoGAN: Fast unsupervised anomaly detection with generative adversarial networks,” *Medical Image Analysis* **54** (2019).

10 T. Xia, A. Chartsias, and S. A. Tsaftaris, “Pseudo-healthy synthesis with pathology disentanglement and adversarial learning,” *Medical Image Analysis* **64**, 101719 (2020).

11 R. Shi, C. Sheng, S. Jin, et al., “Generative adversarial network constrained multiple loss autoencoder: A deep learning-based individual atrophy detection for alzheimer's disease and mild cognitive impairment,” *Human brain mapping* **44**(3), 1129–1146 (2023).

12 C. I. Bercea, B. Wiestler, D. Rueckert, et al., “Reversing the abnormal: Pseudo-healthy generative networks for anomaly detection,” *arXiv preprint arXiv:2303.08452* (2023).

13 Y. Burda, R. B. Grosse, and R. Salakhutdinov, “Importance weighted autoencoders,” in *International Conference on Learning Representations (ICLR)*, (2016).

14 C. P. Burgess, I. Higgins, A. Pal, et al., “Understanding disentangling in β-vae,” *arXiv preprint arXiv:1804.03599* (2018).

15 A. L. Caterini, A. Doucet, and D. Sejdinovic, “Hamiltonian Variational Auto-Encoder,” in *Advances in Neural Information Processing Systems*, **31**, Curran Associates, Inc. (2018).

16 R. T. Chen, X. Li, R. B. Grosse, et al., “Isolating sources of disentanglement in variational autoencoders,” *Advances in Neural Information Processing Systems* **31** (2018).

17 T. R. Davidson, L. Falorsi, N. De Cao, et al., “Hyperspherical variational auto-encoders,” *arXiv:1804.00891* (2018).

18 P. Ghosh, M. S. Sajjadi, A. Vergari, et al., “From variational to deterministic autoencoders,” *arXiv:1903.12436* (2019).

19 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 (ICLR)*, (2017).

20 H. Kim and A. Mnih, “Disentangling by factorising,” in *International conference on machine learning (ICML)*, 2649–2658, PMLR (2018).
21 D. P. Kingma, T. Salimans, R. Jozefowicz, et al., "Improved variational inference with inverse autoregressive flow," *Advances in Neural Information Processing Systems* 29 (2016).

22 A. B. L. Larsen, S. K. Sønderby, H. Larochelle, et al., "Autoencoding beyond pixels using a learned similarity metric," in *International conference on machine learning (ICML)*, 1558–1566, PMLR (2016).

23 A. Makhzani, J. Shlens, N. Jaitly, et al., "Adversarial autoencoders," arXiv:1511.05644 (2015).

24 D. Rezende and S. Mohamed, "Variational inference with normalizing flows," in *International conference on machine learning (ICML)*, 1530–1538, PMLR (2015).

25 J. Snell, K. Ridgeway, R. Liao, et al., "Learning to generate images with perceptual similarity metrics," in *ICIP*, 4277–4281, IEEE (2017).

26 I. Tolstikhin, O. Bousquet, S. Gelly, et al., "Wasserstein auto-encoders," in *International Conference on Learning Representations (ICLR)*, (2018).

27 J. Tomczak and M. Welling, "VAE with a VampPrior," in *International Conference on Artificial Intelligence and Statistics*, 1214–1223, PMLR (2018).

28 A. Van Den Oord, O. Vinyals, et al., "Neural discrete representation learning," *Advances in Neural Information Processing Systems* 30 (2017).

29 S. Zhao, J. Song, and S. Ermon, "Infovae: Balancing learning and inference in variational autoencoders," in *Proc AAAI conference on artificial intelligence*, 33, 5885–5892 (2019).

30 X. Chen and E. Konukoglu, "Unsupervised detection of lesions in brain MRI using constrained adversarial auto-encoders," in *MIDL*, (2018).

31 H. Choi, S. Ha, H. Kang, et al., "Deep learning only by normal brain PET identify unheralded brain anomalies," *EBioMedicine* 43, 447–453 (2019).

32 M. Mostapha, J. Prieto, V. Murphy, et al., "Semi-supervised vae-gan for out-of-sample detection applied to mri quality control," in *International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI)*, 127–136, Springer (2019).

33 H. Uzunova, S. Schultz, H. Handels, et al., "Unsupervised pathology detection in medical images using conditional variational autoencoders," *IJCARS* 14, 451–461 (2019).

34 D. Zimmerer, F. Isensee, J. Petersen, et al., "Unsupervised anomaly localization using variational auto-encoders," in *International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI)*, 289–297, Springer (2019).

35 C. Baur, B. Wiestler, S. Albarqouni, et al., "Deep autoencoding models for unsupervised anomaly segmentation in brain mr images," in *Brainlesion: Glioma, Multiple Sclerosis, Stroke and Traumatic Brain Injuries: 4th International Workshop, BrainLes 2018, Held in Conjunction with MICCAI 2018, Granada, Spain, September 16, 2018, Revised Selected Papers, Part I 4*, 161–169, Springer (2019).

36 X. Chen, S. You, K. C. Tezcan, et al., "Unsupervised lesion detection via image restoration with a normative prior," *Medical image analysis* 64, 101713 (2020).

37 G. Chételat, J. Arbizu, H. Barthel, et al., "Amyloid-PET and 18F-FDG-PET in the diagnostic investigation of Alzheimer's disease and other dementias," *The Lancet Neurology* 19(11), 951–962 (2020).

38 K. Herholz, "FDG PET and differential diagnosis of dementia," *Alzheimer Disease and Associated Disorders* 9(1), 6–16 (1995).

39 H. B. Baydargil, J.-S. Park, and D.-Y. Kang, "Anomaly Analysis of Alzheimer's Disease in PET Images Using an Unsupervised Adversarial Deep Learning Model," *Applied Sciences* 11(5), 2187 (2021).
40 C. Hinge, O. M. Henriksen, U. Lindberg, et al., "A zero-dose synthetic baseline for the personalized analysis of [18F]FDG-PET: Application in Alzheimer's disease," *Frontiers in Neuroscience* **16** (2022).

41 C. R. Jack, D. A. Bennett, K. Blennow, et al., "A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers," *Neurology* **87**(5), 539-547 (2016).

42 N. Burgos, M. J. Cardoso, J. Samper-González, et al., "Anomaly detection for the individual analysis of brain PET images," *J Med Imag* **8**(2), 024003 (2021).

43 C. Chadebec, L. J. Vincent, and S. Allassonniere, "Pythae: Unifying generative autoencoders in python - a benchmarking use case," in *Thirty-sixth Conference on NeurIPS Datasets and Benchmarks Track*, (2022).

44 R. Hassanaly, C. Brianceau, M. Solal, et al., "Evaluation of pseudo-healthy image reconstruction for anomaly detection with deep generative models: Application to brain FDG PET," *Machine Learning for Biomedical Imaging* **2**, 611-656 (2024).

45 R. Hassanaly, C. Brianceau, O. Colliot, et al., "Unsupervised anomaly detection in 3D brain FDG PET: A benchmark of 17 vae-based approaches," in *Deep Generative Models workshop at the 26th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI)*, (2023).

46 Z. Wang, A. C. Bovik, H. R. Sheikh, et al., "Image quality assessment: from error visibility to structural similarity," *IEEE T Image Process* **13**(4), 600-612 (2004).

47 T. Nečasová, N. Burgos, and D. Svoboda, "Validation and evaluation metrics for medical and biomedical image synthesis," in *Biomedical Image Synthesis and Simulation*, N. Burgos and D. Svoboda, Eds., *The MICCAI Society book Series*, 573-600, Academic Press (2022).

48 T. Xia, A. Chartsias, and S. A. Tsaftaris, "Adversarial pseudo healthy synthesis needs pathology factorization," in *International Conference on Medical Imaging with Deep Learning*, 512-526, PMLR (2019).

49 C. I. Bercea, B. Wiestler, D. Rueckert, et al., "Generalizing unsupervised anomaly detection: Towards unbiased pathology screening," in *Medical Imaging with Deep Learning*, (2023).

50 C. T. Lüth, D. Zimmerer, G. Koehler, et al., "Cradl: Contrastive representations for unsupervised anomaly detection and localization," in *Bildverarbeitung für die Medizin* **2023**, (2023).

51 V. Wargnier-Dauchelle, T. Grenier, F. Durand-Dubief, et al., "A weakly supervised gradient attribution constraint for interpretable classification and anomaly detection," *IEEE Transactions on Medical Imaging* (2023).

52 W. H. Pinaya, P.-D. Tudosiu, R. Gray, et al., "Unsupervised brain imaging 3d anomaly detection and segmentation with transformers," *Medical Image Analysis* **79**, 102475 (2022).

53 S. Chatterjee, A. Sciarra, M. Dünnwald, et al., "Strega: Unsupervised anomaly detection in brain mris using a compact context-encoding variational autoencoder," *Computers in Biology and Medicine* **149**, 106093 (2022).

54 G. Luo, W. Xie, R. Gao, et al., "Unsupervised anomaly detection in brain mri: Learning abstract distribution from massive healthy brains," *Computers in Biology and Medicine* **154**, 106610 (2023).

55 M. Bengs, F. Behrendt, M.-H. Laves, et al., "Unsupervised anomaly detection in 3d brain mri using deep learning with multi-task brain age prediction," in *Medical Imaging 2022: Computer-Aided Diagnosis*, **12033**, 291-295, SPIE (2022).

56 R. Hassanaly, S. Bottani, B. Sauty, et al., "Simulation based evaluation framework for deep learning unsupervised anomaly detection on brain fdg-pet," in *Medical Imaging 2023: Image Processing*, **12464**, 511-518, SPIE (2023).
57 W. J. Jagust, D. Bandy, K. Chen, et al., "The Alzheimer's Disease Neuroimaging Initiative positron emission tomography core," *Alzheimer's & Dementia* **6**(3), 221–229 (2010).

58 W. J. Jagust, S. M. Landau, R. A. Koeppe, et al., "The Alzheimer's Disease Neuroimaging Initiative 2 PET Core: 2015," *Alzheimer's & Dementia* **11**(7), 757–771 (2015).

59 A. Routier, N. Burgos, M. Díaz, et al., "Clinica: An open-source software platform for reproducible clinical neuroscience studies," *Front Neuroinform* **15** (2021).

60 V. Fonov, A. Evans, R. McKinstry, et al., "Unbiased nonlinear average age-appropriate brain templates from birth to adulthood," *NeuroImage* **47**, S102 (2009).

61 V. Fonov, A. C. Evans, K. Botteron, et al., "Unbiased average age-appropriate atlases for pediatric studies," *NeuroImage* **54**(1), 313–327 (2011).

62 E. Thibeau-Sutre, M. Díaz, R. Hassanaly, et al., "ClinicaDL: An open-source deep learning software for reproducible neuroimaging processing," *Computer Methods and Programs in Biomedicine* **220** (2022).

63 J. Wen, E. Thibeau-Sutre, M. Diaz-Melo, et al., "Convolutional neural networks for classification of alzheimer's disease: Overview and reproducible evaluation," *Medical Image Analysis* **63**, 101694 (2020).

64 X. Bouthillier, P. Delaunay, M. Bronzi, et al., "Accounting for variance in machine learning benchmarks," *Proceedings of Machine Learning and Systems* **3** (2021).

65 A. Vahdat and J. Kautz, "Nvae: A deep hierarchical variational autoencoder," *Advances in Neural Information Processing Systems* **33**, 19667–19679 (2020).

66 K. He, X. Zhang, S. Ren, et al., "Deep residual learning for image recognition," in *Proceedings of the IEEE conference on computer vision and pattern recognition*, 770–778 (2016).

67 K. Simonyan and A. Zisserman, "Very deep convolutional networks for large-scale image recognition," *arXiv preprint arXiv:1409.1556* (2014).

68 Z. Wang, E. P. Simoncelli, and A. C. Bovik, "Multiscale structural similarity for image quality assessment," in *The Thrifty-Seventh Asilomar Conference on Signals, Systems & Computers*, 2003, **2**, 1398–1402, Ieee (2003).

69 C. Oh, K. Adamczewski, and M. Park, "Radial and directional posteriors for bayesian neural networks," *arXiv preprint arXiv:1902.02603* (2019).

70 R. Cuingnet, "On the computation of the logarithm of the modified bessel function of the second kind," *arXiv preprint arXiv:2308.11964* (2023).

71 M. Solal, R. Hassanaly, and N. Burgos, "Leveraging healthy population variability in deep learning unsupervised anomaly detection in brain FDG PET," in *SPIE Medical Imaging*, (San Diego (California), United States) (2024).

72 S. Kumar, P. R. Payne, and A. Sotiras, "Normative modeling using multimodal variational autoencoders to identify abnormal brain volume deviations in alzheimer's disease," in *Medical Imaging 2023: Computer-Aided Diagnosis*, **12465**, 1246503, SPIE (2023).

73 C. I. Bercea, B. Wiestler, D. Rueckert, et al., "Evaluating normative representation learning in generative ai for robust anomaly detection in brain imaging," *Nature Communications* **16**(1), 1624 (2025).

74 F. Meissen, B. Wiestler, G. Kaissis, et al., "On the pitfalls of using the residual error as anomaly score," in *International Conference on Medical Imaging with Deep Learning*, 914–928, PMLR (2022).

75 F. Angiulli, F. Fassetti, and L. Ferragina, "Improving deep unsupervised anomaly detection by exploiting vae latent space distribution," in *International Conference on Discovery Science*, 596–611, Springer (2020).
76 Z. Sun, P. Wang, W. Zheng, et al., "Dual groupgan: An unsupervised four-competitor (2v2) approach for video anomaly detection," *Pattern Recognition* **153**, 110500 (2024).

77 O. Colliot, E. Thibeau-Sutre, and N. Burgos, "Reproducibility in Machine Learning for Medical Imaging," in *Machine Learning for Brain Disorders*, O. Colliot, Ed., *Neuromethods*, 631-653, Springer US, New York, NY (2023).

78 K. J. Gorgolewski, T. Auer, V. D. Calhoun, et al., "The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments," *Scientific data* **3**, 160044 (2016).

79 R. M. Neal, "Hamiltonian importance sampling," in *talk presented at the Banff International Research Station (BIRS) workshop on Mathematical Issues in Molecular Dynamics*, (2005).

80 M. Germain, K. Gregor, I. Murray, et al., "Made: Masked autoencoder for distribution estimation," in *International conference on machine learning (ICML)*, 881-889, PMLR (2015).