Numerical Uncertainty of Convolutional Neural Networks Inference for
  Structural Brain MRI Analysis
This paper investigates the numerical uncertainty of Convolutional Neural
Networks (CNNs) inference for structural brain MRI analysis. It applies Random
Rounding -- a stochastic arithmetic technique -- to CNN models employed in
non-linear registration (SynthMorph) and whole-brain segmentation (FastSurfer),
and compares the resulting numerical uncertainty to the one measured in a
reference image-processing pipeline (FreeSurfer recon-all). Results obtained on
32 representative subjects show that CNN predictions are substantially more
accurate numerically than traditional image-processing results (non-linear
registration: 19 vs 13 significant bits on average; whole-brain segmentation:
0.99 vs 0.92 S{\o}rensen-Dice score on average), which suggests a better
reproducibility of CNN results across execution environments.
Inés Gonzalez Pepe, Vinuyan Sivakolunthu, Hae Lang Park, Yohan Chatelain, Tristan Glatard
# Numerical Uncertainty of Convolutional Neural Networks Inference for Structural Brain MRI Analysis

###### Abstract

This paper investigates the numerical uncertainty of Convolutional Neural Networks (CNNs) inference for structural brain MRI analysis. It applies Random Rounding--a stochastic arithmetic technique--to CNN models employed in non-linear registration (SynthMorph) and whole-brain segmentation (FastSurfer), and compares the resulting numerical uncertainty to the one measured in a reference image-processing pipeline (FreeSurfer recon-all). Results obtained on 32 representative subjects show that CNN predictions are substantially more accurate numerically than traditional image-processing results (non-linear registration: 19 vs 13 significant bits on average; whole-brain segmentation: 0.99 vs 0.92 Sorensen-Dice score on average), which suggests a better reproducibility of CNN results across execution environments.

Keywords:Numerical Stability Convolutional Neural Networks Nonlinear Registration Whole-Brain Segmentation

## 1 Introduction

A motivating factor to study numerical uncertainty in neuroimaging is to establish measures of reliability in the tools observed, particularly in light of the reproducibility crisis [1, 2, 3]. Numerical uncertainty is key to the robustness of neuroimaging analyses. Small computational perturbations introduced in execution environments-- including operating systems, hardware architecture, and parallelization--may amplify throughout analytical pipelines and result in substantial differences in the final outcome of analyses [4, 5]. Such instabilities have been observed across many different tools and imaging modalities [6, 7], and are likely to impact the reproducibility and robustness of analyses.

Convolutional Neural Networks (CNNs) are increasingly adopted for registration [8, 9, 10] and segmentation [11, 12, 13, 14] of structural MRIs. Once trained, CNNs are orders of magnitude faster than traditional image-processing methods, achieve comparable accuracy, and seem to exhibit better generalizability to image modalities and orientations. However, the numerical uncertainty associated with CNN predictions in neuroimaging remains largely unexplored. While previous
works suggested that CNNs might be subject to numerical instability [15, 16, 17], it is unclear how such instabilities manifest in specific CNN architectures used in structural brain MRI, and how the resulting numerical uncertainty compares to the one of traditional methods.

This paper measures the numerical uncertainty associated with CNN inference in neuroimaging, focusing specifically on non-linear registration and whole-brain segmentation of structural MRIs. To do so, it applies Random Rounding (RR) [18, 19]--a practical stochastic arithmetic technique to estimate numerical uncertainty--to state-of-the-art CNN models SynthMorph [8] and FastSurfer [12], and compare their numerical uncertainty to the one measured from the FreeSurfer [20] "recon-all" reference neuroimaging tool.

## 2 Materials & Methods

We measured the numerical uncertainty of CNN models SynthMorph (non-linear registration) and FastSurfer (whole-brain segmentation) using RR. We applied these models to 35 subjects randomly selected in the CoRR dataset, using the FreeSurfer recon-all pipeline as a baseline for numerical uncertainty comparison.

### Random Rounding

Random Rounding (RR) [18] is a form of Monte-Carlo Arithmetic (MCA) [21] that simulates rounding errors by applying the following perturbation to all floating-point (FP) operations of an application:

\[random\_rounding(x\circ y)=round(inexact(x\circ y))\]

where \(x\) and \(y\) are FP numbers, \(\circ\) is an arithmetic operation, and \(inexact\) is a random perturbation defined at a given virtual precision:

\[inexact(x)=x+2^{e_{x}-t}\xi\]

where \(e_{x}\) is the exponent in the FP representation of \(x\), \(t\) is the virtual precision, and \(\xi\) is a random uniform variable of \((-\frac{1}{2},\frac{1}{2})\). To measure numerical uncertainty, we applied a perturbation of 1 ulp (unit of least precision, a.k.a the spacing between two consecutive FP numbers), which corresponds to a virtual precision of \(t=24\) bits for single-precision and \(t=53\) bits for double-precision.

We applied RR to the CNN models using Verrou [19][22], a tool that implements MCA through dynamic binary instrumentation with Valgrind [23], without needing to modify or recompile the source code. We instrumented the entire executables with RR, additionally using Verrou's custom libmath implementation named Interlibmath to avoid incorrect random perturbations in mathematical functions. We applied RR to FreeSurfer using "fuzzy libmath" [6], a version of the GNU mathematical library instrumented with the Verificarlo [24] compiler following the same principle as Verrou's Interlibmath instrumentation.
### Numerical Uncertainty Metrics

We quantified numerical uncertainty by calculating the number of significant bits across multiple independent RR samples. The number of significant bits is informally defined as the number of bits in common between RR samples for a given FP value. We estimated the number of significant bits using the general non-parametric method described in [25] and implemented in the significant_digits package [26]. Given an RR sample \(X_{i}\) (\(i\leq n\)), this method computes the significance \(S_{i}^{k}\) of the \(k^{th}\) bit in the mantissa of \(X_{i}\) as:

\[S_{i}^{k}=\mathbb{1}_{|Z_{i}|<2^{-k}}\]

where \(Z_{i}=X_{i}-x_{\text{IEEE}}\) and \(x_{\text{IEEE}}\) is the unperturbed result computed with IEEE. The \(k^{th}\) bit in the mantissa is considered significant if the absolute value of \(Z_{i}\) is less than \(2^{-k}\). The number of significant bits across samples, \(\hat{s_{b}}\), is then obtained as the maximal bit index that is significant for all samples:

\[\hat{s_{b}}=\max\left\{k\in\llbracket 1,m\rrbracket\text{ such that }\forall i\in \llbracket 1,n\rrbracket,\ S_{i}^{k}=\mathbb{1}\right\}\]

where \(m\) is the size of the mantissa, i.e., 53 bits for double precision numbers and 24 bits for single precision numbers. A value of 0 significant bits means that \(X\) bears no information while a value of \(m\) means that it has maximal information given the FP format used. The difference between the maximal and achieved values quantifies the information loss resulting from numerical uncertainty.

The number of significant bits is a versatile metric that applies to any program that produces results encoded as FP values. This is, however, not the case of segmentation tools that generally produce categorical variables encoded as integers representing segmentation labels despite the use of intermediate FP operations. Therefore, in order to assess the impact of stochastic rounding in these intermediate FP operations we used the minimum Sorensen-Dice scores computed pairwise across RR samples as uncertainty metric for segmentations. In addition, to have a more local uncertainty metric for segmentation results, we defined an entropy metric at each voxel:

\[E=-\sum_{i=1}^{r}p_{i}\ln p_{i} \tag{1}\]

where \(r\) is the number of segmented regions, and \(p_{i}\) is the probability of region \(i\) at the voxel, computed across \(n\) Random Rounding samples. The entropy is 0 when the voxel is labeled with the same region across all RR samples, and it is maximal when the voxel is labeled differently for each RR sample.

### Non-Linear Registration

SynthMorph[8] is a 3D convolutional U-Net [27] that performs non-linear image registration robustly across MRI contrasts. The encoding section of the U-Net consists of 4 convolutional blocks, while the decoding section consists of 3 blocks
with a skip connection to its respective encoding block. SynthMorph was trained from synthetic label maps created from geometric shapes from which images were generated with various contrasts, deformations, and artifacts. A contrast-invariant loss that measures label overlap was used to optimize model performance. SynthMorph's registration accuracy was shown to outperform state-of-the-art registration methods both within and across contrasts. We applied the SynthMorph "sm-brains" pre-trained model available on GitHub [28] to the linearly registered image produced by FreeSurfer recon-all (output of recon-all's -canorm step), using the MNI305 atlas as reference image. The subject and atlas images were cropped through visual inspection to match the 160 x 192 x 224 dimension required by the model, and intensities were min-max scaled to [0-1]. We applied Verrou directly to the sm-brains model, by instrumenting the entirety of the model and using Verrou's Interlibmath implementation. SynthMorph takes a couple of minutes to run, but when instrumented with Verrou, the runtime increased to a span of 2-3 days.

FreeSurfer "recon-all" [20] is a widely-used neuroimaging pipeline that implements spatial normalization to a brain template, whole-brain segmentation, and cortical surface extraction. The non-linear registration algorithm (mri_ca_register tool) minimizes through gradient descent an error functional that includes an intensity term, a topology constraining one, a metric preservation term, a smoothness term, and a label term [29]. We first ran recon-all on all the subjects with steps --motioncor --talairach --nuitnensitycor --normalization --skullstrip --gcareg --canorm without RR, to obtain the linear registration used as input of SynthMorph. Then we ran the FreeSurfer recon-all step --careg that implements non-linear registration using our FreeSurfer version instrumented with fuzzy libmath. This command typically takes 3 hours to run, but when instrumented with fuzzy libmath, the runtime increased to a span of 6 hours.

### Whole-Brain Segmentation

FastSurfer [12] is a CNN model that performs whole-brain segmentation, cortical surface reconstruction, fast spherical mapping, and cortical thickness analysis. The FastSurfer CNN is inspired from the QuickNAT model [11], which is composed of three 2D fully convolutional neural networks--each associated with a different 2D slice orientation--that each have the same encoder/decoder U-net architecture with skip connections, unpooling layers and dense connections as QuickNAT. The FastSurfer segmentations were shown to surpass state-of-the-art methods, as well as being generalizable to unseen datasets and having better test-retest reliability. We used the pre-trained model from FastSurfer available on GitHub [30] and we applied Verrou directly to this model, in the same way it was applied to SynthMorph. FastSurfer typically takes 30 minutes to run, but when instrumented with Verrou, the runtime increased to a span of 16-17 days.

FreeSurfer recon-all also implements whole-brain segmentation [20], through a maximum a-posteriori estimation of the segmentation based on the non-linear
registration of the subject image to an atlas. Due to time constraints, only the subcortical structures and brain tissues were segmented by FreeSurfer whereas FastSurfer also segmented cortical structures, therefore a mask was applied to FastSurfer's cortical labels to identify them by the super classes "Left/Right Cerebral Cortex". Only regions common to both were further analysed (see list in Fig. 2a). Similar to FreeSurfer recon-all's non-linear registration, RR was applied to FreeSurfer recon-all's whole brain segmentation through Verificarlo's fuzzy libmath. The FreeSurfer recon-all commands --motioncor up to --calabel, the command that specifically performs subcortical segmentation were run. Typically, the segmentation takes around 4 hours to complete, but with FreeSurfer instrumented with Verificarlo, the runtime increased to 10-12 hours.

### Dataset and processing

We used the Consortium for Reliability and Reproducibility (CoRR) dataset [31], a multi-centric, open resource aimed to evaluate test-retest reliability and reproducibility. We randomly selected 35 T1-weighted MRIs from 35 different subjects, one from each CoRR acquisition site, and accessed them through Datalad [32][33]. The selected images included a range of image dimensions, voxel resolutions and data types (Appendix A). We excluded 2 subjects that failed linear registration with FreeSurfer recon-all and a third subject that failed segmentation with FastSurfer. Each pipeline or model was run in Singularity containers over 10 RR samples from which we measured numerical uncertainty. Due to long computing times induced by Verrou instrumentation (\(\approx 17\) days per subject) we were only able to get 4 RR samples for FastSurfer, which we complemented with an IEEE (non-RR) sample conceptually identical to an RR sample.

We processed the data with SynthMorph and FreeSurfer recon-all on the Narval cluster from Ecole de Technologie Superieure (ETS, Montreal), managed by Calcul Quebec and The Digital Alliance of Canada which include AMD Rome 7502, AMD Rome 7532, and AMD Milan 7413 CPUs with 48 to 64 physical cores, 249 GB to 4000 GB of RAM and Linux kernel 3.10. We executed FastSurfer on the slashbin cluster at Concordia University with 8 \(\times\) compute nodes each with an Intel Xeon Gold 6130 CPU, 250 GB of RAM, and Linux kernel 4.18.0-240.1.1.el8_Justre.x86_64.

We used FreeSurfer v7.3.1, SynthMorph v0.2, FastSurfer v2.1.1, Fuzzy v0.9.1, and Singularity/Apptainer v1.1. Verrou v3.21.0 was used for FastSurfer, while Verrou v3.20.0 with a special fix available on GitHub [34] was used for SynthMorph due to compatibility issues between the model and Verrou's Interlibmath. The scripts and Dockerfiles for this experiment can be found on GitHub [35].

## 3 Results

The numerical uncertainty measured for the SynthMorph CNN model was lower than for Freesurfer recon-all (Fig. 1a, as measured both in the resampled images (\(p<10^{-6}\), two-tailed paired t-test) and in the warp fields (\(p<10^{-5}\)) despite
only the libmath libraries in FreeSurfer being instrumented in contrast to the entirety of SynthMorph. The number of significant bits in warp fields was computed as the average number of significant bits across the x,y and z components. On average, out of 24 bits available, the SynthMorph resampled image had 19.56 significant bits while FreeSurfer recon-all's had only 13.43 significant bits; the SynthMorph warp field had 18.55 significant bits while FreeSurfer recon-all's had only 14.12 significant bits. These important differences show a clear superiority of the CNN model compared to FreeSurfer recon-all in terms of numerical uncertainty. Moreover, we also observed a larger variability of the numerical uncertainty across subjects in FreeSurfer recon-all compared to SynthMorph.

The differences in average numerical uncertainty observed between FreeSurfer and SynthMorph were confirmed by visual inspection of the non-linearly registered images and warp fields (Fig. 1b). The numerical uncertainty in registered images was structurally consistent, with higher uncertainty in the gray matter and at the border of the brain than in the white matter, both for SynthMorph and for FreeSurfer recon-all. The numerical uncertainty in warp fields exhibited interesting structural patterns that would benefit from further investigation.

The numerical uncertainty of FastSurfer segmentations was significantly lower than for FreeSurfer recon-all in 31 out of 35 brain regions (Fig. 2a, with very substantial differences in some regions. Here again, a larger variability was observed in FreeSurfer recon-all segmentations than in FastSurfer segmentations. Overall, FastSurfer averages a Sorensen-Dice score of 0.99 across all regions, while FreeSurfer is at 0.92. The differences in Sorensen-Dice scores observed between FreeSurfer recon-all and FastSurfer were confirmed in local entropy maps (Fig. 2b where we visually noted a substantial discrepancy between both methods. For FreeSurfer recon-all, clusters of non-zero entropy values were observed across the brain, whereas for FastSurfer non-zero entropy values were limited to scattered voxels. The entropy maps, in addition to visual inspection, confirm that, despite the relatively high average Sorensen-Dice scores, FreeSurfer recon-all exhibited variability identifying the edges of subcortical structures, while FastSurfer remained certain in its segmentations.

## 4 Conclusion

The numerical uncertainty measured in CNN models SynthMorph and FastSurfer was substantially lower than in FreeSurfer recon-all, amounting to differences in the order of 4 to 6 significant bits in non-linearly registered images, and of up to 0.4 Sorensen-Dice score values in segmentations. We believe that the high numerical uncertainty observed in FreeSurfer recon-all compared to CNN models results from the use of numerical optimization techniques in FreeSurfer recon-all while CNN models only involve low-dimensional convolutions, max-pooling operators, and simple activation functions. The low numerical uncertainty found in CNN models is consistent with previous observations in the very different task of protein function prediction [36]. The numerical uncertainty found in FreeSurfer
recon-all is also consistent with previous observations on FreeSurfer recon-all non-linear registration [6] and segmentation [37].

Our results suggest that neuroimaging CNN models are significantly more robust to small numerical perturbations than traditional image processing approaches. Therefore, we expect CNN results to be more reproducible across execution environments than traditional image processing approaches, implying better portability across software and hardware systems.

Our results report on the numerical uncertainty resulting from CNN _inference_, which is a relevant proxy for the uncertainty experienced by model end-users across different execution environments. However, the numerical uncertainty resulting from CNN _training_ was not measured in our experiments. We speculate that some of the numerical uncertainty observed in FreeSurfer recon-all results is intrinsic to the problems of subject-to-template non-linear registration and whole-brain segmentation, and should therefore manifest during CNN training. Mathematically, training CNN models involves numerical optimization in high-dimensional spaces, which we expect to be less numerically stable than CNN inference, and comparably stable to FreeSurfer recon-all. Should this assumption be accurate, the numerical uncertainty of predictions made by a sample of CNN models trained with Random Rounding should be substantial,

Figure 1: Numerical uncertainty measured in the non-linearly registered images and warp fields produced by FreeSurfer recon-all and the SynthMorph CNN model.
Figure 2: Numerical uncertainty measured in the segmentations produced by FreeSurfer recon-all and the FastSurfer CNN model.
which we plan to leverage in our future work by building efficient ensemble models capturing the numerical variability associated with non-linear registration or segmentation, possibly resulting in improved predictions.

## 5 Acknowledgements

Computations were made on the Narval and Beluga supercomputers from Ecole de Technologie Superieure (ETS, Montreal), managed by Calcul Quebec and The Digital Alliance of Canada. The operation of these supercomputers are funded by the Canada Foundation for Innovation (CFI), le Ministere de l'Economie, des Sciences et de l'Innovation du Quebec (MESI) and le Fonds de recherche du Quebec - Nature et technologies (FRQ-NT).
