<table><caption>Table 13: Average detection accuracy of unseen attacks after training on FGM.</caption><thead><tr><th rowspan="2">Network/Data</th><th rowspan="2">Detector</th><th colspan="7">Attack</th></tr><tr><th>APGD</th><th>BIM</th><th>AA</th><th>DF</th><th>CW</th><th>HSJ</th><th>BA</th></tr></thead><tbody><tr><td rowspan="3">ResNet110 CIFAR10</td><td>TR</td><td>89.32</td><td><strong>96.02</strong></td><td><strong>85.10</strong></td><td><strong>91.02</strong></td><td><strong>93.18</strong></td><td><strong>93.00</strong></td><td><strong>90.92</strong></td></tr><tr><td>MH</td><td>77.34</td><td>77.24</td><td>72.12</td><td>80.12</td><td>79.92</td><td>79.70</td><td>79.32</td></tr><tr><td>NS</td><td><strong>92.08</strong></td><td>93.88</td><td>51.82</td><td>51.40</td><td>50.84</td><td>52.12</td><td>59.88</td></tr><tr><td rowspan="3">ResNeXt50 CIFAR100</td><td>TR</td><td>91.94</td><td><strong>95.02</strong></td><td><strong>73.32</strong></td><td><strong>85.16</strong></td><td><strong>78.18</strong></td><td><strong>85.04</strong></td><td><strong>92.14</strong></td></tr><tr><td>MH</td><td>90.86</td><td>93.16</td><td>73.08</td><td>82.72</td><td>76.44</td><td>82.82</td><td>84.46</td></tr><tr><td>NS</td><td><strong>92.16</strong></td><td>93.88</td><td>51.32</td><td>51.62</td><td>51.02</td><td>52.04</td><td>57.90</td></tr><tr><td rowspan="3">WideResNet TinyImageNet</td><td>TR</td><td>93.26</td><td><strong>94.66</strong></td><td><strong>77.04</strong></td><td><strong>90.62</strong></td><td><strong>91.42</strong></td><td></td><td></td></tr><tr><td>MH</td><td>76.96</td><td>77.02</td><td>60.36</td><td>73.18</td><td>75.52</td><td></td><td></td></tr><tr><td>NS</td><td><strong>94.06</strong></td><td>94.62</td><td>65.60</td><td>72.82</td><td>71.96</td><td></td><td></td></tr></tbody></table>

On our detector, the detection rate of successful adversarial samples from unseen attacks remains higher than 90% in most cases (Table 39 in Appendix H.4) and the FPR is always lower than 10% (Table 41 in Appendix H.5). The AUROC is in Table 43 in Appendix H.6. Our detector almost always outperforms the other detectors on all these metrics. On the attacks that are not variants of FGM, the NS detector almost never predicts that an input is adversarial.

However, this experiment shows that our approach has some limitations. We see in Tables 35, 36 and 37 in Appendix H.3 that LAP training does not improve detection accuracy as much anymore in this setting, and reduces it in some cases. It improves it for the Mahalanobis detector on all attacks on the ResNet110-CIFAR10 and WideResNet-TinyImageNet experiments (by up to 10 percentage points), but not on most attacks on the ResNeXt50-CIFAR100 experiment. It often reduces detection generalization on our detector. But it always does better than RCE training in Tables 35 and 36 in Appendix H.3.

We claim this is because these methods reduce the variance of features extracted during training on the seen attack, thus harming generalization to unseen attacks. This explains why detection of PGD and BIM, which are variants of FGM, still often improves. Fixing this could be an area of future research. That LAP training improves detection of unseen OOD samples in some cases in the out-of-distribution experiment in Tables 14 and 15 in Section 5.4.4 indicates that this problem is not general.

### 5.4.4 Detection of Out-of-Distribution Samples

Since our analysis applies to all out-of-distribution (OOD) samples, we test detection of OOD samples in a similar setting to the Mahalanobis paper [Lee et al., 2018]. We use the same ResNet110 and ResNeXt50 models trained on CIFAR10 and CIFAR100 respectively. Since the detectors need to be trained, we are in the OOD setting where we have a first dataset for training the network (CIFAR10 in Table 14 and CIFAR100 in Table 15) and a second dataset from another distribution that is not the test OOD distribution to train the detector
on. This could be another dataset (CIFAR100 in the first column of Table 14), some images found in the wild, or a perturbation of our dataset that we generate using an adversarial attack (CW on CIFAR10 in the second column of Table 14, and AA and CW on CIFAR100 in Table 15). Detectors can then be used by training them to distinguish between these first two datasets, and then testing them on distinguishing between the first dataset and a third unseen dataset (SVHN [Netzer et al., 2011] in both tables). Results are in Tables 14 and 15.

<table><caption>Table 14: Average OOD detection accuracy and standard deviation over 5 runs using ResNet110 trained on CIFAR10.</caption><thead><tr><th rowspan="2">Detector</th><th colspan="2">OOD Experiment 1</th><th colspan="2">OOD Experiment 2</th></tr><tr><th>CIFAR100 (seen)</th><th>SVHN (unseen)</th><th>CW-CIFAR10 (seen)</th><th>SVHN (unseen)</th></tr></thead><tbody><tr><td>VAN TR</td><td>98.30 ± 0.46</td><td>97.46 ± 0.49</td><td><strong>97.42</strong> ± 0.57</td><td><strong>91.38</strong> ± 0.95</td></tr><tr><td>RCE TR</td><td><strong>98.42</strong> ± 0.40</td><td>98.20 ± 0.39</td><td>91.54 ± 6.06</td><td>77.58 ± 6.72</td></tr><tr><td>LAP TR</td><td>98.30 ± 0.22</td><td><strong>98.50</strong> ± 0.47</td><td>97.28 ± 0.62</td><td>85.46 ± 2.64</td></tr><tr><td>VAN MH</td><td>86.88 ± 1.52</td><td>91.28 ± 0.92</td><td>81.80 ± 1.96</td><td>83.76 ± 1.13</td></tr><tr><td>RCE MH</td><td>94.82 ± 0.45</td><td>92.16 ± 0.57</td><td>76.74 ± 2.75</td><td>54.24 ± 3.46</td></tr><tr><td>LAP MH</td><td>94.84 ± 0.41</td><td>90.46 ± 1.45</td><td>89.68 ± 0.65</td><td>76.72 ± 1.73</td></tr></tbody></table>

<table><caption>Table 15: Average OOD detection accuracy and standard deviation over 5 runs using ResNeXt50 trained on CIFAR100.</caption><thead><tr><th rowspan="2">Detector</th><th colspan="2">OOD Experiment 1</th><th colspan="2">OOD Experiment 2</th></tr><tr><th>AA-CIFAR100 (seen)</th><th>SVHN (unseen)</th><th>CW-CIFAR100 (seen)</th><th>SVHN (unseen)</th></tr></thead><tbody><tr><td>VAN TR</td><td>84.48 ± 0.59</td><td>75.32 ± 0.62</td><td>95.82 ± 0.67</td><td><strong>92.92</strong> ± 1.36</td></tr><tr><td>RCE TR</td><td>50.04 ± 0.07</td><td>55.44 ± 5.76</td><td>76.48 ± 0.75</td><td>75.66 ± 0.68</td></tr><tr><td>LAP TR</td><td><strong>87.10</strong> ± 0.12</td><td>72.98 ± 2.72</td><td><strong>95.94</strong> ± 0.57</td><td>85.94 ± 2.88</td></tr><tr><td>VAN MH</td><td>83.44 ± 0.48</td><td>78.82 ± 0.48</td><td>94.96 ± 0.81</td><td>85.10 ± 1.50</td></tr><tr><td>RCE MH</td><td>50.04 ± 0.07</td><td>58.74 ± 2.36</td><td>76.20 ± 0.72</td><td>72.20 ± 1.47</td></tr><tr><td>LAP MH</td><td>86.04 ± 0.31</td><td><strong>86.84</strong> ± 0.68</td><td>94.82 ± 0.34</td><td>88.92 ± 1.26</td></tr></tbody></table>

Our detector performs very well and better than the MH detector in three of the four experiments, and in the fourth case, the MH detector benefits from LAP training by 8 percentage points (first experiment in Table 15). Without any extra data available, using the CW adversarial attack allows to detect OOD samples from an unseen distribution with more than 90% accuracy and an FPR of less than 10% at a fixed TPR of 95%. The choice of the attack is also important, as CW allows for much better detection of unseen samples from SVHN than AA.

The AUROC and the false positive rate at a fixed true positive rate of 95% are in Appendix H.7, in Tables 44 and 45 for ResNet100, and in Tables 46 and 47 for ResNeXt50. They show the same advantage for our detector and regularization over the Mahalanobis detector and RCE training respectively.
### 5.4.5 Attacking the Detector

We consider here the case where the attacker also attacks the detector (adaptive attacks). We try two such attacks on the TR and MH detectors on ResNet110 trained on CIFAR10. Both attacks are white-box with respect to the network. The first is black-box with respect to the detector. It only knows if an adversarial sample has been detected or not. The second has some knowledge about the detector. It knows what features it uses and can attack it directly to find adversarial features. We test these attacks by looking at the percentage of detected successful adversarial samples that they turn into undetected successful adversarial samples that fool both the network and the detector.

The first attack proceeds as follows. A strong white-box attack (CW) is used on the network on image $x$ that has label $y$. If it finds a successful adversarial image $\tilde{x}$ that fools the network into predicting $\tilde{y} \neq y$ but is detected by the detector, the attacker will attempt to modify this image $\tilde{x}$ so that the network and the detector are both fooled. For this, the image $\tilde{x}$ is used as the initialization for an attack (HSJ with a budget of 50 iterations and 10000 evaluations) on a black-box Network-Detector system. The attacker considers that the Network-Detector behaves as follows: it outputs the class prediction of the network if the detector does not detect an attack and outputs an additional 'detected' class if the detector detects an attack. The attacker attacks this Network-Detector on image $\tilde{x}$ targeting the $\tilde{y}$ label. This way the network makes a mistake and the 'detected' class is avoided. On the vanilla ResNet110, this attack turns 16.5% of 1700 detected successful adversarial samples $\tilde{x}$ into undetected successful adversarial samples on our detector, compared to 25.7% on the Mahalanobis detector. These percentages are lower on the LAP-ResNet110 as they drop to 6.8% on our detector and 12.9% on the Mahalanobis detector. This shows that LAP training improves the robustness of both adversarial detectors to being attacked themselves, and that the Transport detector is more robust than the MH detector.

The second attack is very similar to the adaptive attack used in [Carlini and Wagner, 2017a] to break the Kernel Density detector of [Feinman et al., 2017]. It proceeds as follows. A strong white-box attack (CW) is used on the network on image $x$ that has label $y$. If it finds a successful adversarial image $\tilde{x}$ that fools the network but is detected by the detector, the detection features $\tilde{z}$ that $\tilde{x}$ generates when run through the network are used as the initialization for a black-box attack (HSJ with a budget of 50 iterations and 10000 evaluations) on the detector. If successful adversarial detection features $z^*$ that fool the detector are found, the attacker has to find an adversarial perturbation of $x$ that still fools the network and that generates these features $z^*$ (or close features that also fool the detector) when run through the network. We do this as in [Carlini and Wagner, 2017a] by solving the following optimization problem:

$$ \min_{x^*} -L(N(x^*), y) + c_1 \|D(x^*) - z^*\| + c_2 \|x^* - x\| \quad (42) $$

where $L$ is the cross-entropy loss, $N$ is the network, and $D$ is the (differentiable)
function that returns the detection features of its input. This optimization problem is differentiable and we try differentiable optimization algorithms such as BFGS and NR to solve it. The initial detected successful adversarial image $\tilde{x}$ is used as initialization as in [Carlini and Wagner, 2017a]. This attack turns 14% of detected successful adversarial samples $\tilde{x}$ into undetected successful adversarial samples on our detector on the LAP-ResNet110.

Given that initial detection rates of successful adversarial samples are very often close to 100% (see Tables 38 and 39 in Appendix H.4), this shows that adaptive attacks do not (at least not easily) circumvent the detector, as detection rates drop to 85% at worst.

The second attack is stronger than the first one, but it can probably still be improved by using a white-box attack that is specific to random forests for attacking the detector such as [Kantchelian et al., 2016] or [Zhang et al., 2020], or a different loss than cross-entropy such as the one used in the CW attack. The choice of the differentiable term representing the detector (here $||D(x^*) - z^*||$ in (42)) is also important. However, the difficulty of combining the attack on the network with that on the detector remains. It is the non-differentiability of the random forest that forces either this separate treatment of network and detector then the use of a proxy differentiable term for the detector in (42) to combine both, or the use of a black-box method as in the first attack.

Also, we did not consider here the ensemble of the class-conditional detector and the general detector, which is the best performing version of the detector (see Section 5.4.2), and should be more robust to adaptive attacks, as the attacker will have to fool two random forest detectors at once and target a particular label, constraining further the optimization problem he solves. This is similar to [Eniser et al., 2020] where a random detector is selected at test time from a pool of trained detectors to render adaptive attacks on the detector more difficult.

### 5.4.6 Time Comparison

With a ResNeXt50 on CIFAR100 and a Tesla V100 GPU, it takes our method (including the time to generate FGM attacks) 66 seconds to extract its features from both the clean and the adversarial samples, while it takes the Mahalanobis method 110 seconds. Mahalanobis also extracts some statistics from the training set prior to adversarial training, which takes an additional 35 seconds. Our feature vector is of size 32, compared to 5 for the Mahalanobis detector. So our random forest takes only 4 more seconds to train than the Mahalanobis one (7 vs 3 seconds). Computation of the features our detector uses (norms and cosines) is in $O(MD)$, where $M$ is the number of residual blocks and $D$ is the largest embedding dimension inside the network.

## 5.5 Discussion and Conclusion

To conclude, we proposed a method for detecting abnormal samples, based on the dynamical viewpoint of neural networks. The method examines the discrete
vector field moving the inputs to distinguish clean and adversarial samples. The detector requires minimal processing and computation to extract the features it needs for detection and achieves state-of-the-art detection performances on 3 different architectures and datasets, on seen and unseen attacks, and also on out-of-distribution detection.

The low dimension of the features we extract (2 per residual block) allows to keep the features from all blocks as opposed to other detection methods that only experimentally select which layers to use. This could explain the superior performance of our detector, as we closely monitor the evolution of the input's embeddings inside the network.

We have also used the transport regularization that was previously shown to improve test accuracy to improve the detection accuracy of adversarial detectors. This regularization provably makes the embeddings closer to each other on the support of the input distribution, thus making them more distinguishable from those of abnormal samples. The main difference with other regularizations that make the network more Lipschitz by regularizing its weights directly is that we only endue the network with additional regularity on the support of the data distribution and not on the entire space.

A possible simple extension of our detection method would be to add cosines to other orthogonal vectors as detection features. A trade-off between the accuracy of the detector and the size and computation time of the features would then have to be taken into account.
<figure><img src="image_29.png" /></figure>
# 6 Conclusion

## 6.1 Summary

Residual connections are among the workhorses that make neural networks successful. We have expanded on interpretations of their functioning, namely, the iterative refinement bias (i.e. the bias to remain close to the identity function) and the discrete dynamical system interpretation. We use optimal transport theory to theoretically analyze their functioning. In optimal transport terms, these two interpretations are unified to show that residual connections induce a desirable bias towards networks with low kinetic energy and transport cost. This analysis allowed us to propose improvements to the accuracy and stability of residual networks and applications to module-wise training that make it competitive with end-to-end training with only a fraction of the memory requirements and to detection of adversarial and out-of-distribution samples.

In a first work [Karkar et al., 2020] detailed in Section 3, we empirically motivated and formulated the following learning problem: among all the networks that solve the task, find the one that transforms the data with the lowest cost. This amounts to regularizing the transport cost of the network. We proved existence and regularity of models trained this way, studied their behavior in low-dimensional settings and shown their efficiency on standard classification tasks. Training is stabilized in an adaptive fashion without being slowed down.

In a second work [Karkar et al., 2023a] detailed in Section 4, we used this regularization for module-wise training of neural networks without end-to-end backpropagation, including on Transformers, saving up to 60% memory usage. We proved that this makes every module apply a proximal step in the Wasserstein space to maximize the separability of the data.

In a third work [Karkar et al., 2023b] detailed in Section 5, we proposed a detector of adversarial and out-of-distribution samples based on the dynamic viewpoint of neural networks. We also used our detector successfully on detection of unseen out-of-distribution samples. The transport regularization provably makes the activations of out-of-distribution samples more distinguishable from those of clean samples, thus improving the accuracy of adversarial detectors.

## 6.2 Future Directions

Besides the direct possible extensions mentioned in the conclusions to Sections 3 to 5, other avenues of research follow from our analysis and point of view.

### 6.2.1 Adaptive Depth During Training

#### 6.2.1.1 Introduction

The dynamical point of view of neural networks, initiated by the analogy between residual networks and the Euler scheme, has led to new architectures [Lu et al.,
2018, Haber et al., 2019]. While these new architectures often imitate numerical solvers of higher order than the Euler scheme through skip connections, adaptive solvers have only yet been explored in a few papers. We propose to use them to adjust the step size and the depth of the network during training, possibly reducing training cost. Given that adaptive methods are made up of multiple imbedded schemes, this also allows to ensemble different pre-classification forward passes at inference, without an additional training cost.

### 6.2.1.2 Method

This is currently a work in progress with positive results on the CIFAR100 dataset. We build a neural network that imitates the structure of an adaptive Runge-Kutta scheme (see Appendices A.2 and A.3). The block we use contains 4 convolutional functions denoted $f_i$ and does the following computations for an input $x$ with a step size $h$:

$$ \left\{ \begin{array}{l} p_1 = f_1(x) \\ y_{\text{euler}} = x + hp_1 \\ p_2 = f_2(x + \frac{1}{2}hp_1) \\ y_{\text{midpoint}} = x + hp_2 \\ p_3 = f_3(x + \frac{3}{4}hp_2) \\ y_{\text{order3}} = x + \frac{2}{9}hp_1 + \frac{1}{3}hp_2 + \frac{4}{9}hp_3 \\ p_4 = f_4(y_{\text{order3}}) \\ y_{\text{order2}} = x + \frac{7}{24}hp_1 + \frac{1}{4}hp_{n,2} + \frac{1}{3}hp_3 + \frac{1}{8}hp_4 \end{array} \right. \qquad (43) $$

Depending on the numerical scheme we want to use, the block forwards one of the outputs $y$ to the next block. Only layers $f_1$ and $f_2$ will be trained for example if we want to use the midpoint scheme. This allows to compare the efficiency of the different schemes for solving classification tasks, i.e. is it better to train 1 block of order 3 or 3 blocks of order 1 (since the order 3 method contains 3 functions $f_i$)?

If all the functions $f_i$ are learned, and since all the values $y$ approximate the same differential equation, it makes sense to average them after the last block before the classification layer as an ensembling method.

The depth is adaptively adjusted during training as follows. First, a tolerance hyperparameter $\epsilon > 0$ is fixed. Then, on the validation set, every $S$ epochs, the average error $err = \|y_{\text{order3}} - y_{\text{order2}}\|_2$ for each block is computed. If it is larger than $\epsilon$ then the step size is halved and the block is duplicated, i.e. a block with the same weights is added right after it. If the error is smaller than $\epsilon/3$, then the step size is doubled and the following block is deleted.

To summarize, the contributions we aim to verify on larger datasets are:

- Using neural blocks of higher order is more efficient than using normal Euler residual blocks, meaning that even though the blocks are larger,
fewer are required to reach higher accuracy with less memory usage.

- Since adaptive scheme are made-up of interwoven schemes of different order, different forward passes representing the schemes of different order can be ensembled at test time, through averaging the final embeddings before classification, to improve accuracy.

- Adaptively adjusting the step-sizes and depth of the network allows to reduce, since the first training epochs, the size of a large network, while still obtaining the same performance at the end of training.

- This architecture also benefits from LAP training as in Section 3.

### 6.2.1.3 Related Work

[Yang et al., 2020a] introduce a time-step controller that is an LSTM [Hochreiter and Schmidhuber, 1997] which is dependent on the convolution weights of the current and previous layers and that is jointly optimized with the network. In contrast, our time-step is adjusted as in adaptive numerical schemes, therefore requiring no extra training, and also adjusts the depth of the network during training. [Chang et al., 2018b] use the analogy between ResNets and numerical schemes for differential equations to accelerate training by starting the training with a shallow network with a large step-size and then doubling the depth and halving the step-size during training. This is not an adaptive method to adjust the depth and step-size, as it is not dependent on the stability of the approximated differential equation. [Zhu et al., 2022] simply implement an architecture that imitates a Runge-Kutta scheme. However, the coefficients are not fixed to the theoretical values of the scheme and are instead learned inside the convolutions. The relations between them necessary to the theoretical validity of the scheme are therefore not enforced. There is also no modification of the step-size and depth. Likewise, [Lu et al., 2018] use an architecture that imitates the linear multi-step numerical method and [Haber et al., 2019] imitate the implicit-explicit method. Modifying the step-size has shown its benefits in [Zhang et al., 2019, Zhang and Wynter, 2018], where fixing it to a small value around 0.1 makes training more stable.

Somewhat similar to our proposed method is the adaptive computation time method applied to residual networks [Figurnov et al., 2017], inspired by a similar idea applied to recurrent neural networks [Graves, 2017, Fojo et al., 2018, Jernite et al., 2017]. In these methods, secondary neural networks are trained to choose (minimize) the number of times the state transition is applied for an input step in a recurrent network, and the number of residual blocks applied for an input in a residual network. The goal in [Figurnov et al., 2017] is faster computation time and not the adjustment during training of the size of the network to the dataset at hand, as the number of residual blocks applied is chosen for each input (or even patch of the input), so the entire network still has to be trained.
## 6.2.2 Similarity of Neural Network Representations

We propose to use the study of the dynamics inside a neural network, especially as used in Section 5 to detect adversarial and out-of-distribution samples, to consider the problem of measuring similarity between the representations learned by neural networks. This area of research [Raghu et al., 2017, Kornblith et al., 2019, Laakso and Cottrell, 2000, Li et al., 2015, Morcos et al., 2018, Wang et al., 2018] tries to answer questions such as:

- Do different architectures learn similar representations when trained on the same dataset?

- Does the same architecture tend to learn similar representations when trained on similar datasets ?

- How much does changing the initialization change the learned representations?

- Can we establish correspondences between layers of different networks?

The available methods usually propose similarity measures (e.g. SVCCA [Raghu et al., 2017] and CKA [Kornblith et al., 2019]) between activations that have certain desirable invariances such as invariance to orthogonal transformations.

Optimal transport theory was used in [Singh and Jaggi, 2020] to establish correspondences between the neurons of two networks, which is then useful for fusing these networks into a network that keeps as much of the capabilities of each network as possible. However, [Singh and Jaggi, 2020] do this only layer by layer and not across layers, which could be an extension of their method.

## 6.2.3 Optimal Transport between Different Dimensions

Another avenue of investigation is to use optimal transport between distributions that are supported in spaces of different dimensions [McCann and Pass, 2019, Pass, 2010] to regularize pooling and embedding layers that change the dimension of their input. This could extend the work done here to different architectures, and could be particularly useful to regularizing generative network where we would like to start generating in lower resolutions in the earlier layers before upsampling to use less resources and computation time.

## 6.2.4 Simulating Continuity Partial Differential Equations

In Appendix C.3, the minimizing movement scheme defined in Definition 3 in Section 4.3.2 as

$$ \rho_{k+1}^{\tau} \in \arg \min_{\rho \in \mathcal{P}(\Omega)} \mathcal{Z}(\rho) + \frac{1}{2\tau} W_2^2(\rho, \rho_k^{\tau}) $$
for $\mathcal{Z}: W_2(\Omega) \to \mathbb{R}$, is shown to converge, under some conditions on $\mathcal{Z}$, as $\tau$ goes to 0, to a curve $\tilde{\rho}_t$ in $W_2(\Omega)$ that solves the continuity equation

$$ \partial_t \rho_t(x) - \nabla \cdot (\rho_t(x) \nabla(\frac{\delta \mathcal{Z}}{\delta \rho}(\rho_t))(x)) = 0 $$

where $\frac{\delta \mathcal{Z}}{\delta \rho}$ is the first variation of $\mathcal{Z}$ defined in Definition 20 in Appendix C.3.

For example, if $\mathcal{Z}(\rho) = \int g(\rho(x)) dx$ with $g(x) = x \log x$, we obtain the heat equation $\partial_t \rho_t - \Delta \rho = 0$. And if $\mathcal{Z}(\rho) = \int g(\rho(x)) dx + \int V d\rho$, we obtain the Fokker-Planck equation $\partial_t \rho_t - \Delta \rho - \nabla \cdot (\rho \nabla V) = 0$.

This suggests that solving the problems in the minimizing movement scheme by approximating optimal transport maps through residual neural networks could be a way to discretize and numerically solve such continuity partial differential equations.

## 6.3 Conclusion

This principle of least action is central in physics. It can be found in classical and relativistic mechanics, thermodynamics and quantum mechanics [García-Morales et al., 2008, Feynman, 1942, Gray, 2009]. It states that the dynamical trajectory of a system is one that makes a certain action locally stationary. Since essentially all the laws of physics can be derived through some version of this principle, this suggests that nature operates in an optimal and efficient way. The principle of least action was therefore extended to derive laws of mathematical economy [Kombarov, 2021] and of evolution and natural selection [Kaila and Annila, 2008]. Neural networks that respect conservation laws have been proposed to approximate physical systems where conservation of energy is important [Greydanus et al., 2019].

Here, we have found that this principle is also central to finding good solutions to machine learning problems outside of physics, and to doing so while minimizing the resources and computations required. In particular, it explains why many techniques such residual connections, batch normalizations and orthogonal initialization have been so successful in deep learning, and opens the door to new considerations when designing novel neural networks and their training methods.
<figure><img src="image_30.png" /></figure>
References

[Aldahdooh et al., 2022] Aldahdooh, A., Hamidouche, W., Fezza, S. A., and Déforges, O. (2022). Adversarial example detection for dnn models: a review and experimental comparison. *Artificial Intelligence Review*.

[Alemany and Pissinou, 2022] Alemany, S. and Pissinou, N. (2022). The dilemma between data transformations and adversarial robustness for time series application systems. In *Proceedings of the Workshop on Artificial Intelligence Safety 2022 (SafeAI 2022) co-located with the Thirty-Sixth AAAI Conference on Artificial Intelligence (AAAI2022), Virtual, February, 2022*, volume 3087 of *CEUR Workshop Proceedings*. CEUR-WS.org.

[Alvarez-Melis and Fusi, 2021] Alvarez-Melis, D. and Fusi, N. (2021). Dataset dynamics via gradient flows in probability space. *Proceedings of the International Conference on Machine Learning (ICML 2021)*.

[Amari, 1967] Amari, S. (1967). A theory of adaptive pattern classifiers. *IEEE Transactions on Electronic Computers*, EC-16(3):299-307.

[Ambrosio et al., 2005] Ambrosio, L., Gigli, N., and Savare, G. (2005). *Gradient Flows in Metric Spaces and in the Space of Probability Measures*. Birkhäuser Basel.

[Amodei et al., 2016] Amodei, D., Olah, C., Steinhardt, J., Christiano, P., Schulman, J., and Mané, D. (2016). Concrete problems in ai safety. *arXiv*.

[Amsaleg et al., 2015] Amsaleg, L., Chelly, O., Furon, T., Girard, S., Houle, M. E., Kawarabayashi, K.-i., and Nett, M. (2015). Estimating local intrinsic dimensionality. In *Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, KDD '15, page 29-38, New York, NY, USA. Association for Computing Machinery.

[Andriushchenko et al., 2020] Andriushchenko, M., Croce, F., Flammarion, N., and Hein, M. (2020). Square attack: a query-efficient black-box adversarial attack via random search. In *Computer Vision - ECCV 2020*. Springer.

[Ansari et al., 2021] Ansari, A. F., Ang, M. L., and Soh, H. (2021). Refining deep generative models via discriminator gradient flow. In *International Conference on Learning Representations (ICLR 2021)*.

[Arbel et al., 2019] Arbel, M., Korba, A., SALIM, A., and Gretton, A. (2019). Maximum mean discrepancy gradient flow. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Arora et al., 2018] Arora, R., Basu, A., Mianjy, P., and Mukherjee, A. (2018). Understanding deep neural networks with rectified linear units. In *International Conference on Learning Representations (ICLR 2018)*.
[Arora et al., 2014] Arora, S., Bhaskara, A., Ge, R., and Ma, T. (2014). Prov- able bounds for learning some deep representations. In *Proceedings of the International Conference on Machine Learning (ICML 2014)*.

[Ascher and Petzold, 1998] Ascher, U. and Petzold, L. (1998). *Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations*. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104).

[Ascher et al., 1997] Ascher, U. M., Ruuth, S. J., and Spiteri, R. J. (1997). Implicit-explicit runge-kutta methods for time-dependent partial differential equations. *Applied Numerical Mathematics*, 25(2):151–167. Special Issue on Time Integration.

[Ascher et al., 1995] Ascher, U. M., Ruuth, S. J., and Wetton, B. T. R. (1995). Implicit-explicit methods for time-dependent partial differential equations. *SIAM Journal on Numerical Analysis*, 32(3):797–823.

[Athans and Falb, 2007] Athans, M. and Falb, P. (2007). *Optimal Control: An Introduction to the Theory and Its Applications*. Dover Books on Engineering. Dover Publications.

[Ba et al., 2016] Ba, J. L., Kiros, J. R., and Hinton, G. E. (2016). Layer normalization. *arXiv*.

[Bach, 2017] Bach, F. (2017). Breaking the curse of dimensionality with convex neural networks. *Journal of Machine Learning Research (JMLR 18)*, 18(19):1–53.

[Baluzzi et al., 2017] Baluzzi, D., Frean, M., Leary, L., Lewis, J., Ma, K. W.-D., and McWilliams, B. (2017). The shattered gradients problem: If resnets are the answer, then what is the question? In *Proceedings of the International Conference on Machine Learning (ICML 2017)*.

[Baluja and Fischer, 2017] Baluja, S. and Fischer, I. (2017). Adversarial transformation networks: Learning to generate adversarial examples. *arXiv*.

[Baydin et al. 2022] Baydin, A. G., Pearlmutter, B. A., Syme, D., Wood, F., and Torr, P. (2022). Gradients without backpropagation. *arXiv*.

[Belilovsky et al., 2019] Belilovsky, E., Eickenberg, M., and Oyallon, E. (2019). Greedy layerwise learning can scale to imagenet. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[Belilovsky et al., 2020] Belilovsky, E., Eickenberg, M., and Oyallon, E. (2020). Decoupled greedy learning of cnns. In *Proceedings of the International Conference on Machine Learning (ICML 2020)*.

[Belkin et al., 2018] Belkin, M., Ma, S., and Mandal, S. (2018). To understand deep learning we need to understand kernel learning. In *Proceedings of the International Conference on Machine Learning (ICML 2018)*.
[Belkin et al., 2019] Belkin, M., Ma, S., and Mandal, S. (2019). Reconciling modern machine-learning practice and the classical bias-variance trade-off. In *Proceedings of the National Academy of Sciences (PNAS)*, volume 116, pages 15849–15854.

[Benamou and Brenier, 2000] Benamou, J. and Brenier, Y. (2000). A computational fluid mechanics solution to the monge-kantorovich mass transfer problem. *Numerische Mathematik*.

[Bengio et al., 2006] Bengio, Y., Lamblin, P., Popovici, D., and Larochelle, H. (2006). Greedy layer-wise training of deep networks. In Schölkopf, B., Platt, J., and Hoffman, T., editors, *Advances in Neural Information Processing Systems (NeurIPS 2006)*, volume 19. MIT Press.

[Bernd Illing, 2020] Bernd Illing, Wulfram Gerstner, G. B. (2020). Towards truly local gradients with clapp: Contrastive, local and predictive plasticity. *arXiv*.

[Bertasius et al., 2021] Bertasius, G., Wang, H., and Torresani, L. (2021). Is space-time attention all you need for video understanding? In Meila, M. and Zhang, T., editors, *Proceedings of the 38th International Conference on Machine Learning (ICML 2021)*, volume 139 of *Proceedings of Machine Learning Research*, pages 813–824. PMLR.

[Bertsekas, 2016] Bertsekas, D. (2016). *Nonlinear Programming*. Athena scientific optimization and computation series. Athena Scientific.

[Bhowmick et al., 2021] Bhowmick, A., D'Souza, M., and Raghavan, G. S. (2021). Lipbab: Computing exact lipschitz constant of relu networks. *arXiv*, abs/2105.05495.

[Bolley, 2008] Bolley, F. (2008). Separability and completeness for the wasserstein distance. In *Séminaire de Probabilités XLI*. Springer, Berlin, Heidelberg.

[Bottou, 2012] Bottou, L. (2012). *Stochastic Gradient Descent Tricks*, pages 421–436. Springer Berlin Heidelberg, Berlin, Heidelberg.

[Bottou et al., 2018] Bottou, L., Curtis, F. E., and Nocedal, J. (2018). Optimization methods for large-scale machine learning. *SIAM Review*, 60(2):223–311.

[Brendel et al., 2018] Brendel, W., Rauber, J., and Bethge, M. (2018). Decision-based adversarial attacks: Reliable attacks against black-box machine learning models. In *International Conference on Learning Representations (ICLR 2018)*.

[Carlini and Wagner, 2017a] Carlini, N. and Wagner, D. (2017a). *Adversarial Examples Are Not Easily Detected: Bypassing Ten Detection Methods*, page 3–14. Association for Computing Machinery, New York, NY, USA.

[Carlini and Wagner, 2017b] Carlini, N. and Wagner, D. (2017b). Towards evaluating the robustness of neural networks. In *2017 IEEE Symposium*
on Security and Privacy (SP), pages 39–57, Los Alamitos, CA, USA. IEEE Computer Society.

[Carrara et al., 2019] Carrara, F., Becarelli, R., Caldelli, R., Falchi, F., and Amato, G. (2019). Adversarial examples detection in features distance spaces. In Leal-Taixé, L. and Roth, S., editors, *Workshop at European Conference on Computer Vision (ECCV 2018)*, pages 313–327, Cham. Springer International Publishing.

[Chang et al., 2018a] Chang, B., Meng, L., Haber, E., Ruthotto, L., Begert, D., and Holtham, E. (2018a). Reversible architectures for arbitrarily deep residual neural networks. In *The AAAI Conference on Artificial Intelligence (AAAI-18)*.

[Chang et al., 2018b] Chang, B., Meng, L., Haber, E., Tung, F., and Begert, D. (2018b). Multi-level residual networks from dynamical systems view. In *International Conference on Learning Representations (ICLR 2018)*.

[Chattopadhyay et al., 2019] Chattopadhyay, N., Chattopadhyay, A., Gupta, S. S., and Kasper, M. (2019). Curse of dimensionality in adversarial examples. In *International Joint Conference on Neural Networks (IJCNN 2019)*.

[Chen et al., 2020a] Chen, J., Jordan, M. I., and Wainwright, M. J. (2020a). Hopskipjumpattack: A query-efficient decision-based attack. In *2020 IEEE Symposium on Security and Privacy*, pages 1277–1294. IEEE.

[Chen et al., 2017a] Chen, L.-C., Papandreou, G., Kokkinos, I., Murphy, K., and Yuille, A. L. (2017a). Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs. *IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI 2017)*, 40.

[Chen et al., 2017b] Chen, P.-Y., Zhang, H., Sharma, Y., Yi, J., and Hsieh, C.-J. (2017b). Zoo: Zeroth order optimization based black-box attacks to deep neural networks without training substitute models. In *Proceedings of the 10th ACM Workshop on Artificial Intelligence and Security*, AISec ’17, page 15–26, New York, NY, USA. Association for Computing Machinery.

[Chen et al., 2018] Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. (2018). Neural ordinary differential equations. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Chen et al., 2020b] Chen, T., Kornblith, S., Norouzi, M., and Hinton, G. (2020b). A simple framework for contrastive learning of visual representations. In III, H. D. and Singh, A., editors, *Proceedings of the 37th International Conference on Machine Learning (ICML 2020)*, volume 119 of *Proceedings of Machine Learning Research*, pages 1597–1607. PMLR.
[Chen et al., 2022] Chen, Y., Guo, Z., Yin, Q., Chen, H., and Huang, K. (2022). Layer-wisely supervised learning for one-shot neural architecture search. In *2022 International Joint Conference on Neural Networks (IJCNN)*, pages 1–8.

[Chernousko and Lyubushin, 1982] Chernousko, F. L. and Lyubushin, A. A. (1982). Method of successive approximations for solution of optimal control problems. *Optimal Control Applications and Methods*, 3(2):101–114.

[Chizat, 2020] Chizat, L. (2020). Sparse optimization on measures with over-parameterized gradient descent. *arXiv*.

[Chizat and Bach, 2018] Chizat, L. and Bach, F. (2018). On the global convergence of gradient descent for over-parameterized models using optimal transport. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Ciccone et al., 2018] Ciccone, M., Gallieri, M., Masci, J., Osendorfer, C., and Gomez, F. (2018). Nais-net: Stable deep networks from non-autonomous differential equations. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Cisse et al., 2017] Cisse, M., Bojanowski, P., Grave, E., Dauphin, Y., and Usunier, N. (2017). Parseval networks: Improving robustness to adversarial examples. In *Proceedings of the 34th International Conference on Machine Learning (ICML 2017)*, page 854–863. PMLR.

[Coate et al., 2011] Coate, A., Lee, H., and Ng, A. Y. (2011). An analysis of single layer networks in unsupervised feature learning. In *International Conference on Artificial Intelligence and Statistics (AISTATS 2011)*.

[Croce and Hein, 2020a] Croce, F. and Hein, M. (2020a). Minimally distorted adversarial examples with a fast adaptive boundary attack. In *Proceedings of the 37th International Conference on Machine Learning (ICML 2020)*. PMLR.

[Croce and Hein, 2020b] Croce, F. and Hein, M. (2020b). Reliable evaluation of adversarial robustness with an ensemble of diverse parameter-free attacks. In *Proceedings of the 37th International Conference on Machine Learning (ICML 2020)*. PMLR.

[Cybenko, 1989] Cybenko, G. (1989). Approximation by superpositions of a sigmoidal function. *Mathematics of Control, Signals and Systems*, 2:303–314.

[Czarnecki et al., 2017] Czarnecki, W. M., Świrszcz, G., Jaderberg, M., Osindero, S., Vinyals, O., and Kavukcuoglu, K. (2017). Understanding synthetic gradients and decoupled neural interfaces. In *Proceedings of the International Conference on Machine Learning (ICML 2017)*.

[Dahlquist, 1963] Dahlquist, G. G. (1963). A special stability problem for linear multistep methods. *BIT Numerical Mathematics*, 3:27–43.
[De and Smith, 2020] De, S. and Smith, S. (2020). Batch normalization biases residual blocks towards the identity function in deep networks. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 19964–19975. Curran Associates, Inc.

[de Bezenac et al., 2018] de Bezenac, E., Pajot, A., and Gallinari, P. (2018). Deep learning for physical processes: Incorporating prior scientific knowledge. In *International Conference on Learning Representations (ICLR 2018)*.

[de Bézenac et al., 2021] de Bézenac, E., Ayed, I., and Gallinari, P. (2021). Cyclegan through the lens of (dynamical) optimal transport. In *European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD 2021)*.

[De Palma et al., 2019] De Palma, G., Kiani, B., and Lloyd, S. (2019). Random deep neural networks are biased towards simple functions. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Dohmatob, 2019] Dohmatob, E. (2019). Limitations of adversarial robustness: Strong no free lunch theorem. In *Proceedings of the 36th International Conference on Machine Learning (ICML 2019)*, Proceedings of Machine Learning Research. PMLR.

[Dosovitskiy et al., 2021] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weis-senborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N. (2021). An image is worth 16x16 words: Transformers for image recognition at scale. In *International Conference on Learning Representations (ICLR 2021)*.

[Du and Goel, 2018] Du, S. S. and Goel, S. (2018). Improved learning of one-hidden-layer convolutional neural networks with overlaps. *arXiv*.

[Du et al., 2021] Du, X., Farrahi, K., and Niranjan, M. (2021). Information bottleneck theory based exploration of cascade learning. *Entropy*, 23(10).

[Duan and Príncipe, 2022] Duan, S. and Príncipe, J. C. (2022). Training deep architectures without end-to-end backpropagation: A survey on the provably optimal methods. *IEEE Computational Intelligence Magazine*, 17(4):39–51.

[Duan et al., 2022] Duan, S., Yu, S., and Príncipe, J. C. (2022). Modularizing deep learning via pairwise learning with kernels. *IEEE Transactions on Neural Networks and Learning Systems*, 33(4):1441–1451.

[Dun et al., 2021] Dun, C., Wolfe, C. R., Jermaine, C. M., and Kyrillidis, A. (2021). Resist: Layer-wise decomposition of resnets for distributed training. *arXiv*.
[Engstrom et al., 2019] Engstrom, L., Tran, B., Tsipras, D., Schmidt, L., and Madry, A. (2019). Exploring the landscape of spatial robustness. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 1802–1811. PMLR.

[Eniser et al., 2020] Eniser, H. F., Christakis, M., and Wüstholz, V. (2020). *Raid: Randomized adversarial-input detection for neural networks*. arXiv.

[et al., 2016] et al., Y. W. (2016). Google's neural machine translation system: *Bridging the gap between human and machine translation*. ArXiv.

[Eykholt et al., 2018] Eykholt, K., Evtimov, I., Fernandes, E., Li, B., Rahmati, A., Xiao, C., Prakash, A., Kohno, T., and Song, D. (2018). Robust physical-world attacks on deep learning visual classification. In *2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 1625–1634, Los Alamitos, CA, USA. IEEE Computer Society.

[Feinman et al., 2017] Feinman, R., Curtin, R. R., Shintre, S., and Gardner, A. B. (2017). Detecting adversarial samples from artifacts. arXiv.

[Feynman, 1942] Feynman, R. P. (1942). *The principle of least action in quantum mechanics*. PhD thesis, Princeton U.

[Figalli, 2017] Figalli, A. (2017). *The Monge-Ampere Equation and Its Applications*. Zurich lectures in advanced mathematics. European Mathematical Society.

[Figurnov et al., 2017] Figurnov, M., Collins, M. D., Zhu, Y., Zhang, L., Huang, J., Vetrov, D., and Salakhutdinov, R. (2017). Spatially adaptive computation time for residual networks. In *2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 1790–1799.

[Finlay et al., 2020] Finlay, C., Jacobsen, J.-H., Nurbekyan, L., and Oberman, A. M. (2020). How to train your neural ode. In *Proceedings of the International Conference on Machine Learning (ICML 2020)*.

[Fojo et al., 2018] Fojo, D., Campos, V., and i Nieto, X. G. (2018). Comparing fixed and adaptive computation time for recurrent neural networks. In *Workshop at International Conference on Learning Representations (ICLR 2018)*.

[Fournier et al., 2023] Fournier, L., Rivaud, S., Belilovsky, E., Eickenberg, M., and Oyallon, E. (2023). Can forward gradient match backpropagation? In *Proceedings of the International Conference on Machine Learning (ICML 2023)*.

[Freund and Schapire, 1997] Freund, Y. and Schapire, R. E. (1997). A decision-theoretic generalization of on-line learning and an application to boosting. *Journal of Computer and System Sciences*, 55(1):119–139.
[Fukushima, 1975] Fukushima, K. (1975). Cognitron: A self-organizing multi-layered neural network. *Biological Cybernetics*, 20.

[Gai and Zhang, 2021] Gai, K. and Zhang, S. (2021). A mathematical principle of deep learning: Learn the geodesic curve in the wasserstein space. *arXiv*.

[Gao et al., 2019] Gao, Y., Jiao, Y., Wang, Y., Wang, Y., Yang, C., and Zhang, S. (2019). Deep generative learning via variational gradient flow. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[García-Morales et al., 2008] García-Morales, V., Pellicer, J., and Manzanares, J. A. (2008). Thermodynamics based on the principle of least abbreviated action: Entropy production in a network of coupled oscillators. *Annals of Physics*, 323(8):1844–1858.

[Ge et al., 2018] Ge, R., Lee, J. D., and Ma, T. (2018). Learning one-hidden-layer neural networks with landscape design. In *International Conference on Learning Representations (ICLR 2018)*.

[Gianazza et al., 1994] Gianazza, U., Gobbino, M., and Savarè, G. (1994). Evolution problems and minimizing movements. *Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche, Matematiche e Naturali. Rendiconti Lincei. Matematica e Applicazioni*, 5:289-296.

[Gilmer et al., 2018] Gilmer, J., Metz, L., Faghri, F., Schoenholz, S. S., Raghu, M., Wattenberg, M., and Goodfellow, I. (2018). Adversarial spheres: The relationship between high-dimensional geometry and adversarial examples. *arXiv*.

[Giorgi et al., 1980] Giorgi, E. D., Marino, A., and Tosques, M. (1980). Problemi di evoluzione in spazi metrici e curve di massima pendenza. *Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche, Matematiche e Naturali. Rendiconti*, 68(3):180-187.

[Gomez et al., 2020] Gomez, A. N., Key, O., Gou, S., Frosst, N., Dean, J., and Gal, Y. (2020). Interlocking backpropagation: Improving depthwise model-parallelism. *Journal of Machine Learning Research 23* (JMLR 23), 23:171:1–171:28.

[Gomez et al., 2017] Gomez, A. N., Ren, M., Urtasun, R., and Grosse, R. B. (2017). The reversible residual network: Backpropagation without storing activations. In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2017)*, volume 30. Curran Associates, Inc.

[Goodfellow et al., 2016] Goodfellow, I., Bengio, Y., and Courville, A. (2016). *Deep Learning*. MIT Press.

[Goodfellow et al., 2015a] Goodfellow, I. J., Shlens, J., and Szegedy, C. (2015a). Explaining and harnessing adversarial examples. In *International Conference on Learning Representations (ICLR 2015)*.
[Goodfellow et al., 2015b] Goodfellow, I. J., Shlens, J., and Szegedy, C. (2015b). Explaining and harnessing adversarial examples. In *International Conference on Learning Representations (ICLR 2015)*.

[Graves, 2017] Graves, A. (2017). Adaptive computation time for recurrent neural networks. *arXiv*.

[Gray, 2009] Gray, C. G. (2009). Principle of least action. *Scholarpedia*, 4(12):8291.

[Greff et al., 2016] Greff, K., Srivastava, R. K., and Schmidhuber, J. (2016). Highway and residual networks learn unrolled iterative estimation. *arXiv*.

[Gretton et al., 2012] Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. (2012). A kernel two-sample test. *Journal of Machine Learning Research* 13 (JMLR 13), 13.

[Greydanus et al., 2019] Greydanus, S., Dzamba, M., and Yosinski, J. (2019). Hamiltonian neural networks. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Grosse et al., 2017] Grosse, K., Manoharan, P., Papernot, N., Backes, M., and McDaniel, P. (2017). On the (statistical) detection of adversarial examples. *arXiv*.

[Gu and Rigazio, 2015] Gu, S. and Rigazio, L. (2015). Towards deep neural network architectures robust to adversarial examples. In *Workshop at International Conference on Learning Representations (ICLR 2015)*.

[Guo et al., 2022a] Guo, W., Fouda, M. E., Eltawil, A. M., and Salama, K. N. (2022a). Backlink: Supervised local training with backward links. *arXiv*.

[Guo et al., 2022b] Guo, W., Fouda, M. E., Eltawil, A. M., and Salama, K. N. (2022b). Efficient training of spiking neural networks with temporally-truncated local backpropagation through time. *arXiv*.

[Gupta, 2020] Gupta, S. K. (2020). A more biologically plausible local learning rule for anns. *arXiv*.

[Haber et al., 2019] Haber, E., Lensink, K., Treister, E., and Ruthotto, L. (2019). Imexnet - a forward stable deep neural network. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[Haber and Ruthotto, 2017] Haber, E. and Ruthotto, L. (2017). Stable architectures for deep neural networks. *CoRR*, abs/1705.03341.

[Harder et al., 2021] Harder, P., Pfreundt, F.-J., Keuper, M., and Keuper, J. (2021). Spectraldefense: Detecting adversarial attacks on cnns in the fourier domain. *arXiv*.
[Hastie et al., 2001] Hastie, T., Tibshirani, R., and Friedman, J. (2001). *The Elements of Statistical Learning*. Springer-Verlag New York.

[Hauser, 2019] Hauser, M. (2019). On residual networks learning a perturbation from identity. *arXiv*.

[Hayou and Ayed, 2021] Hayou, S. and Ayed, F. (2021). Regularization in resnet with stochastic depth. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, *Advances in Neural Information Processing Systems (NeurIPS 2021)*, volume 34, pages 15464–15474. Curran Associates, Inc.

[He et al., 2020] He, K., Fan, H., Wu, Y., Xie, S., and Girshick, R. (2020). Momentum contrast for unsupervised visual representation learning. In *2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 9726–9735, Los Alamitos, CA, USA. IEEE Computer Society.

[He et al., 2016a] He, K., Zhan, X., Ren, S., and Sun, J. (2016a). Identity mappings in deep residual networks. In *European Conference on Computer Vision (ECCV 2016)*.

[He et al., 2015] He, K., Zhang, X., Ren, S., and Sun, J. (2015). Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In *2015 IEEE International Conference on Computer Vision (ICCV)*, pages 1026–1034.

[He et al., 2016b] He, K., Zhang, X., Ren, S., and Sun, J. (2016b). Deep residual learning for image recognition. In *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2016)*.

[Hein and Andriushchenko, 2017] Hein, M. and Andriushchenko, M. (2017). Formal guarantees on the robustness of a classifier against adversarial manipulation. In *Advances in Neural Information Processing Systems (NeurIPS 2017)*, page 2263–2273, Red Hook, NY, USA. Curran Associates Inc.

[Hinton et al., 1995] Hinton, G. E., Dayan, P., Frey, B. J., and Neal, R. M. (1995). The "wake-sleep" algorithm for unsupervised neural networks. *Science*, 268(5214):1158–1161.

[Hochreiter and Schmidhuber, 1997] Hochreiter, S. and Schmidhuber, J. (1997). Long short-term memory. *Neural computation*, 9:1735–80.

[Hornik et al., 1989] Hornik, K., Stinchcombe, M., and White, H. (1989). Multilayer feedforward networks are universal approximators. *Neural Networks*, 2(5):359–366.

[Houle, 2013] Houle, M. E. (2013). Dimensionality, discriminability, density and distance distributions. In *Workshop at 2013 IEEE 13th International Conference on Data Mining*, pages 468–473.
[Hu and Tan, 2022] Hu, W. and Tan, Y. (2022). Generating adversarial malware examples for black-box attacks based on gan. In Tan, Y. and Shi, Y., editors, *Data Mining and Big Data*, pages 409–423, Singapore. Springer Nature Singapore.

[Huang et al., 2018] Huang, F., Ash, J., Langford, J., and Schapire, R. (2018). Learning deep resnet blocks sequentially using boosting theory. In *Proceedings of the International Conference on Machine Learning (ICML 2018)*.

[Huang et al., 2016] Huang, G., Sun, Y., Liu, Z., Sedra, D., and Weinberger, K. (2016). Deep networks with stochastic depth. *arXiv*.

[Huo et al., 2018a] Huo, Z., Gu, B., and Huang, H. (2018a). Training neural networks using features replay. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Huo et al., 2018b] Huo, Z., Gu, B., Yang, Q., and Huang, H. (2018b). Decoupled parallel backpropagation with convergence guarantee. In *Proceedings of the International Conference on Machine Learning (ICML 2018)*.

[Huster et al., 2019] Huster, T., Chiang, C.-Y. J., and Chadha, R. (2019). Limitations of the lipschitz constant as a defense against adversarial examples. In *Workshop at European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD 2019)*, pages 16–29. Springer International Publishing.

[Ioffe and Szegedy, 2015] Ioffe, S. and Szegedy, C. (2015). Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Bach, F. and Blei, D., editors, *Proceedings of the 32nd International Conference on Machine Learning (ICML 2015)*, volume 37 of *Proceedings of Machine Learning Research*, pages 448–456, Lille, France. PMLR.

[Jacobsen et al., 2018] Jacobsen, J.-H., Smeulders, A., and Oyallon, E. (2018). i-revnet: Deep invertible networks. In *International Conference on Learning Representations (ICLR 2018)*.

[Jacot et al., 2018] Jacot, A., Gabriel, F., and Hongler, C. (2018). Neural tangent kernel: Convergence and generalization in neural networks. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Jaderberg et al., 2017] Jaderberg, M., Czarnecki, W. M., Osindero, S., Vinyals, O., Graves, A., Silver, D., and Kavukcuoglu, K. (2017). Decoupled neural interfaces using synthetic gradients. In *Proceedings of the International Conference on Machine Learning (ICML 2017)*.
[Janner et al., 2022] Janner, M., Du, Y., Tenenbaum, J., and Levine, S. (2022). Planning with diffusion for flexible behavior synthesis. In *Proceedings of the 39th International Conference on Machine Learning (ICML)*.

[Janzamin et al., 2016] Janzamin, M., Sedghi, H., and Anandkumar, A. (2016). Beating the perils of non-convexity: Guaranteed training of neural networks using tensor methods. *arXiv*.

[Jastrzebski et al., 2018] Jastrzebski, S., Arpit, D., Ballas, N., Verma, V., Che, T., and Bengio, Y. (2018). Residual connections encourage iterative inference. In *International Conference on Learning Representations (ICLR 2018)*.

[Jernite et al., 2017] Jernite, Y., Grave, E., Joulin, A., and Mikolov, T. (2017). Variable computation in recurrent neural networks. In *International Conference on Learning Representations (ICLR 2017)*.

[Jordan et al., 1998] Jordan, R., Kinderlehrer, D., and Otto, F. (1998). The variational formulation of the fokker-planck equation. *SIAM Journal on Mathematical Analysis*, 29(1):1–17.

[Jumper et al., 2021] Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Žídek, A., Potapenko, A., Bridgland, A., Meyer, C., Kohl, S. A. A., Ballard, A. J., Cowie, A., Romera-Paredes, B., Nikolov, S., Jain, R., Adler, J., Back, T., Petersen, S., Reiman, D., Clancy, E., Zielinski, M., and Hassabis, D. (2021). Highly accurate protein structure prediction with alphafold. *Nature*.

[Kaila and Annila, 2008] Kaila, V. R. and Annila, A. (2008). Natural selection for least action. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 464(2099):3055–3070.

[Kantchelian et al., 2016] Kantchelian, A., Tygar, J. D., and Joseph, A. D. (2016). Evasion and hardening of tree ensemble classifiers. In *Proceedings of the 33rd International Conference on Machine Learning (ICML 2016)*. PMLR.

[Karkar et al., 2020] Karkar, S., Ayed, I., de Bézenac, E., and Gallinari, P. (2020). A principle of least action for the training of neural networks. In *European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD 2020)*.

[Karkar et al., 2023a] Karkar, S., Ayed, I., de Bézenac, E., and Gallinari, P. (2023a). Module-wise training of neural networks via the minimizing movement scheme. In *Advances in Neural Information Processing Systems (NeurIPS 2023)*.

[Karkar et al., 2023b] Karkar, S., Gallinari, P., and Rakotomamonjy, A. (2023b). Adversarial sample detection through neural network transport dynamics. In *European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD 2023)*.
[Kherchouche et al., 2022] Kherchouche, A., Fezza, S. A., and Hamidouche, W. (2022). Detect and defense against adversarial examples in deep learning using natural scene statistics and adaptive denoising. *Neural Computing and Applications*, 34(24):21567–21582.

[Kherchouche et al., 2020] Kherchouche, A., Fezza, S. A., Hamidouche, W., and Déforges, O. (2020). Detection of adversarial examples in deep neural networks with natural scene statistics. In *2020 International Joint Conference on Neural Networks (IJCNN)*, pages 1–7.

[Khosla et al., 2020] Khosla, P., Teterwak, P., Wang, C., Sarna, A., Tian, Y., Isola, P., Maschinot, A., Liu, C., and Krishnan, D. (2020). Supervised contrastive learning. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 18661–18673. Curran Associates, Inc.

[Khoury and Hadfield-Menell, 2018] Khoury, M. and Hadfield-Menell, D. (2018). On the geometry of adversarial examples. *arXiv*.

[Kolouri et al., 2017] Kolouri, S., Park, S., Thorpe, M., Slepčev, D., and Rohde, G. K. (2017). Optimal mass transport: Signal processing and machine-learning applications. *IEEE signal processing magazine*.

[Kombarov, 2021] Kombarov, S. (2021). Action in economics: Mathematical derivation of laws of economics from the principle of least action in physics. In *Proceedings of International Conference of Eurasian Economies*.

[Kornblith et al., 2019] Kornblith, S., Norouzi, M., Lee, H., and Hinton, G. (2019). Similarity of neural network representations revisited. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning (ICML 2019)*, volume 97 of *Proceedings of Machine Learning Research*, pages 3519–3529. PMLR.

[Kotyan and Vargas, 2020] Kotyan, S. and Vargas, D. V. (2020). Adversarial robustness assessment: Why both $l_0$ and $l_\infty$ attacks are necessary. *arXiv*.

[Krizhevsky, 2009] Krizhevsky, A. (2009). Learning multiple layers of features from tiny images. *University of Toronto Technical Report*.

[Krizhevsky et al., 2012] Krizhevsky, A., Sutskever, I., and Hinton, G. E. (2012). Imagenet classification with deep convolutional neural networks. In Pereira, F., Burges, C., Bottou, L., and Weinberger, K., editors, *Advances in Neural Information Processing Systems (NeurIPS 2012)*, volume 25. Curran Associates, Inc.

[Kurakin et al., 2017] Kurakin, A., Goodfellow, I., and Bengio, S. (2017). Adversarial examples in the physical world. In *Workshop at International Conference on Learning Representations (ICLR 2017)*.
[Laakso and Cottrell, 2000] Laakso, A. and Cottrell, G. (2000). Content and cluster analysis: Assessing representational similarity in neural systems. *Philosophical Psychology*, 13(1):47–76.

[Larsson et al., 2017] Larsson, G., Maire, M., and Shakhnarovich, G. (2017). Fractalnet: Ultra-deep neural networks without residuals. In *International Conference on Learning Representations (ICLR 2017)*.

[Laskin et al., 2021] Laskin, M., Metz, L., Nabarro, S., Saroufim, M., Nouné, B., Luschi, C., Sohl-Dickstein, J., and Abbeel, P. (2021). Parallel training of deep networks with local updates. *arXiv*.

[Latorre et al., 2020] Latorre, F., Rolland, P., and Cevher, V. (2020). Lipschitz constant estimation of neural networks via sparse polynomial optimization. In *International Conference on Learning Representations (ICLR 2020)*.

[LeCun et al., 2010] LeCun, Y., Cortes, C., and Burges, C. J. (2010). MNIST handwritten digit database. *yann.lecun.com/exdb/mnist*.

[Lee et al., 2018] Lee, K., Lee, K., Lee, H., and Shin, J. (2018). A simple unified framework for detecting out-of-distribution samples and adversarial attacks. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Lengellé and Denœux, 1996] Lengellé, R. and Denœux, T. (1996). Training mlps layer by layer using an objective function for internal representations. *Neural Networks*, 9(1):83–97.

[Li et al., 2018] Li, Q., Chen, L., Tai, C., and E, W. (2018). Maximum principle based algorithms for deep learning. *Journal of Machine Learning Research (JMLR 2018)*.

[Li et al., 2022] Li, Y., Cheng, M., Hsieh, C.-J., and Lee, T. C. M. (2022). A review of adversarial attack and defense for classification methods. *The American Statistician*, 76(4):329–345.

[Li et al., 2015] Li, Y., Yosinski, J., Clune, J., Lipson, H., and Hopcroft, J. (2015). Convergent learning: Do different neural networks learn the same representations? In Storcheus, D., Rostamizadeh, A., and Kumar, S., editors, *Proceedings of the 1st International Workshop on Feature Extraction: Modern Questions and Challenges at Advances in Neural Information Processing Systems (NeurIPS 2015)*, volume 44 of *Proceedings of Machine Learning Research*, pages 196–212, Montreal, Canada. PMLR.

[Liao et al., 2016] Liao, Q., Leibo, J., and Poggio, T. (2016). How important is weight symmetry in backpropagation? *Proceedings of the AAAI Conference on Artificial Intelligence*, 30(1).

[Liao and Poggio, 2016] Liao, Q. and Poggio, T. (2016). Bridging the gaps between residual learning, recurrent neural networks and visual cortex. *arXiv*.
[Liu et al., 2017] Liu, Y., Chen, X., Liu, C., and Song, D. (2017). Delving into transferable adversarial examples and black-box attacks. In *International Conference on Learning Representations (ICLR 2017)*. OpenReview.net.

[Liu et al., 2021] Liu, Z., Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang, Z., Lin, S., and Guo, B. (2021). Swin transformer: Hierarchical vision transformer using shifted windows. In *Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)*.

[Liutkus et al., 2019] Liutkus, A., Imşekli, U., Majewski, S., Durmus, A., and Stoter, F.-R. (2019). Sliced-wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[Löwe et al., 2019] Löwe, S., O' Connor, P., and Veeling, B. (2019). Putting an end to end-to-end: Gradient-isolated learning of representations. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Lu et al., 2018] Lu, Y., Zhong, A., Li, Q., and Dong, B. (2018). Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. In *Proceedings of the International Conference on Machine Learning (ICML 2018)*.

[Ma et al., 2020] Ma, W.-D. K., Lewis, J., and Kleijn, W. B. (2020). The hsic bottleneck: Deep learning without back-propagation. In *The AAAI Conference on Artificial Intelligence (AAAI-20)*.

[Ma et al., 2018] Ma, X., Li, B., Wang, Y., Erfani, S. M., Wijewickrema, S., Schoenebeck, G., Song, D., Houle, M. E., and Bailey, J. (2018). Characterizing adversarial subspaces using local intrinsic dimensionality. In *International Conference on Learning Representations (ICLR 2018)*.

[Madry et al., 2018] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A. (2018). Towards deep learning models resistant to adversarial attacks. In *International Conference on Learning Representations (ICLR 2018)*.

[Mahajan et al., 2018] Mahajan, D., Girshick, R., Ramanathan, V., He, K., Paluri, M., Li, Y., Bharambe, A., and van der Maaten, L. (2018). Exploring the limits of weakly supervised pretraining. In *European Conference on Computer Vision (ECCV 2018)*.

[Malach and Shalev-Shwartz, 2018] Malach, E. and Shalev-Shwartz, S. (2018). A provably correct algorithm for deep learning that actually works. *arXiv*.

[Mandar Kulkarni, 2016] Mandar Kulkarni, S. K. (2016). Layer-wise training of deep networks using kernel similarity. In *DLPR workshop at the 25th International Conference on Pattern Recognition (ICPR 2016)*.
[Marquez et al., 2018] Marquez, E. S., Hare, J. S., and Niranjan, M. (2018). Deep cascade learning. *IEEE Transactions on Neural Networks and Learning Systems*.

[Masoumian et al., 2022] Masoumian, A., Rashwan, H. A., Cristiano, J., Asif, M. S., and Puig, D. (2022). Monocular depth estimation using deep learning: A review. *Sensors*, 22(14).

[McCann and Pass, 2019] McCann, R. J. and Pass, B. (2019). Optimal transportation between unequal dimensions. *arXiv*.

[Mehrish et al., 2023] Mehrish, A., Majumder, N., Bhardwaj, R., Mihalcea, R., and Poria, S. (2023). A review of deep learning techniques for speech processing. *arXiv*.

[Meng and Chen, 2017] Meng, D. and Chen, H. (2017). Magnet: A two-pronged defense against adversarial examples. In *Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security*, CCS '17, page 135-147, New York, NY, USA. Association for Computing Machinery.

[Metzen et al., 2017] Metzen, J. H., Genewein, T., Fischer, V., and Bischoff, B. (2017). On detecting adversarial perturbations. In *ICLR*.

[Minaee et al., 2022] Minaee, S., Boykov, Y., Porikli, F., Plaza, A., Kehtarnavaz, N., and Terzopoulos, D. (2022). Image segmentation using deep learning: A survey. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 44(7):3523-3542.

[Mittal et al., 2012] Mittal, A., Moorthy, A. K., and Bovik, A. C. (2012). No-reference image quality assessment in the spatial domain. *IEEE Transactions on Image Processing*, 21(12):4695-4708.

[Moosavi-Dezfooli et al., 2017a] Moosavi-Dezfooli, S., Fawzi, A., Fawzi, O., and Frossard, P. (2017a). Universal adversarial perturbations. In *2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 86-94, Los Alamitos, CA, USA. IEEE Computer Society.

[Moosavi-Dezfooli et al., 2017b] Moosavi-Dezfooli, S.-M., Fawzi, A., Fawzi, O., and Frossard, P. (2017b). Universal adversarial perturbations. In *2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 86-94.

[Moosavi-Dezfooli et al., 2016] Moosavi-Dezfooli, S.-M., Fawzi, A., and Frossard, P. (2016). Deepfool: A simple and accurate method to fool deep neural networks. *2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 2574-2582.

[Morcos et al., 2018] Morcos, A., Raghu, M., and Bengio, S. (2018). Insights on representational similarity in neural networks with canonical correlation. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N.,
and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Mostafa et al., 2018] Mostafa, H., Ramesh, V., and Cauwenberghs, G. (2018). Deep supervised learning using local errors. *Frontiers in Neuroscience*.

[Nakkiran et al., 2020] Nakkiran, P., Kaplun, G., Bansal, Y., Yang, T., Barak, B., and Sutskever, I. (2020). Deep double descent: Where bigger models and more data hurt. In *International Conference on Learning Representations (ICLR 2020)*.

[Netzer et al., 2011] Netzer, Y., Wang, T., Coates, A., Bissacco, A., Wu, B., and Ng, A. Y. (2011). Reading digits in natural images with unsupervised feature learning. In *Workshop on Deep Learning and Unsupervised Feature Learning at Advances in Neural Information Processing Systems (NeurIPS 2011)*.

[Nguyen et al., 2015] Nguyen, A., Yosinski, J., and Clune, J. (2015). Deep neural networks are easily fooled: High confidence predictions for unrecognizable images. In *2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 427–436.

[Nguyen and Choi, 2019] Nguyen, T. T. and Choi, J. (2019). Layer-wise learning of stochastic neural networks with information bottleneck. *Entropy*, 21.

[Ni et al., 2022] Ni, Z., Wang, Y., Yu, J., Jiang, H., Cao, Y., and Huang, G. (2022). Deep model assembling. *arXiv*.

[Nicolae et al., 2018] Nicolae, M.-I., Sinn, M., Tran, M. N., Buesser, B., Rawat, A., Wistuba, M., Zantedeschi, V., Baracaldo, N., Chen, B., Ludwig, H., Molloy, I. M., and Edwards, B. (2018). Adversarial robustness toolbox v1.0.0. *arXiv*.

[Novak et al., 2018] Novak, R., Bahri, Y., Abolafia, D. A., Pennington, J., and Sohl-Dickstein, J. (2018). Sensitivity and generalization in neural networks: an empirical study. In *International Conference on Learning Representations (ICLR 2018)*.

[Nøkland and Eidnes, 2019] Nøkland, A. and Eidnes, L. H. (2019). Training neural networks with local error signals. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[Pang et al., 2018] Pang, T., Du, C., Dong, Y., and Zhu, J. (2018). Towards robust detection of adversarial examples. In *Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Papernot et al., 2016a] Papernot, N., McDaniel, P., and Goodfellow, I. (2016a). Transferability in machine learning: from phenomena to black-box attacks using adversarial samples. *arXiv*.
[Papernot et al., 2017] Papernot, N., McDaniel, P., Goodfellow, I., Jha, S., Celik, Z. B., and Swami, A. (2017). Practical black-box attacks against machine learning. In *Proceedings of the 2017 ACM on Asia Conference on Computer and Communications Security*, ASIA CCS '17, page 506–519, New York, NY, USA. Association for Computing Machinery.

[Papernot et al., 2016b] Papernot, N., McDaniel, P., Jha, S., Fredrikson, M., Celik, Z. B., and Swami, A. (2016b). The limitations of deep learning in adversarial settings. In *2016 IEEE European Symposium on Security and Privacy (EuroSP)*, pages 372–387.

[Papernot et al., 2016c] Papernot, N., McDaniel, P., Wu, X., Jha, S., and Swami, A. (2016c). Distillation as a defense to adversarial perturbations against deep neural networks. In *2016 IEEE Symposium on Security and Privacy (SP)*, pages 582–597.

[Parzen, 1962] Parzen, E. (1962). On Estimation of a Probability Density Function and Mode. *The Annals of Mathematical Statistics*, 33(3):1065 – 1076.

[Pass, 2010] Pass, B. (2010). Regularity of optimal transportation between spaces with different dimensions. *arXiv*.

[Patel et al., 2022] Patel, A., Eickenberg, M., and Belilovsky, E. (2022). Local learning with neuron groups. *From Cells to Societies: Collective Learning Across Scales International Conference on Learning Representations (ICLR 2022)*.

[Pathak et al., 2022] Pathak, P., Zhang, J., and Samaras, D. (2022). Local learning on transformers via feature reconstruction. *arXiv*.

[Pedregosa et al., 2011] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011). Scikit-learn: Machine learning in Python. *Journal of Machine Learning Research (JMLR 12)*, 12:2825–2830.

[Peyre and Cuturi, 2019] Peyre, G. and Cuturi, M. (2019). *Computational Optimal Transport*. Now Publishers.

[Pohlen et al., 2017] Pohlen, T., Hermans, A., Mathias, M., and Leibe, B. (2017). Full-resolution residual networks for semantic segmentation in street scenes. In *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2017)*.

[Pontryagin, 1987] Pontryagin, L. (1987). *Mathematical Theory of Optimal Processes*. Classics of Soviet Mathematics. Taylor & Francis.

[Purwins et al., 2019] Purwins, H., Li, B., Virtanen, T., Schlüter, J., Chang, S.-Y., and Sainath, T. (2019). Deep learning for audio signal processing. *IEEE Journal of Selected Topics in Signal Processing*, 13(2):206–219.
[Pyeon et al., 2021] Pyeon, M., Moon, J., Hahn, T., and Kim, G. (2021). {SE-DONA}: Search for decoupled neural networks toward greedy block-wise learning. In *International Conference on Learning Representations (ICLR 2021)*.

[Quarteroni et al., 2007] Quarteroni, A., Sacco, R., and Saleri, F. (2007). *Numerical Mathematics*. Springer Berlin, Heidelberg.

[Radford et al., 2018] Radford, A., Narasimhan, K., Salimans, T., and Sutskever, I. (2018). Improving language understanding with unsupervised learning. *openai.com/research/language-unsupervised*.

[Raghu et al., 2017] Raghu, M., Gilmer, J., Yosinski, J., and Sohl-Dickstein, J. (2017). Svcca: Singular vector canonical correlation analysis for deep learning dynamics and interpretability. In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2017)*, volume 30. Curran Associates, Inc.

[Raghunathan et al., 2020] Raghunathan, A., Xie, S. M., Yang, F., Duchi, J., and Liang, P. (2020). Understanding and mitigating the tradeoff between robustness and accuracy. In III, H. D. and Singh, A., editors, *Proceedings of the 37th International Conference on Machine Learning (ICML 2020)*, volume 119 of *Proceedings of Machine Learning Research*, pages 7909–7919. PMLR.

[Rahaman et al., 2019] Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., and Courville, A. (2019). On the spectral bias of neural networks. In *Proceedings of the International Conference on Machine Learning (ICML 2019)*.

[Ren et al., 2021] Ren, J., Fort, S., Liu, J., Roy, A. G., Padhy, S., and Lakshmi-narayanan, B. (2021). A simple fix to mahalanobis distance for improving near-ood detection. *arXiv*.

[Ren et al., 2023] Ren, M., Kornblith, S., Liao, R., and Hinton, G. (2023). Scaling forward gradient with local losses. In *The Eleventh International Conference on Learning Representations (ICLR 2023)*.

[Rosenblatt, 1956] Rosenblatt, M. (1956). Remarks on Some Nonparametric Estimates of a Density Function. *The Annals of Mathematical Statistics*, 27(3):832 – 837.

[Ruiz-Balet and Zuazua, 2021] Ruiz-Balet, D. and Zuazua, E. (2021). Neural ode control for classification, approximation and transport. *arXiv*.

[Rumelhart et al., 1986] Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986). Learning representations by back-propagating errors. *Nature*, 323.

[Russell and Norvig, 2020] Russell, S. and Norvig, P. (2020). *Artificial Intelligence: A Modern Approach, 4th edition*. Pearson Series in Artificial Intelligence., Pearson.
[Ruthotto and Haber, 2018] Ruthotto, L. and Haber, E. (2018). Deep neural networks motivated by partial differential equations. *arXiv*.

[Sabour et al., 2016] Sabour, S., Cao, Y., Faghri, F., and Fleet, D. J. (2016). Adversarial manipulation of deep representations. In *International Conference on Learning Representations (ICLR 2016)*.

[Samangouei et al., 2018] Samangouei, P., Kabkab, M., and Chellappa, R. (2018). Defense-GAN: Protecting classifiers against adversarial attacks using generative models. In *International Conference on Learning Representations (ICLR 2018)*.

[Sandler et al., 2019] Sandler, M., Baccash, J., Zhmoginov, A., and Howard, A. (2019). Non-discriminative data or weak model? on the relative importance of data and model resolution. In *Workshop on Real-World Recognition from Low-Quality Images and Videos at International Conference on Computer Vision (ICCV 2019)*.

[Sandler et al., 2018] Sandler, M., Howard, A. G., Zhu, M., Zhmoginov, A., and Chen, L. (2018). Mobilenetv2: Inverted residuals and linear bottlenecks. In *2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018*, pages 4510–4520. Computer Vision Foundation / IEEE Computer Society.

[Santambrogio, 2015] Santambrogio, F. (2015). *Optimal Transport for Applied Mathematicians*. Birkhäuser.

[Santambrogio, 2016] Santambrogio, F. (2016). Euclidean, metric, and wasserstein gradient flows: an overview. *arXiv*.

[Saxe et al., 2014] Saxe, A. M., Mcclelland, J. L., and Ganguli, S. (2014). Exact solutions to the nonlinear dynamics of learning in deep linear neural network. In *International Conference on Learning Representations (ICLR 2014)*.

[Schmidt et al., 2018] Schmidt, L., Santurkar, S., Tsipras, D., Talwar, K., and Madry, A. (2018). Adversarially robust generalization requires more data. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Shannon, 1948] Shannon, C. E. (1948). A mathematical theory of communication. *The Bell System Technical Journal*, 27(4):623–656.

[Sharma et al., 2021] Sharma, V., Gupta, M., Kumar, A., and Mishra, D. (2021). Video processing using deep learning techniques: A systematic literature review. *IEEE Access*, 9:139489–139507.

[Shashkin, 1991] Shashkin, I. (1991). *Fixed Points*. Mathematical world. American Mathematical Soc.
[Silver et al., 2016] Silver, D., Huang, A., Maddison, C. J., Guez, A., Sifre, L., van den Driessche, G., Schrittwieser, J., Antonoglou, I., Panneershelvam, V., Lanctot, M., Dieleman, S., Grewe, D., Nham, J., Kalchbrenner, N., Sutskever, I., Lillicrap, T., Leach, M., Kavukcuoglu, K., Graepel, T., and Hassabis, D. (2016). Mastering the game of go with deep neural networks and tree search. *nature. Nature*.

[Sim et al., 2020] Sim, B., Oh, G., Kim, J., Jung, C., and Ye, J. C. (2020). Optimal transport driven cyclegan for unsupervised learning in inverse problems. *SIAM Journal on Imaging Sciences*, 13(4):2281-2306.

[Simonyan and Zisserman, 2014] Simonyan, K. and Zisserman, A. (2014). Very deep convolutional networks for large-scale image recognition. In *International Conference on Learning Representations (ICLR 2014)*.

[Singh and Jaggi, 2020] Singh, S. P. and Jaggi, M. (2020). Model fusion via optimal transport. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 22045-22055. Curran Associates, Inc.

[Song et al., 2018] Song, Y., Kim, T., Nowozin, S., Ermon, S., and Kushman, N. (2018). Pixeldefend: Leveraging generative models to understand and defend against adversarial examples. In *International Conference on Learning Representations (ICLR 2018)*.

[Sonoda and Murata, 2019] Sonoda, S. and Murata, N. (2019). Transport analysis of infinitely deep neural network. *Journal of Machine Learning Research (JMLR 2019)*.

[Srivastava et al., 2014] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014). Dropout: A simple way to prevent neural networks from overfitting. *Journal of Machine Learning Research (JMLR 15)*, 15(56):1929-1958.

[Stich, 2019] Stich, S. U. (2019). Local sgd converges fast and communicates little. In *International Conference on Learning Representations (ICLR 2019)*.

[Su et al., 2019] Su, J., Vargas, D. V., and Sakurai, K. (2019). One pixel attack for fooling deep neural networks. *IEEE Transactions on Evolutionary Computation*, 23(5):828-841.

[Sun et al., 2021] Sun, Q., Dong, H., Chen, Z., Sun, J., Li, Z., and Dong, B. (2021). Layer-parallel training of residual networks with auxiliary-variable networks. *arXiv*.

[Szegedy et al., 2017] Szegedy, C., Ioffe, S., Vanhoucke, V., and Alemi, A. A. (2017). Inception-v4, inception-resnet and the impact of residual connections on learning. In *Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence*, AAAI’17, page 4278-4284. AAAI Press.
[Szegedy et al., 2013] Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I., and Fergus, R. (2013). Intriguing properties of neural networks. *arXiv*.

[Söderlind, 2006] Söderlind, G. (2006). Time-step selection algorithms: Adaptivity, control, and signal processing. *Applied Numerical Mathematics*, 56(3):488–502. Selected Papers, The Third International Conference on the Numerical Solutions of Volterra and Delay Equations.

[Tan and Le, 2019] Tan, M. and Le, Q. (2019). EfficientNet: Rethinking model scaling for convolutional neural networks. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning (ICML 2019)*, volume 97 of *Proceedings of Machine Learning Research*, pages 6105–6114. PMLR.

[Tanay and Griffin, 2016] Tanay, T. and Griffin, L. (2016). A boundary tilting perspective on the phenomenon of adversarial examples. *arXiv*.

[Tang et al., 2021] Tang, Y., Teng, Q., Zhang, L., Min, F., and He, J. (2021). Layer-wise training convolutional neural networks with smaller filters for human activity recognition using wearable sensors. *IEEE Sensors Journal*.

[Tempczyk et al., 2022] Tempczyk, P., Michaluk, R., Garncarek, L., Spurek, P., Tabor, J., and Golinski, A. (2022). LIDL: Local intrinsic dimension estimation using approximate likelihood. In Chaudhuri, K., Jegelka, S., Song, L., Szepesvari, C., Niu, G., and Sabato, S., editors, *Proceedings of the 39th International Conference on Machine Learning (ICML 2022)*, volume 162 of *Proceedings of Machine Learning Research*, pages 21205–21231. PMLR.

[Teng et al., 2020] Teng, Q., Wang, K., Zhang, L., and He, J. (2020). The layer-wise training convolutional neural networks using local loss for sensor-based human activity recognition. *IEEE Sensors Journal*.

[Tran et al., 2018] Tran, D., Wang, H., Torresani, L., Ray, J., LeCun, Y., and Paluri, M. (2018). A closer look at spatiotemporal convolutions for action recognition. In *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2018)*.

[van den Oord et al., 2016] van den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A., and Kavukcuoglu, K. (2016). Wavenet: A generative model for raw audio. *arXiv*.

[Veit et al., 2016] Veit, A., Wilber, M. J., and Belongie, S. (2016). Residual networks behave like ensembles of relatively shallow networks. In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2016)*, volume 29. Curran Associates, Inc.

[Villani, 2008] Villani, C. (2008). *Optimal Transport: Old and New*. Springer-Verlag.
[Villani, 2021] Villani, C. (2021). *Topics in Optimal Transportation*. Graduate Studies in Mathematics. American Mathematical Society.

[Virmaux and Scaman, 2018] Virmaux, A. and Scaman, K. (2018). Lipschitz regularity of deep neural networks: analysis and efficient estimation. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Wang et al., 2019a] Wang, B., Shi, Z., and Osher, S. (2019a). Resnets ensemble via the feynman-kac formalism to improve natural and robust accuracies. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Wang et al., 2019b] Wang, B., Shi, Z., and Osher, S. (2019b). Resnets ensemble via the feynman-kac formalism to improve natural and robust accuracies. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2019)*, volume 32. Curran Associates, Inc.

[Wang et al., 2018] Wang, L., Hu, L., Gu, J., Hu, Z., Wu, Y., He, K., and Hopcroft, J. (2018). Towards understanding learning representations: To what extent do different neural networks learn the same representation. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, *Advances in Neural Information Processing Systems (NeurIPS 2018)*, volume 31. Curran Associates, Inc.

[Wang et al., 2021] Wang, Y., Ni, Z., Song, S., and Le Yang, G. H. (2021). Revisiting locally supervised learning: an alternative to end-to-end training. In *International Conference on Learning Representations (ICLR 2021)*.

[Weinan, 2017] Weinan, E. (2017). A proposal on machine learning via dynamical systems. *Commun. Math. Stat.*

[Weng et al., 2018] Weng, T.-W., Zhang, H., Chen, P.-Y., Yi, J., Su, D., Gao, Y., Hsieh, C.-J., and Daniel, L. (2018). Evaluating the robustness of neural networks: An extreme value theory approach. In *International Conference on Learning Representations (ICLR 2018)*.

[Wightman et al., 2021] Wightman, R., Touvron, H., and Jégou, H. (2021). Resnet strikes back: An improved training procedure in timm. *arXiv*.

[Wu et al., 2020] Wu, D., Wang, Y., Xia, S.-T., Bailey, J., and Ma, X. (2020). Skip connections matter: On the transferability of adversarial examples generated with resnets. In *International Conference on Learning Representations (ICLR 2020)*.

[Wu et al., 2016] Wu, Y., Schuster, M., Chen, Z., Le, Q. V., Norouzi, M., Macherey, W., Krikun, M., Cao, Y., Gao, Q., Macherey, K., Klingner, J.,
Shah, A., Johnson, M., Liu, X., Łukasz Kaiser, Gouws, S., Kato, Y., Kudo, T., Kazawa, H., Stevens, K., Kurian, G., Patil, N., Wang, W., Young, C., Smith, J., Riesa, J., Rudnick, A., Vinyals, O., Corrado, G., Hughe, M., and Dean, J. (2016). Google's neural machine translation system: Bridging the gap between human and machine translation. *arXiv*.

[Xie et al., 2017] Xie, S. et al. (2017). Aggregated residual transformations for deep neural networks. In *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2017)*.

[Xiong et al., 2017] Xiong, W., Droppo, J., Huang, X., Seide, F., Seltzer, M., Stolcke, A., Yu, D., and Zweig, G. (2017). The microsoft 2016 conversational speech recognition system. In *IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2017)*.

[Xiong et al., 2020] Xiong, Y., Ren, M., and Urtasun, R. (2020). Loco: Local contrastive representation learning. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 11142-11153. Curran Associates, Inc.

[Xu and Principe, 1999] Xu, D. and Principe, J. (1999). Training mlps layer-by-layer with the information potential. In *International Joint Conference on Neural Networks (IJCNN'99)*, volume 3, pages 1716-1720 vol.3.

[Yan et al., 2020] Yan, H., Du, J., Tan, V., and Feng, J. (2020). On robustness of neural ordinary differential equations. *International Conference on Learning Representations (ICLR 2020)*.

[Yang et al., 2020a] Yang, Y., Wu, J., Li, H., Li, X., Shen, T., and Lin, Z. (2020a). Dynamical system inspired adaptive time stepping controller for residual network families. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 34, pages 6648-6655.

[Yang et al., 2020b] Yang, Y.-Y., Rashtchian, C., Zhang, H., Salakhutdinov, R. R., and Chaudhuri, K. (2020b). A closer look at accuracy vs. robustness. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 8588-8601. Curran Associates, Inc.

[Ypma, 1995] Ypma, T. J. (1995). Historical development of the newton-raphson method. *SIAM Review*, 37(4):531-551.

[Zagoruyko and Komodakis, 2016] Zagoruyko, S. and Komodakis, N. (2016). Wide residual networks. In *Proceedings of the British Machine Vision Conference (BMVC 2016)*.

[Zhang et al., 2017a] Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O. (2017a). Understanding deep learning requires rethinking generalization. In *International Conference on Learning Representations (ICLR 2017)*.
[Zhang et al., 2020] Zhang, C., Zhang, H., and Hsieh, C.-J. (2020). An efficient adversarial attack for tree ensembles. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, *Advances in Neural Information Processing Systems (NeurIPS 2020)*, volume 33, pages 16165–16176. Curran Associates, Inc.

[Zhang et al., 2019] Zhang, J. et al. (2019). Towards robust resnet: A small step but a giant leap. In *Proceedings of the Twenty-Eight International Joint Conference on Artificial Intelligence (IJCAI-19)*.

[Zhang and Wynter, 2018] Zhang, J. and Wynter, L. (2018). Smooth inter-layer propagation of stabilized neural networks for classification. *arXiv*.

[Zhang et al., 2022] Zhang, J., Zhang, X., Ma, K., Gupta, R., Saltz, J., Vakalopoulou, M., and Samaras, D. (2022). Gigapixel whole-slide images classification using locally supervised learning. In *Medical Image Computing and Computer Assisted Intervention – MICCAI 2022*, pages 192–201, Cham. Springer Nature Switzerland.

[Zhang et al., 2017b] Zhang, X., Li, Z., Loy, C., and Lin, D. (2017b). Polynet: A pursuit of structural diversity in very deep networks. In *2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2017)*, pages 3900–3908, Los Alamitos, CA, USA. IEEE Computer Society.

[Zhu et al., 2017] Zhu, J.-Y., Park, T., Isola, P., and Efros, A. A. (2017). Unpaired image-to-image translation using cycle-consistent adversarial networks. In *International Conference on Computer Vision (ICCV 2017)*.

[Zhu et al., 2022] Zhu, M., Chang, B., and Fu, C. (2022). Convolutional neural networks combined with runge-kutta methods. *Neural Computing and Applications*.

[Zhuang et al., 2021] Zhuang, H., Weng, Z., Luo, F., Kar-Ann, T., Li, H., and Lin, Z. (2021). Accumulated decoupled learning with gradient staleness mitigation for convolutional neural networks. In *Meila, M. and Zhang, T., editors, Proceedings of the 38th International Conference on Machine Learning (ICML 2021)*, volume 139 of *Proceedings of Machine Learning Research*, pages 12935–12944. PMLR.
<figure><img src="image_31.png" /></figure>
List of Acronyms

AA AutoAttack.

ADL Accumulated Decoupled Learning.

APGD Auto Projected Gradient Descent.

ATN Adversarial Transformation Networks.

AUROC Area Under Receiver Operating Curve.

BA Boundary Attack.

BFGS Broyden-Fletcher-Goldfarb-Shanno.

BIM Basic Iterative Method.

BN Batch-Normalization.

BS Bogacki-Shampine.

CKA Centered Kernel Alignment.

Conv Convolution.

CW Carlini-Wagner.

DDG Decoupled Delayed Gradients.

DF DeepFool.

DFT Discrete Fourier Transform.

DGL Decoupled Greedy Learning.

DNI Decoupled Neural Interface.

FA Feature Adversary.

FAB Fast Adaptive Boundary.

FGM Fast Gradient Method.

FPR False Positive Rate.

FR Features Replay.

GGD Generalized Gaussian Distribution.

HSJ HopSkipJump.

IMEX Implicit-Explicit.

InfoPro Information Propagation.
**JSMA** Jacobian Saliency Map Attack.

**KD** Kernel Density.

**LAP** Least Action Principle.

**LID** Local Intrinsic Dimensionality.

**LM** Linear Multi-step.

**LN** Layer Normalization.

**MFS** Magnitude Fourier Spectrum.

**MH** Mahalanobis.

**MLE** Maximum Likelihood Estimator.

**MLP** Multi-Layer Perception.

**MLS** Multi-Lap Sequential.

**MMD** Maximum Mean Discrepancy.

**MSA** Method of Successive Approximations.

**MSCN** Mean Subtracted Contrast Normalized.

**NODE** Neural Ordinary Differential Equation.

**NR** Newton-Raphson.

**NSS** Natural Scene Statistics.

**ODE** Ordinary Differential Equation.

**OOD** Out-Of-Distribution.

**OT** Optimal Transport.

**PDE** Partial Differential Equation.

**PFS** Phase Fourier Spectrum.

**PGD** Projected Gradient Descent.

**PMP** Pontryagin's Maximum Principle.

**PredSim** Prediction Similarity.

**RCE** Reverse Cross Entropy.

**ReLU** Rectified Linear Unit.

**ResBlock** Residual Block.

**ResNet** Residual Network.

**RK** Runge-Kutta.
**RMD** Relative Mahalanobis Distance.

**ROC** Receiver Operating Characteristic.

**SA** Square Attack.

**Sedona** Searching for Decoupled Neural Architectures.

**SGD** Stochastic Gradient Descent.

**ST** Spatial Transformation.

**SVCCA** Singular Vector Canonical Correlation Analysis.

**SW-MSA** Shifted Window Multi-head Self-Attention.

**TPR** True Positive Rate.

**TR** Transport.

**TRGL** Transport-Regularized Greedy Learning.

**UAP** Universal Adversarial Perturbation.

**VanGL** Vanilla Greedy Learning.

**W-MSA** Window Multi-head Self-Attention.

**ZOO** Zeroth Order Optimization.
<figure><img src="image_32.png" /></figure>
List of Tables

<table><tr><td>1</td><td>Accuracy of LAP-ResNet-9 on MNIST</td><td>31</td></tr><tr><td>2</td><td>Accuracy of LAP-ResNet-9 on CIFAR10</td><td>31</td></tr><tr><td>3</td><td>Accuracy of LAP-ResNeXt50 on CIFAR100</td><td>32</td></tr><tr><td>4</td><td>Accuracy of parallel TRGL with 4 modules on TinyImageNet</td><td>50</td></tr><tr><td>5</td><td>Accuracy of parallel TRGL with 2 modules on CIFAR100</td><td>50</td></tr><tr><td>6</td><td>Accuracy of parallel TRGL with K modules on STL10</td><td>51</td></tr><tr><td>7</td><td>Accuracy of parallel TRGL with 4 modules using a Swin Transformer</td><td>51</td></tr><tr><td>8</td><td>Accuracy of sequential and MLS TRGL on TinyImageNet</td><td>51</td></tr><tr><td>9</td><td>Memory savings of module-wise training on STL10</td><td>52</td></tr><tr><td>10</td><td>Accuracy of 10-1 ResNet with sequential, MLS and parallel TRGL</td><td>54</td></tr><tr><td>11</td><td>Last block accuracy of 10-1 ResNet with sequential TRGL</td><td>54</td></tr><tr><td>12</td><td>Detection accuracy of seen attacks</td><td>72</td></tr><tr><td>13</td><td>Detection accuracy of unseen attacks</td><td>73</td></tr><tr><td>14</td><td>OOD detection accuracy on ResNet110</td><td>74</td></tr><tr><td>15</td><td>OOD detection accuracy on ResNeXt50</td><td>74</td></tr><tr><td>16</td><td>Accuracy of LAP-ResNeXt50 on CIFAR100</td><td>159</td></tr><tr><td>17</td><td>Accuracy of LAP-ResNet-9 on CIFAR100</td><td>159</td></tr><tr><td>18</td><td>Accuracy of LAP-ResNet-9 on CIFAR10</td><td>159</td></tr><tr><td>19</td><td>Accuracy on the circles dataset with 1 block</td><td>162</td></tr><tr><td>20</td><td>Accuracy on the circles dataset with 100 blocks</td><td>162</td></tr><tr><td>21</td><td>Accuracy on the circles dataset with 50 training points</td><td>163</td></tr><tr><td>22</td><td>Accuracy on the circles dataset with a large initialization</td><td>163</td></tr><tr><td>23</td><td>Accuracy of 2-7 ResNet with parallel TRGL</td><td>169</td></tr><tr><td>24</td><td>Accuracy of parallel TRGL with 2 modules</td><td>170</td></tr><tr><td>25</td><td>Last block accuracy of 20-1 ResNet with parallel TRGL</td><td>170</td></tr><tr><td>26</td><td>Last block accuracy of 20-1 ResNet with parallel TRGL</td><td>170</td></tr><tr><td>27</td><td>Accuracy of 10-1 ResNet with sequential TRGL</td><td>171</td></tr><tr><td>28</td><td>Last block accuracy of 50-1 ResNet with sequential and MLS TRGL</td><td>171</td></tr><tr><td>29</td><td>Accuracy of ResNeXt50 with sequential, MLS and parallel TRGL</td><td>171</td></tr><tr><td>30</td><td>Accuracy of parallel TRGL* with K modules on STL10</td><td>172</td></tr><tr><td>31</td><td>Memory savings of module-wise training with 4 modules</td><td>172</td></tr><tr><td>32</td><td>Detection accuracy of seen attacks on ResNet110</td><td>175</td></tr><tr><td>33</td><td>Detection accuracy of seen attacks on ResNeXt50</td><td>176</td></tr><tr><td>34</td><td>Detection accuracy of seen attacks on WideResNet</td><td>176</td></tr><tr><td>35</td><td>Detection accuracy of unseen attacks on ResNet110</td><td>176</td></tr><tr><td>36</td><td>Detection accuracy of unseen attacks on ResNeXt50</td><td>177</td></tr><tr><td>37</td><td>Detection accuracy of unseen attacks on WideResNet</td><td>177</td></tr><tr><td>38</td><td>Detection rate of successful samples from seen attacks</td><td>177</td></tr><tr><td>39</td><td>Detection rate of successful samples from unseen attacks</td><td>178</td></tr><tr><td>40</td><td>False positive rate on seen attacks</td><td>178</td></tr><tr><td>41</td><td>False positive rate on unseen attacks</td><td>178</td></tr><tr><td>42</td><td>AUROC of seen attacks</td><td>179</td></tr><tr><td>43</td><td>AUROC of unseen attacks</td><td>179</td></tr></table>
<table><tr><td>44</td><td>AUROC of OOD detection on ResNet110</td><td>180</td></tr><tr><td>45</td><td>FPR at 95% TPR of OOD detection on ResNet110</td><td>180</td></tr><tr><td>46</td><td>AUROC of OOD detection on ResNeXt50</td><td>180</td></tr><tr><td>47</td><td>FPR at 95% TPR of OOD detection on ResNeXt50</td><td>181</td></tr></table>
List of Figures

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Basic residual block architecture</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Bottleneck residual block and ResNeXt block architectures</td>
      <td>7</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Bottleneck residual block and inverted residual block architectures</td>
      <td>7</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Consecutive basic residual blocks</td>
      <td>8</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Swin Transformer architecture with 4 residual stages</td>
      <td>8</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Ratio of residue norm to input norm and cosine loss of ResNets</td>
      <td>10</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Architectures inspired by numerical schemes</td>
      <td>14</td>
    </tr>
    <tr>
      <td>8</td>
      <td>The Kantorovitch problem from μ to ν</td>
      <td>16</td>
    </tr>
    <tr>
      <td>9</td>
      <td>The Monge problem from μ to ν</td>
      <td>16</td>
    </tr>
    <tr>
      <td>10</td>
      <td>Dynamic transport of μ<sub>0</sub> to μ<sub>1</sub> according to velocity field v</td>
      <td>17</td>
    </tr>
    <tr>
      <td>11</td>
      <td>Test transport cost against test accuracy of ResNet-9 models</td>
      <td>23</td>
    </tr>
    <tr>
      <td>12</td>
      <td>Circles test set after blocks 1, 5 and 9</td>
      <td>24</td>
    </tr>
    <tr>
      <td>13</td>
      <td>Decomposition of a neural network</td>
      <td>25</td>
    </tr>
    <tr>
      <td>14</td>
      <td>Accuracy of ResNets without batch-normalization</td>
      <td>33</td>
    </tr>
    <tr>
      <td>15</td>
      <td>Accuracy during training of ResNeXt50 models</td>
      <td>34</td>
    </tr>
    <tr>
      <td>16</td>
      <td>Module-wise training</td>
      <td>38</td>
    </tr>
    <tr>
      <td>17</td>
      <td>Early overfitting of vanilla module-wise training</td>
      <td>39</td>
    </tr>
    <tr>
      <td>18</td>
      <td>Accuracy of each module after module-wise training</td>
      <td>55</td>
    </tr>
    <tr>
      <td>19</td>
      <td>Example of an adversarial attack</td>
      <td>59</td>
    </tr>
    <tr>
      <td>20</td>
      <td>Boundary tilting perspective of adversarial samples</td>
      <td>61</td>
    </tr>
    <tr>
      <td>21</td>
      <td>Transport cost histogram for clean and attacked test samples</td>
      <td>69</td>
    </tr>
    <tr>
      <td>22</td>
      <td>Circles test set and OOD points after blocks 6 and 9</td>
      <td>69</td>
    </tr>
    <tr>
      <td>23</td>
      <td>Adversarial detection dataset creation</td>
      <td>70</td>
    </tr>
    <tr>
      <td>24</td>
      <td>The Monge problem from μ to ν</td>
      <td>125</td>
    </tr>
    <tr>
      <td>25</td>
      <td>The Kantorovitch problem from μ to ν</td>
      <td>126</td>
    </tr>
    <tr>
      <td>26</td>
      <td>Dynamic transport of μ<sub>0</sub> to μ<sub>1</sub> according to velocity field v</td>
      <td>136</td>
    </tr>
    <tr>
      <td>27</td>
      <td>Early overfitting of vanilla module-wise training</td>
      <td>145</td>
    </tr>
    <tr>
      <td>28</td>
      <td>Decoding the internal representations of a network.</td>
      <td>155</td>
    </tr>
    <tr>
      <td>29</td>
      <td>Internal representations of ResNet-9 on MNIST</td>
      <td>156</td>
    </tr>
    <tr>
      <td>30</td>
      <td>Internal representations of LAP-ResNet-9 on MNIST</td>
      <td>156</td>
    </tr>
    <tr>
      <td>31</td>
      <td>Internal representations of ResNet-9 on CIFAR10</td>
      <td>157</td>
    </tr>
    <tr>
      <td>32</td>
      <td>Internal representations of LAP-ResNet-9 on CIFAR10</td>
      <td>157</td>
    </tr>
    <tr>
      <td>33</td>
      <td>Internal representations of ResNet-9 on CIFAR100</td>
      <td>158</td>
    </tr>
    <tr>
      <td>34</td>
      <td>Internal representations of LAP-ResNet-9 on CIFAR100</td>
      <td>158</td>
    </tr>
    <tr>
      <td>35</td>
      <td>Weight distribution of the classifier of a vanilla ResNet-9</td>
      <td>160</td>
    </tr>
    <tr>
      <td>36</td>
      <td>Weight distribution of the classifier of a LAP-ResNet-9</td>
      <td>160</td>
    </tr>
    <tr>
      <td>37</td>
      <td>Distances between successive point clouds in a ResNet-9</td>
      <td>161</td>
    </tr>
    <tr>
      <td>38</td>
      <td>Circles test set embeddings with LAP</td>
      <td>164</td>
    </tr>
    <tr>
      <td>39</td>
      <td>Circles test set embeddings with large initialization</td>
      <td>165</td>
    </tr>
    <tr>
      <td>40</td>
      <td>Circles test set embeddings with large initialization and BN</td>
      <td>166</td>
    </tr>
    <tr>
      <td>41</td>
      <td>Circles test set embeddings with large initialization, BN and LAP</td>
      <td>167</td>
    </tr>
    <tr>
      <td>42</td>
      <td>Accuracy of each block after block-wise training</td>
      <td>171</td>
    </tr>
    <tr>
      <td>43</td>
      <td>Sensitivity of TRGL to hyperparameter τ</td>
      <td>173</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_33.png" /></figure>
# A Background on Numerical Methods for Differential Equations

We refer to [Quarteroni et al., 2007] for this quick background on numerical methods for solving ordinary differential equations.

Given an interval $I = [t_0, t_0 + T]$ and a function $f: I \times \mathbb{R}^d \to \mathbb{R}^d$, we consider the *initial value problem* of finding a function $x: I \to \mathbb{R}^d$ that satisfies

$$x'(t) = f(t, x(t))$$

with initial condition $x(t_0) = x_0$, or equivalently in integral form, a function $x$ that satisfies

$$x(t) = x_0 + \int_{t_0}^{t} f(t, x(t)) \, dt$$

It is often difficult to find closed-form solutions to initial value problems. We wish then to approximate the values $x$ takes at certain time points. We consider a subdivision $t_0 < t_1 < \dots < t_N = t_0 + T$ of $I$. Denote the *time-steps* (or *step-sizes*) $h_n := t_{n+1} - t_n$ for $0 \le n < N$ and define $h_{\max} := \max h_n$.

## A.1 The Euler Method

In a one-step method, an approximation of $x(t_n)$ is $x_n$ given iteratively starting from a known initial condition $x_0$ by

$$\frac{x_{n+1} - x_n}{h_n} = \phi(t_n, x_n, t_{n+1}, x_{n+1}, h_n)$$

If $f(x_{n+1})$ does not appear in the definition of $\phi$, the method is called explicit. Otherwise, it is called implicit, as $x_{n+1}$ depends then implicitly on itself through $f$ and an equation needs to be solved to find it, for example, using a fixed-point [Shashkin, 1991] or a Newton-Raphson [Ypma, 1995] method. So for $\phi(t_n, x_n, t_{n+1}, x_{n+1}, h_n) = f(t_n, x_n)$, we get the explicit Euler method:

$$x_{n+1} = x_n + h_n f(t_n, x_n)$$

And for $\phi(t_n, x_n, t_{n+1}, x_{n+1}, h_n) = f(t_{n+1}, x_{n+1})$, we get the implicit Euler method

$$x_{n+1} = x_n + h_n f(t_{n+1}, x_{n+1})$$

We only need the convergence analysis of explicit methods, so in the following we will limit ourselves to methods given by

$$\frac{x_{n+1} - x_n}{h_n} = \phi(t_n, x_n, h_n)$$
**Definition 4.** For a one-step method, the consistency errors $e_n$, for $0 \le n < N$, are

$$e_n := \frac{x(t_{n+1}) - \Phi(t_n, x(t_n), h_n)}{h_n} = \frac{x(t_{n+1}) - x(t_n)}{h_n} - \phi(t_n, x(t_n), h_n)$$

where $x$ is solution. The local (truncation) errors are $h_n e_n$. The method is consistent if $\max|e_n|$ goes to zero as $h_{\max}$ goes to zero. For $p \in \mathbb{N}^*$, the method has order $p$ if $\max|e_n| \le Ch_{\max}^p$ for a constant $C$ that depends on $f, t_0$ and $T$.

**Theorem 9.** If $f$ and $\phi$ are continuous, then the one-step method is consistent if and only if $\phi(t, x, 0) = f(t, x)$ for all $(t, x)$.

**Theorem 10.** If $f$ is $C^p$ and $\phi$ is $C^p$ in $h$ then the one-step method is of order $p$ if and only if $\partial_h^k \phi(t, x, 0) = \frac{1}{k+1} f^{[k]}(t, x)$ for all $(t, x)$ and $0 \le k < p$ where $f^{[0]} = f$ and $f^{[k]} = \partial_t f^{[k-1]} + f \partial_x f^{[k-1]}$.

**Corollary 11.** If $f$ is continuous then Euler's method is consistent. If $f$ is $C^1$ then Euler's method has order 1.

**Definition 5.** A one-step method is zero-stable (or stable) if $\exists S > 0$ such that for all $(x_n)_{0 \le n \le N}$, $(\tilde{x}_n)_{0 \le n \le N}$ and $(\varepsilon_n)_{0 \le n < N}$ satisfying

$$\frac{x_{n+1} - x_n}{h_n} = \phi(t_n, x_n, h_n) \text{ and } \frac{\tilde{x}_{n+1} - \tilde{x}_n}{h_n} = \phi(t_n, \tilde{x}_n, h_n) + \varepsilon_n \quad (44)$$

for $0 \le n < N$, we have $\max_n \|\tilde{x}_n - x_n\| \le S(\|\tilde{x}_0 - x_0\| + T \max_n |\varepsilon_n|)$, where $\varepsilon_n = \varepsilon_n/h_n$. The constant $S$ is the stability constant of the method.

**Theorem 12.** If $\phi$ is uniformly $L$-Lipschitz in its second variable, then the one-step method is stable with constant $e^{LT}$.

**Corollary 13.** If $f$ is Lipschitz in its second variable, then Euler's method is stable.

**Definition 6.** A numerical method converges if its global error $\max_n \|x(t_n) - x_n\|$ goes to zero as $h_{\max}$ goes to zero.

**Theorem 14.** If a method is consistent and stable with stability constant $S$, then it converges and $\max_n \|x(t_n) - x_n\| \le ST \max_n |e_n|$. If the method is of order $p$ with constant $C$, then $\max_n \|x(t_n) - x_n\| \le STCh_{\max}^p$

**Corollary 15.** Euler's method converges if $f$ is $C^0$ and Lipschitz in $x$. If $f$ is also $C^1$ then it converges with speed $O(h_{\max})$.

## A.2 Runge-Kutta Methods

**Definition 7.** Given $q \in \mathbb{N}^*$, a Runge-Kutta method with $q$ stages generates the estimates $(x_n)_{0 \le n \le N}$ through, for $0 \le n < N$,

$$t_{n,i} = t_n + c_i h_n, \quad x_{n,i} = x_n + h_n \sum_{j=1}^{i-1} a_{i,j} p_{n,j}, \quad p_{n,i} = f(t_{n,i}, x_{n,i}), \quad 0 \le i \le q$$
$$t_{n+1} = t_n + h_n, \quad x_{n+1} = x_n + h_n \sum_{i=1}^{q} b_i p_{n,i}$$

where $c_i, b_i, a_{i,j} \in [0, 1]$ for all $0 \le i \le q$ and $1 \le j < i$ with

$c_1 = 0$ so that $t_{n,1} = t_n, x_{n,1} = x_n$ and $p_{n,1} = f(t_n, x_n)$ for $0 \le n < N$

$$\sum_{j=1}^{i-1} a_{i,j} = c_i \quad \text{for } 0 \le i \le q \quad \text{and} \quad \sum_{i=1}^{q} b_i = 1$$

This is an explicit one-step method with $\phi$ given by

$$\left\{
\begin{array}{l}
\phi(t, x, h) = \sum_{i=1}^{q} b_i f(t + c_i h, x_i) \\
\\
x_i = x + h \sum_{j=1}^{i-1} a_{i,j} f(t + c_j h, x_j) \quad \text{for } 1 \le i \le q
\end{array}
\right.
\qquad (45)$$

**Remark 16.** The conditions on $(a_{i,j})$ and $(b_i)$ are so that the numerical integration methods used are of order 1, i.e. exact on constant functions, which is equivalent to the coefficients summing up to the length of the integration interval.

**Remark 17.** For $q = 1$ we find Euler's method.

**Theorem 18.** If $f$ is continuous then Runge-Kutta methods are consistent.

**Definition 8.** The Runge-Kutta method with $q = 2$ stages, $c_2 = a_{2,1} = 1$ and $b_1 = b_2 = 1/2$ is called Heun's method and is given by $\phi(t, x, h) = \frac{1}{2}f(t, x) + \frac{1}{2}f(t+h, x+hf(t, x))$.

**Proposition 4.** If $f$ is continuous then Heun's method is consistent. If $f$ is $C^2$ then Heun's method has order 2.

**Definition 9.** The Runge-Kutta method with $q = 4$ stages, $c_2 = c_3 = a_{2,1} = a_{3,2} = 1/2$, $c_4 = a_{4,3} = 1$, $b_1 = b_4 = 1/6$, $b_2 = b_3 = 1/3$ and all other coefficients equal to zero is the classic RK4 method. It is given by

$$\left\{
\begin{array}{l}
t_{n,1} = t_n \\
t_{n,2} = t_n + \frac{h_n}{2} \\
t_{n,3} = t_n + \frac{h_n}{2} \\
t_{n,4} = t_n + h_n
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
x_{n,1} = x_n \\
x_{n,2} = x_n + \frac{h_n}{2} p_{n,1} \\
x_{n,3} = x_n + \frac{h_n}{2} p_{n,2} \\
x_{n,4} = x_n + h_n p_{n,3}
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
p_{n,1} = f(t_{n,1}, x_{n,1}) \\
p_{n,2} = f(t_{n,2}, x_{n,2}) \\
p_{n,3} = f(t_{n,3}, x_{n,3}) \\
p_{n,4} = f(t_{n,4}, x_{n,4})
\end{array}
\right.$$

And finally $t_{n+1} = t_n + h_n$ and $x_{n+1} = x_n + \frac{h_n}{6}(p_{n,1} + 2p_{n,2} + 2p_{n,3} + p_{n,4})$

**Proposition 5.** If $f$ is continuous then the RK4 method is consistent. If $f$ is $C^4$ then the RK4 method has order 4.
**Theorem 19.** If $f$ is continuous and $L$-Lipschitz in its second variable and the time-steps are bounded then a Runge-Kutta method with $q$ stages is stable with stability constant $S = e^{KT}$ where

$$K = L \sum_{i=1}^{q} b_i \sum_{j=0}^{i-1} L^j h_{\max}^j$$

**Corollary 20.** Runge-Kutta methods converge if $f$ is $C^0$ and Lipschitz in $x$.

**Corollary 21.** Heun's method converges if $f$ is $C^0$ and Lipschitz in $x$. If $f$ is also $C^2$ then it converges with speed $O(h_{\max}^2)$.

**Corollary 22.** The RK4 method converges if $f$ is $C^0$ and Lipschitz in $x$. If $f$ is also $C^4$ then it converges with speed $O(h_{\max}^4)$.

## A.3 Adaptive Runge-Kutta Methods

Adaptive methods [Söderlind, 2006] estimate the local truncation error of a single Runge-Kutta step and decrease the step size if the error estimate is large (which indicates that the differential equation is not approximated well at that point and more conservative steps need to be taken) and increase it if it is low. This is done by using two interwoven Runge-Kutta methods that have common intermediate stages, meaning that their coefficients $p_{n,i}$ in Definition 7 are the same but the weights $b_i$ in front of them differ. The local error is estimated by the difference between the predictions of the two methods. This way, estimating the error has little or negligible computational cost compared to a step of the higher-order method which already requires the computation of the lower-order method. For example, Heun's method includes Euler's method.

During the integration, the step size is adapted such that the estimated error stays below a chosen threshold. If the estimated error is too high, a step is repeated with a lower step size. If the estimated error is much smaller than the threshold, the step size is increased to save time. This results in an (almost) optimal step size, which saves computation time.

For example, a step of the adaptive Bogacki-Shampine adaptive Runge-Kutta method starting from $x_n$ at time $t_n$ is computed as follows. The gradient evaluation are given by

$$\left\{
\begin{aligned}
p_{n,1} &= f(t_n, x_n) \\
p_{n,2} &= f(t_n + \frac{1}{2}h_n, x_n + \frac{1}{2}h_n p_{n,1}) \\
p_{n,3} &= f(t_n + \frac{3}{4}h_n, x_n + \frac{3}{4}h_n p_{n,2}) \\
p_{n,4} &= f(t_n + h_n, x_n + \frac{2}{9}h_n p_{n,1} + \frac{1}{3}h_n p_{n,2} + \frac{4}{9}h_n p_{n,3})
\end{aligned}
\right.
\qquad (46)$$
And the predictions for $x(t_{n+1})$ are given by

$$ 
\left\{
\begin{aligned}
x_{n+1}^{\text{euler}} &= x_n + h_n p_{n,1} \\
x_{n+1}^{\text{midpoint}} &= x_n + h_n p_{n,2} \\
x_{n+1}^{\text{order3}} &= x_n + \frac{2}{9} h_n p_{n,1} + \frac{1}{3} h_n p_{n,2} + \frac{4}{9} h_n p_{n,3} \\
x_{n+1}^{\text{order2}} &= x_n + \frac{7}{24} h_n p_{n,1} + \frac{1}{4} h_n p_{n,2} + \frac{1}{3} h_n p_{n,3} + \frac{1}{8} h_n p_{n,4}
\end{aligned}
\right.
\qquad (47) $$

The difference between $x_{n+1}^{\text{order3}}$ and $x_{n+1}^{\text{order2}}$ is used to adapt the step-size. Different rules exist. They all require the user to fix a tolerance $\varepsilon > 0$. The simplest [Quarteroni et al., 2007] is to simply halve the step-size if the error $\text{err} = \|x_{n+1}^{\text{order3}} - x_{n+1}^{\text{order2}}\|$ is larger than $\varepsilon$, to double it if the error is smaller than $\varepsilon/K$, and to keep it as is if the error is between $\varepsilon$ and $\varepsilon/K$. Usually $K = 2^{p+1}$, where $p$ is the order of the method. Other rules [Söderlind, 2006] are

$$ h_{n+1} = \left(\frac{\varepsilon}{\text{err}}\right)^p h_n $$

and

$$ \text{err} = \sqrt{\frac{1}{d} \sum_{i=1}^{d} \left( \frac{x_{n+1,i}^{(\text{Order3})} - x_{n+1}^{(\text{Order2})}}{\text{sc}_i} \right)^2} $$

where

$$ \text{sc}_i = \varepsilon(1 + \max(|x_{n+1,i}^{(\text{Order3})}|, |x_{n+1,i}^{(\text{Order2})}|)) $$

and then

$$ h_{\text{new}} = h_n \min(\text{facmax}, \max(\text{facmin}, \alpha(\frac{1}{\text{err}})^{\frac{1}{p}})) $$

where usually $\text{facmax} \in [1.5, 5]$ and $\alpha \in \{0.8, 0.9, 0.25^{1/3}, 0.38^{1/3}\}$. Then, if $\text{err} \le 1$, the step-size $h_n$ is kept as is and $h_{n+1}$ is set to $h_{\text{new}}$. And if $\text{err} > 1$, then $h_n$ is set to $h_{\text{new}}$ and the step is redone.
<figure><img src="image_34.png" /></figure>
# B Background on Optimal Transport Theory

## B.1 Monge and Kantorovich Problems

We follow [Santambrogio, 2015, Villani, 2008, Villani, 2021] in this overview of optimal transport. Denote $\mathcal{P}(E)$ the space of probability measures on a space $E$. Given probability measures $\mu \in \mathcal{P}(X)$ and $\nu \in \mathcal{P}(Y)$ and a cost function $c: X \times Y \to [0, \infty]$, the starting point of optimal transport is the following problem

**Definition 10.** The Monge problem is

$$ \begin{array}{ll} \inf_{T: X \to Y} & \displaystyle\int_X c(x, T(x)) \, d\mu(x) \\ \text{subject to} & T_{\#} \mu = \nu \end{array} \qquad (48) $$

Here $T_{\#}\mu$ is the push-forward measure of $\mu$ through $T$ defined by the formula $T_{\#}\mu(A) := \mu(T^{-1}(A)) \forall A$ measurable. For this definition and the constraint to make sense, $T$ has to be measurable with respect to the $\sigma$-algebras of $\mu$ and $\nu$. The Monge problem expresses the problem of moving mass distributed according to $\mu$ so that it becomes distributed according to $\nu$, while minimizing the sum (so the integral) of efforts made to transport each point, the effort being encoded in the cost function $c$, which is usually a distance between points.

An admissible $T$ is called a transport map between $\mu$ and $\nu$ and the objective of the problem is called the transport cost of $T$.

<figure><img src="image_35.png" /><figcaption>Figure 24: The Monge problem from $\mu$ to $\nu$. From [Kolouri et al., 2017].</figcaption></figure>

The Monge problem is difficult because of the constraint on $T$ which is not closed under weak convergence and a $T$ satisfying it need not even exist (e.g. when $\mu$ is a Dirac measure and $\nu$ is not). A relaxation exists
**Definition 11.** The Kantorovitch problem is

$$ \inf_{\gamma \in \Pi(\mu, \nu)} \int_{X \times Y} c \, d\gamma \qquad (49) $$

where $\Pi(\mu, \nu) = \{\gamma \in \mathcal{P}(X \times Y) | (\pi_x)_{\#}\gamma = \mu, (\pi_y)_{\#}\gamma = \nu\}$, $\pi_x$ is the projection from $X \times Y$ on $X$ and $\pi_y$ is the projection from $X \times Y$ on $Y$.

The constraint set $\Pi(\mu, \nu)$ of the Kantorovich problem is never empty as it always contains $\mu \otimes \nu$. This relaxation means that instead of moving the entire mass at $x$ to $T(x)$, $\gamma(A \times B)$ denotes the amount of mass moving from $A$ to $B$, so that two particles located at $x \in X$ can move to different destination points in $Y$ (see Figure 25 below). An admissible $\gamma \in \Pi(\mu, \nu)$ is called a transport plan between $\mu$ and $\nu$.

<figure><img src="image_36.png" /><figcaption>Figure 25: The Kantorovitch problem from $\mu$ to $\nu$. From [Peyre and Cuturi, 2019].</figcaption></figure>

The Kantorovich problem has a smaller infimum than the Monge problem, since if $T_{\#}\mu = \nu$ then $\gamma_T := (\mathrm{id}, T)_{\#}\mu \in \Pi(\mu, \nu)$ and the transport costs are equal

$$ \int_{X \times Y} c \, d\gamma_T = \int_X c(x, T(x)) \, d\mu(x) \qquad (50) $$

If $\gamma \in \Pi(\mu, \nu)$ is concentrated on a graph in $X \times Y$, i.e. on a set of the form $\{(x, T(x)), x \in X\}$, then all the mass in $x$ goes to $T(x)$ and $\gamma = (\mathrm{id}, T)_{\#}\mu$.

The Kantorovich problem is a relaxation of the Monge problem in the following sense: denote the Kantorovich objective $J(\gamma) := \int c \, d\gamma$ for $\gamma \in \Pi(\mu, \nu)$. The Monge objective can also be written as a functional $K$ on $\Pi(\mu, \nu)$ defined by $K(\gamma) := J(\gamma)$ if $\gamma = \gamma_T$ for a transport map $T$ and $K(\gamma) := \infty$ otherwise. The Monge problem is equivalent to minimizing $K$ over $\Pi(\mu, \nu)$ because of
(50). Then under some conditions the Kantorovich cost $J$ is the largest lower semi-continuous functional on $\Pi(\mu, \nu)$ that is smaller than the Monge cost $K$.

For example, if $X = Y = \Omega \subset \mathbb{R}^d$ with $\Omega$ compact and $c$ is continuous and $\mu$ atomless, then problems (48) and (49) have the same infimum and the set $\{\gamma_T = (\mathrm{id}, T)_{\#}\mu \mid T_{\#}\mu = \nu\}$ is dense (in a weak-$^*$ sense) in $\Pi(\mu, \nu)$ so any transport plan can be approximated by transport maps. This allows to prove that $J$ is a relaxation of $K$.

We have general existence results for the Kantorovich problem. We will prove the following one.

**Theorem 23.** If $X$ and $Y$ are separable compact metric spaces and $c$ is continuous (or lower semi-continuous), then the Kantorovich problem admits a solution.

*Proof.* The proof uses the direct method of calculus of variations. First, we know by Riesz representation theorem that, for a separable locally compact metric space $E$, the space of finite signed Borel measures on $E$ is the dual of $C_0(E)$ (the set of continuous real valued functions on $E$ that vanish at infinity, which is equal to the set of continous bounded functions $C_b(E)$ and to the set of continuous functions $C(E)$ when $E$ is compact). The norm of an element of the dual is the total mass of the measure. $\Pi(\mu, \nu)$ is then included in the dual of $C(X \times Y)$ and is bounded since all its elements have mass equal to 1. By the Banach-Alaoglu theorem the unit ball of the dual space is weak-$^*$ compact in duality with $C(X \times Y)$.

Denote $J(\gamma) := \int_{X \times Y} c \, d\gamma$ and $l$ the value of the infimum in (49). Take a minimizing sequence $(\gamma_n)_n$ of problem (49), i.e. a sequence such that $J(\gamma_n) \to l$ and $\gamma_n \in \Pi(\mu, \nu)$. By compacity, we can assume up to sub-sequencing that $\gamma_n \xrightarrow{\ast} \gamma$. This implies for $c$ continuous that $J(\gamma_n) \to J(\gamma)$ and therefore that $J(\gamma) = l$. It remains to show that $\gamma \in \Pi(\mu, \nu)$. First, $\gamma$ is clearly a probability measure since by taking the constant function equal to 1 everywhere we have $\gamma(X \times Y) = \int_{X \times Y} 1 \, d\gamma = \lim_{X \times Y} \int 1 \, d\gamma_n = 1$. Then take $\varphi \in C(X)$. We have $\int_{X \times Y} \varphi(x) \, d\gamma(x, y) = \lim_{X \times Y} \int \varphi(x) \, d\gamma_n(x, y) = \int \varphi \, d\mu$, because $\gamma_n$ is in $\Pi(\mu, \nu)$. Likewise for $\psi \in C(Y)$, $\int_{X \times Y} \psi(y) \, d\gamma(x, y) = \int \psi \, d\nu$.

We used that $(\pi_x)_{\#}\gamma = \mu$ is equivalent to $\int_{X \times Y} \varphi(x) \, d\gamma(x, y) = \int \varphi \, d\mu$ for all test functions $\varphi$ and that $\gamma_n \xrightarrow{\ast} \gamma$ by definition means $\int_{X \times Y} f \, d\gamma_n \to \int_{X \times Y} f \, d\gamma$ for all $f \in C(X \times Y)$. Finally $\gamma \in \Pi(\mu, \nu)$ and so $\gamma$ is optimal. $\square$

Existence for the Kantorovich problem can be extended to the case of complete and separable metric spaces (Polish spaces) $X$ and $Y$ via Prokhorov's theorem. For the Monge problem, we have the following result which we will need later.

**Theorem 24.** Consider $X = Y = \Omega \subset \mathbb{R}^d$ with $\Omega$ compact and $c(x, y) = h(x-y)$ with $h$ strictly convex. If $c$ is continuous we have the existence of an optimal transport plan $\gamma$. If $\mu$ is absolutely continuous with respect to the Lebesgue
measure and $\partial\Omega$ is $\mu$-negligible, then the optimal transport plan $\gamma$ is unique and is of the form $(\text{id}, T)_{\#}\mu$ for a transport map $T$, meaning that it is constructed from a unique optimal transport map $T$ for the Monge problem.

A final result that we mention is the following, which says that the inverse of the optimal transport map between $\mu$ and $\nu$ is the optimal transport map from $\nu$ to $\mu$,

**Theorem 25.** If $\mu$ and $\nu$ are absolutely continuous measures supported respectively on compact subsets $X$ and $Y$ of $\mathbb{R}^d$ with negligible boundaries, and if $c(x, y) = h(x - y)$ with $h$ strictly convex, then there exists a unique couple $(T, S)$ of functions such that the five following points hold

$$-T: X \to Y \text{ and } S: Y \to X$$

$$-T_{\#}\mu = \nu \text{ and } S_{\#}\nu = \mu$$

$$-T \text{ is optimal for the Monge problem from } \mu \text{ to } \nu \text{ with the cost } c$$

$$-S \text{ is optimal for the Monge problem from } \nu \text{ to } \mu \text{ with the cost } c$$

$$-T \circ S \stackrel{\nu-\text{a.s.}}{=} \text{id} \text{ and } S \circ T \stackrel{\mu-\text{a.s.}}{=} \text{id}$$

## B.2 Duality

The Kantorovitch problem is a linear optimization problem with convex constraints. Duality theory was hence used to study it. We simply state here the duality result, which will be needed to differentiate Wasserstein distances later.

The first step is to rewrite the constraint $\gamma \in \Pi(\mu, \nu)$ as follows

$$\sup_{\varphi, \psi} \int_X \varphi \, d\mu + \int_Y \psi \, d\nu - \int_{X \times Y} \varphi(x) + \psi(y) \, d\gamma(x, y)$$

as this quantity is equal to 0 if $\gamma \in \Pi(\mu, \nu)$ and to $\infty$ otherwise. The Kantorovitch problem can then be rewritten

$$\inf_{\gamma} \int_{X \times Y} c \, d\gamma + \sup_{\varphi, \psi} \int_X \varphi \, d\mu + \int_Y \psi \, d\nu - \int_{X \times Y} (\varphi(x) + \psi(y)) \, d\gamma(x, y) \quad (51)$$

Interchanging the infimum and supremum in the above expression, we get

$$\sup_{\varphi, \psi} \int_X \varphi \, d\mu + \int_Y \psi \, d\nu + \inf_{\gamma} \int_{X \times Y} (c(x, y) - \varphi(x) - \psi(y)) \, d\gamma(x, y) \quad (52)$$

The dual of the Kantorovitch problem comes from showing that the maximum of (52) is equal to the minimum of the Kantorovitch problem, and from rewriting $\inf_{\gamma} \int_{X \times Y} (c(x, y) - \varphi(x) - \psi(y)) \, d\gamma(x, y)$ as the condition $\varphi(x) + \psi(y) \le c(x, y)$ first, and then as a property known as *c-concavity*.
**Definition 12.** Given $\varphi: X \to \bar{\mathbb{R}}$, its $c$-transform (c-conjugate) is $\varphi^c: Y \to \bar{\mathbb{R}}$ given by

$$ \varphi^c(y) = \inf_{x \in X} c(x, y) - \varphi(x) $$

Similarly, the $\bar{c}$-transform of $\psi: Y \to \bar{\mathbb{R}}$ is $\psi^c: X \to \bar{\mathbb{R}}$ given by

$$ \psi^{\bar{c}}(x) = \inf_{y \in Y} c(x, y) - \psi(x) $$

**Definition 13.** A function $\varphi: X \to \bar{\mathbb{R}}$ is said to be $c$-concave if there exists $\psi: Y \to \bar{\mathbb{R}}$ such that $\varphi = \psi^{\bar{c}}$. We denote $c$-conc($X$) the set of $c$-concave functions on $X$.

**Theorem 26.** If $X$ and $Y$ are separable compact metric spaces and $c$ is uniformly continuous and bounded, then

$$ \min_{\gamma \in \Pi(\mu, \nu)} \int_{X \times Y} c \, d\gamma = \max_{\varphi \in c\text{-conc}(X)} \int_X \varphi \, d\mu + \int_Y \varphi^c \, d\nu $$

The functions $\varphi$ realizing the maximum in the above equation are called the *Kantorovitch potentials* for the transport from $\mu$ to $\nu$.

**Theorem 27.** Under the conditions of Theorem 24, the optimal transport map between $\mu$ and $\nu$ is given by $T(x) = x - (\nabla h)^{-1}(\nabla \varphi(x))$ where $\varphi$ is a Kantorovitch potential for the transport from $\mu$ to $\nu$.

## B.3 The Wasserstein Space

From now on, we take $X = Y = \Omega \subset \mathbb{R}^d$, with $\Omega$ complete (i.e. closed), and $c(x, y) = \|x - y\|^p$ for $p \in [1, \infty[$ so that an optimal transport plan for the Kantorovich problem exists (by the extension of Theorem 23 to Polish spaces). Note that if $p > 1$ then $c$ is of the form $c(x, y) = h(x - y)$ with $h$ strictly convex and Theorem 24 stands if $\Omega$ is compact and $\mu$ absolutely continuous with respect to the Lebesgue measure. Compactness is not necessary if $c(x, y) = \|x - y\|^2$.

A distance between probability measures can be deduced from the Kantorovich problem. For $p \in [1, \infty[$, we restrict ourselves to the set

$$ \mathcal{P}_p(\Omega) := \{\mu \in \mathcal{P}(\Omega) \mid \int \|x\|^p d\mu(x) < \infty\} \quad (53) $$

where $\|\cdot\|$ is the euclidean norm. For $\mu, \nu \in \mathcal{P}_p(\Omega)$ we define $W_p(\mu, \nu)$ to be the $p$-th root of the minimal value of the Kantorovich problem between $\mu$ and $\nu$ for the cost $c(x, y) = \|x - y\|^p$, i.e.

$$ W_p^p(\mu, \nu) := \min_{\gamma \in \Pi(\mu, \nu)} \int_{\Omega \times \Omega} \|x - y\|^p d\gamma(x, y) \quad (54) $$

It can be shown that $W_p$ is a distance over $\mathcal{P}_p(\Omega)$ called the *Wasserstein distance*. One way to prove this is through the *gluing lemma*:
**Lemma 1.** Let $(X_1, \mu_1)$, $(X_2, \mu_2)$ and $(X_3, \mu_3)$ be Polish probability spaces. Take $\gamma_1 \in \Pi(\mu_1, \mu_2) \subset \mathcal{P}(X_1 \times X_2)$ and $\gamma_2 \in \Pi(\mu_2, \mu_3) \subset \mathcal{P}(X_2 \times X_3)$. Then there exists $\gamma \in \mathcal{P}(X_1 \times X_2 \times X_3)$ such that $(\pi_{1,2})_{\#} \gamma = \gamma_1$ and $(\pi_{2,3})_{\#} \gamma = \gamma_2$, where $\pi_{i,j}$ is the projection from $X_1 \times X_2 \times X_3$ on $X_i \times X_j$.

**Proposition 6.** $W_p$ is a distance over $\mathcal{P}_p(\Omega)$.

*Proof.* $W_p$ is clearly positive and finite because $c$ is positive and finite, the constraint set is never empty and because of the restriction to $\mathcal{P}_p(\Omega)$. $W_p$ is also clearly symmetric: define $F(x, y) := (y, x)$ then $F_{\#} \gamma \in \Pi(\nu, \mu)$ if and only if $\gamma \in \Pi(\mu, \nu)$ and the transport costs are the same because $c$ is symmetric.

Taking $\gamma = (\mathrm{id}, \mathrm{id})_{\#} \mu \in \Pi(\mu, \mu)$ gives $W_p(\mu, \mu) = \int c(x, x) \, d\mu(x) = 0$. Conversely, $W_p(\mu, \nu) = 0$ implies the existence of $\gamma \in \Pi(\mu, \nu)$ such that $\int c \, d\gamma = 0$. This means $\gamma$ has to be concentrated on $c^{-1}(\{0\}) = \{(x, y) \mid x = y\}$. This means that $\gamma = (\mathrm{id}, \mathrm{id})_{\#} \mu$ and that $\nu = \mathrm{id}_{\#} \mu = \mu$.

To prove the triangle inequality, consider $\mu_1, \mu_2$ and $\mu_3 \in \mathcal{P}_p(\Omega)$. There exists $\gamma_1 \in \Pi(\mu_1, \mu_2)$ an optimal transport plan between $\mu_1$ and $\mu_2$, and $\gamma_2 \in \Pi(\mu_2, \mu_3)$ an optimal transport plan between $\mu_2$ and $\mu_3$. Then by Lemma 1, there exists $\gamma \in \mathcal{P}(\Omega^3)$ with $(\pi_{1,2})_{\#} \gamma = \gamma_1$ and $(\pi_{2,3})_{\#} \gamma = \gamma_2$. In particular, this means that $(\pi_{1,3})_{\#} \gamma \in \Pi(\mu_1, \mu_3)$ because

$$ (\pi_{1,3})_{\#} \gamma(A \times \Omega) = (\pi_1)_{\#} \gamma(A) = (\pi_1)_{\#} ((\pi_{1,2})_{\#} \gamma)(A) = (\pi_1)_{\#} \gamma_1(A) = \mu_1(A) $$

where $\pi_1$ denotes the projection on the first coordinate (from $\Omega^3$ or $\Omega^2$ depending on the context) and we used the general formulas $(f \circ g)_{\#} \gamma = f_{\#}(g_{\#} \gamma)$ and $\pi_1 = \pi_1 \circ \pi_{1,2}$.

Likewise $(\pi_{1,3})_{\#} \gamma(\Omega \times B) = \mu_3(B)$. Since $(\pi_{1,3})_{\#} \gamma \in \Pi(\mu_1, \mu_3)$ but is not necessarily optimal we can write (denoting $d(x, y) := \|x - y\|$)

$$ \begin{aligned} W_p(\mu_1, \mu_3) &\le \left( \int d(x, z)^p \, d(\pi_{1,3})_{\#} \gamma(x, z) \right)^{\frac{1}{p}} \\ &= \left( \int d(x, z)^p \, d\gamma(x, y, z) \right)^{\frac{1}{p}} \\ &\le \left( \int (d(x, y) + d(y, z))^p \, d\gamma(x, y, z) \right)^{\frac{1}{p}} \\ &\le \left( \int d(x, y)^p \, d\gamma(x, y, z) \right)^{\frac{1}{p}} + \left( \int d(y, z)^p \, d\gamma(x, y, z) \right)^{\frac{1}{p}} \\ &= \left( \int d(x, y)^p \, d\gamma_1(x, y) \right)^{\frac{1}{p}} + \left( \int d(y, z)^p \, d\gamma_2(y, z) \right)^{\frac{1}{p}} \\ &= W_p(\mu_1, \mu_2) + W_p(\mu_2, \mu_3) \end{aligned} $$

where the last equality comes from the optimality of $\gamma_1$ and $\gamma_2$ and the third inequality comes from the Minkowski inequality $\|f + g\|_p \le \|f\|_p + \|g\|_p$. $\square$
We denote $W_p(\Omega)$ the metric space $P_p(\Omega)$ endowed with the distance $W_p$, and call it the *Wasserstein space*.

**Definition 14.** A curve $w$ in a metric space $(E, d)$ is a continuous function $w: [0, 1] \to E$. Its *metric derivative* at time $t$ is

$$|w'(t)| := \lim_{h \to 0} \frac{d(w(t+h), w(t))}{|h|} \qquad (55)$$

provided the limit exists.

**Definition 15.** A curve $w$ is *absolutely continuous* if there exists $g \in L^1([0, 1])$ such that

$$d(w(t_0), w(t_1)) \le \int_{t_0}^{t_1} g(s) \, ds \quad \forall t_0 < t_1 \in [0, 1] \qquad (56)$$

We denote $\text{AC}(E)$ the set of absolutely continuous curves valued in $E$. From the bound in (56) we see that absolutely continuous curves admit a metric derivative at all times.

**Definition 16.** The *length* of a curve $w$ is

$$\text{Length}(w) := \sup \left\{ \sum_{k=0}^{n-1} d(w(t_k), w(t_{k+1})) \mid 0 = t_0 < \dots < t_n = 1, n \ge 1 \right\} \qquad (57)$$

From Definition 15 of absolutely continuous curves and Definition 16 of the length of curve, we can see that every absolutely continuous curve $w$ has finite length, as $\text{Length}(w) \le \int_0^1 g(s) \, ds < \infty$ for some $g \in L^1([0, 1])$. Actually, we have a formula for the length of an absolutely continuous curve.

**Proposition 7.** For $w \in \text{AC}(E)$

$$\text{Length}(w) = \int_{0}^{1} |w'(t)| \, dt \qquad (58)$$

It is immediate from the definition of length and the triangle inequality that for any $x, y \in E$ we have $d(x, y) \le \text{Length}(w)$ for every curve $w$ connecting $x$ to $y$.

**Definition 17.** A metric space $(E, d)$ is called a *geodesic space* if $\forall x, y \in E$

$$d(x, y) = \min \left\{ \text{Length}(w) \mid w \in \text{AC}(E), w(0) = x, w(1) = y \right\} \qquad (59)$$

which means that an absolutely continuous curve must realize the minimum and that its length must be equal to $d(x, y)$. Such a curve is called a *geodesic* between $x$ and $y$.

**Definition 18.** A curve $w$ in $(E, d)$ is called a *constant-speed geodesic* between $w(0)$ and $w(1)$ if for all $t, s \in [0, 1]$

$$d(w(t), w(s)) = |t - s| \, d(w(0), w(1)) \qquad (60)$$
If we take $g(t) = d(w(0), w(1))$ for all $t \in [0, 1]$, it is immediate that a constant-speed geodesic $w$ is absolutely continuous. Its metric derivative is clearly equal to $d(w(0), w(1))$ at all times. This implies via Proposition 7 that its length is $d(w(0), w(1))$ and therefore that it is a geodesic between its endpoints. In fact, we have two other equivalent definitions of constant-speed geodesics.

**Proposition 8.** Fix $p > 1$ and consider a curve $w$ connecting $x$ to $y$ in a metric space $(E, d)$, then the following three points are equivalent

- $w$ is a constant-speed geodesic

- $w$ is absolutely continuous and $|w'|(t) = d(w(0), w(1))$ for almost every $t$

- $w$ solves $\min\{\int_0^1 |v'|^p(t) dt \mid v \in AC(E), v(0) = x, v(1) = y\}$

Now we prove that $W_p(\Omega)$ is a geodesic space when $\Omega$ is convex.

**Theorem 28.** $W_p(\Omega)$ for $\Omega$ convex is a geodesic space, i.e. there exists a geodesic curve between every two points in $W_p(\Omega)$.

More precisely, for $t \in [0, 1]$, define $\pi_t: \Omega \times \Omega \to \Omega$ through $\pi_t(x, y) = (1-t)x+ty$. For any $\mu, \nu \in W_p(\Omega)$, there exists by the previous section an optimal transport plan $\gamma \in \Pi(\mu, \nu)$ between $\mu$ and $\nu$ for the cost $c(x, y) = \|x - y\|^p$, i.e. $\gamma$ is optimal for the problem in (54). The curve defined by $\mu_t := (\pi_t)_{\#}\gamma$ for $t \in [0, 1]$ is then a constant-speed geodesic between $\mu$ and $\nu$.

*Proof.* Take $\mu, \nu \in P_p(\Omega)$ and $\gamma \in \Pi(\mu, \nu)$ optimal for problem (54). Define the curve $(\mu_t)_t$ through $\mu_t := (\pi_t)_{\#}\gamma$. It is immediate that $\mu_0 = \mu$ and $\mu_1 = \nu$. We start by proving that for every $t < s$, $W_p(\mu_t, \mu_s) \le |t-s| W_p(\mu, \nu)$. First, take $\gamma_t^s := (\pi_t, \pi_s)_{\#}\gamma$. Then $\gamma_t^s \in \Pi(\mu_t, \mu_s)$ because for every measurable set $A \subset \Omega$ we have

$$ \gamma_t^s(\Omega \times A) = \gamma((\pi_t, \pi_s)^{-1}(\Omega \times A)) = \gamma(\pi_s^{-1}(A)) = ((\pi_s)_{\#}\gamma)(A) = \mu_s(A) $$

and likewise $\gamma_t^s(B \times \Omega) = \mu_t(B)$. Since $\gamma_t^s$ is an admissible transport plan between $\mu_t$ and $\mu_s$ but not necessarily an optimal one we can write

$$ \begin{align*} W_p(\mu_t, \mu_s) &\le \left( \int \|x - y\|^p d\gamma_t^s(x, y) \right)^{\frac{1}{p}} \\ &= \left( \int \|\pi_t(x, y) - \pi_s(x, y)\|^p d\gamma(x, y) \right)^{\frac{1}{p}} \\ &= |t - s| \left( \int \|x - y\|^p d\gamma(x, y) \right)^{\frac{1}{p}} \\ &= |t - s| W_p(\mu, \nu) \end{align*} $$

The first equality comes from the definition of $\gamma_t^s$ as a push-forward of $\gamma$, i.e from the formula $\int F(x, y) d\gamma_t^s(x, y) = \int F(\pi_t(x, y), \pi_s(x, y)) d\gamma(x, y)$ for any function $F$. The second equality comes from the definitions of $\pi_t(x, y)$ and
$\pi_s(x, y)$ and the relation $\|(1-t)x + ty - (1-s)x - sy\| = |t-s||x-y|$. The third equality comes from the fact that $\gamma$ is an optimal transport plan between $\mu$ and $\nu$.

The inequality $W_p(\mu_t, \mu_s) \le |t-s| W_p(\mu, \nu)$ for $t < s$ is actually an equality, because if it were a strict inequality we could write

$$ 
\begin{aligned}
W_p(\mu, \nu) &\le W_p(\mu, \mu_t) + W_p(\mu_t, \mu_s) + W_p(\mu_s, \nu) \\
&< W_p(\mu, \nu)(t + (s-t) + (1-s)) \\
&= W_p(\mu, \nu)
\end{aligned}
$$

and so we would have $W_p(\mu, \nu) < W_p(\mu, \nu)$ which is impossible. The equality $W_p(\mu_t, \mu_s) = |t-s| W_p(\mu, \nu)$ for all $t, s$ shows that $(\mu_t)_t$ is a constant-speed geodesic between $\mu$ and $\nu$ according to Definition 18 and concludes the proof that $W_p(\Omega)$ for $\Omega$ convex is a geodesic space. Convexity of $\Omega$ is needed so that $\pi_t(\Omega \times \Omega) \subset \Omega$. $\square$

If $\Omega$ is convex and $p > 1$, then every constant-speed geodesic in $W_p(\Omega)$ is of the form $((\pi_t)_{\#}\gamma)_t$ for an optimal transport plan $\gamma$.

## B.4 The Dynamic Formulation of Optimal Transport

Chapter 5 of [Santambrogio, 2015] also shows a link between geodesics in $W_p(\Omega)$ (so optimal transport plans) and solutions of the continuity equation

$$ \partial_t \mu_t + \nabla \cdot (\mu_t v_t) = 0, \quad \mu_0 = \mu \qquad (61) $$

where for every $t \in [0, 1]$, $\mu_t$ is a probability measure on $\Omega$ and $v_t : \Omega \to \mathbb{R}^d$ is a vector field. When looking for solutions $(\mu_t, v_t)_t$ to the continuity equation, $\mu_t$ does not have to be a density as we need weak solutions where we only integrate against $\mu_t$.

More precisely, we say that $(\mu_t, v_t)_t$ solves the continuity equation in the weak sense if for any test function $\psi \in C_c^1(\bar{\Omega})$, the function $t \mapsto \int \psi d\mu_t$ is absolutely continuous in $t$ and for almost every $t \in [0, 1]$ we have

$$ \frac{d}{dt} \int_{\Omega} \psi \, d\mu_t = \int_{\Omega} \nabla \psi \cdot v_t \, d\mu_t \qquad (62) $$

If $\mu$ is Lipschitz continuous in $(t, x)$ and $v$ is Lipschitz continuous in $x$ then being a weak solution is equivalent to being a solution almost everywhere.

If is well known that under some conditions $(\Omega$ bounded and $v$ Lipschitz and bounded in $x$ uniformly in $t$ with a flow that remains in $\Omega$), the continuity equation describes the evolution of the density $\mu_t$ of particles distributed initially according to $\mu_0 = \mu$ and moving according to the velocity field $v_t$.
More precisely, given a time-dependent vector field $v$, if points $x \in \Omega$ are distributed according to $\mu_0 = \mu$ and move tangentially to $v$ at all times, i.e satisfy

$$ \partial_t \phi_t(x) = v_t(\phi_t(x)), \phi_0(x) = x \qquad (63) $$

where $\phi_t(x)$ is the position of point $x$ at time $t$, then their distribution across time $((\phi_t)_{\#\mu})_t$ solves the continuity equation (61) for the given vector field. This is the unique solution if we restrict ourselves to distributions that are absolutely continuous with respect to the Lebesgue measure. A solution $\phi$ to (63) is called the flow of the differential equation (63) determined by $v$.

The link with optimal transport comes from the equivalence between absolutely continuous curves in $W_p(\Omega)$ and solutions of the continuity equation.

**Theorem 29.** Take $p > 1$ and $\Omega$ compact.

- If $(\mu_t)_t$ is an absolutely continuous curve in $W_p(\Omega)$, then there exists a vector field $v_t \in L^p(\mathbb{R}^d; \mu_t)$ for every $t \in [0, 1]$ such that $\partial_t \mu_t + \nabla \cdot (\mu_t v_t) = 0$ is satisfied in a weak sense.

- Conversely, if $(\mu_t)_t$ is a family of measure in $P_p(\Omega)$ and we have for every $t \in [0, 1]$ a vector field $v_t \in L^p(\mathbb{R}^d; \mu_t)$ with $\int_0^1 \|v_t\|_{L^p(\mu_t)} dt < \infty$ such that $\partial_t \mu_t + \nabla \cdot (\mu_t v_t) = 0$ is satisfied in a weak sense, then $(\mu_t)_t$ is an absolutely continuous curve in $W_p(\Omega)$.

- Furthermore, $(\mu_t)_t$ and $(v_t)_t$ solving the continuity equation together must satisfy $\|v_t\|_{L^p(\mu_t)} = |\mu'|(t)$ for almost every $t$.

If $\Omega$ is convex and compact and $p > 1$, we can combine this result with the fact that $W_p(\Omega)$ is a geodesic space to obtain another formula for the Wasserstein distance.

**Theorem 30.** On a convex and compact set $\Omega$, for all $\mu, \nu \in P_p(\Omega)$ with $p > 1$ we have

$$ W_p^p(\mu, \nu) = \min \left\{ \int_0^1 \|v_t\|_{L^p(\rho_t)}^p dt \mid \partial_t \rho_t + \nabla \cdot (\rho_t v_t) = 0, \rho_0 = \mu, \rho_1 = \nu \right\} $$

This is known as the *Benamou-Brenier formula*.

*Proof.* First, given Definition 17 of geodesic spaces and formula (58) for the length of absolutely continuous curves we can write for any $\mu, \nu \in P_p(\Omega)$

$$ W_p^p(\mu, \nu) = \left( \min \left\{ \int_0^1 |\rho'|(t) dt \mid \rho \in AC(W_p(\Omega)), \rho_0 = \mu, \rho_1 = \nu \right\} \right)^p $$

Since we proved the existence of constant-speed geodesics between every two points of $W_p(\Omega)$, the minimum can be restricted to these geodesics. Given the
second equivalent definition we gave of constant-speed geodesics when $p > 1$ in Proposition 8, we can write

$$W_p^p(\mu, \nu) = \min \left\{ \int_0^1 |\rho'|^p(t) dt \mid \rho \in \mathrm{AC}(W_p(\Omega)), \rho_0 = \mu, \rho_1 = \nu \right\}$$

Finally, the equivalence between absolutely continuous curves in $W_p(\Omega)$ and solutions of the continuity equation allows us to minimize over these solutions

$$W_p^p(\mu, \nu) = \min \left\{ \int_0^1 \|v_t\|_{L^p(\rho_t)}^p dt \mid \partial_t \rho_t + \nabla \cdot (\rho_t v_t) = 0, \rho_0 = \mu, \rho_1 = \nu \right\}$$

where we also used the relation $\|v_t\|_{L^p(\mu_t)} = |\rho'|(t)$ between $(\rho_t)_t$ and $(v_t)_t$ that jointly solve the continuity equation. □

In this dynamic formulation of the Kantorovich problem we solve for pairs $(\rho_t, v_t)_t$ satisfying the continuity equation with border conditions $\rho_0 = \mu$ and $\rho_1 = \nu$ that minimize the energy $\int_0^1 \|v_t\|_{L^p(\mu_t)}^p$ which is the integral across time of the integral across points of the norms of the tangent vectors to the movement of points distributed initially according to $\mu$, so another way of expressing the total distance traveled by all points.

From an optimal transport plan $\gamma$ between $\mu$ and $\nu$ we can construct the geodesic $((\pi_t)_{\#}\gamma)_t$ and the associated velocity field $v_t$ for $t \in ]0, 1]$ without solving the continuity equation. In fact, if an optimal transport map $T$ exists, then $v_t = (T - \mathrm{id}) \circ T^{-1}$, where $T_t = (1 - t)\mathrm{id} + tT$ and is invertible. Since the continuity equation needs to hold only weakly, mass can still be split between different destinations.

However, we want to solve directly for $(\mu_t, v_t)_t$. As we have mentioned above, if we restrict ourselves to time-dependent vector fields $v$ that are Lipschitz and bounded in $x$ uniformly in $t$ and have a flow $\phi$ that remains in $\Omega$ and to distributions that are absolutely continuous with respect to the Lebesgue measure, then the continuity equation has a unique solution $(\mu_t)$ for a given $v$. The minimization is therefore only over vector fields and we have the problem

$$\begin{array}{ll} \inf_{v} & \displaystyle \int_{0}^{1} \|v_t\|_{L^p(\mu_t)}^{p} dt & (64) \\ \text{subject to} & \partial_t \mu_t + \nabla \cdot (\mu_t v_t) = 0, \mu_0 = \mu, \mu_1 = \nu & \end{array}$$

Since in our machine learning setting distributions $\mu$ and $\nu$ are always discrete samples, we prefer to use the constraint (63) on the movement of points:

$$\begin{array}{ll} \inf_{v} & \displaystyle \int_{0}^{1} \|v_t\|_{L^p((\phi_t)_{\#}\mu)}^{p} dt & (65) \\ \text{subject to} & \partial_t \phi_t(x) = v_t(\phi_t(x)) \quad \forall x \in \Omega, \phi_0 = \mathrm{id}, (\phi_1)_{\#}\mu = \nu & \end{array}$$
where we used the fact that $((\phi_t)_{\#}\mu)_t$ with $\phi$ solution to (63) is a solution to the continuity equation for a given $v$. Again, since in all generality we don't have uniqueness of $\phi$ solution to (63), then the optimization is actually over $(v_t, \phi_t)_t$. Since we are going to discretize the problem and use neural networks to approximate the velocity field, we can limit ourselves to uniformly Lipschitz velocity fields $v$, so that the solution $\phi$ exists and is unique and the problem as expressed in (65) is well defined.

If we want (63) to hold for all $t \in [0, 1]$, then this formulation is not exactly equivalent to the first dynamic formulation (64) and to the Kantorovich problem, as all the mass in a position $x$ will follow the velocity field to a unique point $y = \phi_1(x)$. We are then back to the Monge problem, but here we look for a velocity field that points distributed according to $\mu$ will follow so that they end up distributed according to $\nu$, while minimizing the total energy of this field.

<figure><img src="image_37.png" /><figcaption>Figure 26: Dynamic transport of $\mu_0$ to $\mu_1$ according to velocity field $v$ ($\phi_t^x$ is $\phi_t(x)$ in the text). From [de Bézenac et al., 2021].</figcaption></figure>

## B.5 Regularity of Optimal Transport Maps

Optimal transport maps have some regularity properties under some further assumptions. In particular, Theorem 6.27 of [Ambrosio et al., 2005] gives a classical almost-everywhere regularity result:

**Theorem 31.** If $c(x, y) = \|x-y\|^p$ for $p > 1$, and $\mu$ and $\nu$ have compact supports with $d(\text{supp}(\mu), \text{supp}(\nu)) > 0$, then the optimal transportation map $T$ between $\mu$ and $\nu$ is $\nu$-a.e. differentiable and its Jacobian $\nabla T(x)$ has non-negative eigenvalues $\nu$-almost-surely.

More recently, results summarized below, which correspond to Theorems 4.23, 4.24 and Remark 4.25 of [Figalli, 2017], state that the optimal transportation map has one degree of regularity more than the initial transported density:

**Theorem 32.** Take $c(x, y) = \|x - y\|^p$ for $p > 1$. Suppose there are $X, Y$, bounded open sets, such that the densities of $\mu$ and $\nu$ are null in their respective complements and bounded away from zero and infinity over them respectively.

Then, if $Y$ is convex, there exists $\eta > 0$ such that the optimal transport map $T$
between $\mu$ and $\nu$ is $C^{0,\eta}$ over $X$. If $Y$ is not convex, there exist two relatively
closed sets $A, B$ in $X, Y$ respectively such that $T \in C^{0,\eta}(X \setminus A, Y \setminus B)$, where
$A$ and $B$ are of null Lebesgue measure.

Moreover, if the densities are in $C^{k,\eta}$, then $C^{0,\eta}$ can be replaced by $C^{k+1,\eta}$
in the conclusions above. In particular, if the densities are smooth, then the
transport map is a diffeomorphism.
<figure><img src="image_38.png" /></figure>
# C Background on Gradient Flows

We follow [Santambrogio, 2015, Santambrogio, 2016, Ambrosio et al., 2005, Jordan et al., 1998] for this background on Wasserstein gradient flows.

## C.1 Gradient Flows in Euclidean Space

Given $F: \mathbb{R}^d \to \mathbb{R}$ and $x_0 \in \mathbb{R}^d$, a gradient flow is a curve $x$ solution of the Cauchy problem

$$ \begin{cases} x'(t) = -\nabla F(x(t)) \\ x(0) = x_0 \end{cases} \qquad (66) $$

A solution exists and is unique if $\nabla F$ is Lipschitz or $F$ is convex. Given $\tau > 0$ and $x_0 = x_0$ define a sequence $(x_k^\tau)_k$ through the minimizing movement scheme

$$ x_{k+1}^\tau \in \arg \min_{x \in \mathbb{R}^d} F(x) + \frac{1}{2\tau} \|x - x_k^\tau\|^2 \qquad (67) $$

$F$ lower semi-continuous and $F(x) \ge C_1 - C_2 \|x\|^2$ guarantees existence of a solution of (67) for $\tau$ small enough. $F$ $\lambda$-convex meets these conditions and also provides uniqueness of the solution because of strict convexity of the target.

The optimality condition of problem (67) is

$$ \frac{x_{k+1}^\tau - x_k^\tau}{\tau} = -\nabla F(x_{k+1}^\tau) \qquad (68) $$

which is exactly the formula of the implicit Euler scheme (see Appendix A.1) with step size $\tau$.

We interpret the point $x_k^\tau$ as the value of a curve $x$ at time $k\tau$. We then construct a curve $x^\tau$ as the piecewise constant interpolation of the points $x_k^\tau$

$$ x^\tau(t) = x_k^\tau \text{ for } t \in (k-1)\tau, k\tau], k \in \{1, \dots, K\} \qquad (69) $$

We can also construct a curve $\tilde{x}^\tau$ as the affine interpolation of the points $x_k^\tau$

$$ \tilde{x}^\tau(t) = x_{k-1}^\tau + (t - (k-1)\tau) \frac{x_k^\tau - x_{k-1}^\tau}{\tau} \text{ for } t \in (k-1)\tau, k\tau], k \in \{1, \dots, K\} \quad (70) $$

If $F(x_0) < \infty$ and $\inf F > -\infty$ then $x^\tau$ and $\tilde{x}^\tau$ converge uniformly to the same curve $x$ as $\tau$ goes to zero (up to extracting a subsequence). If $F$ is $C^1$, then the limit curve $x$ is a solution of (66) (i.e. a gradient flow of $F$). If $F$ is not differentiable then $x$ is solution of the problem defined using the subdifferential of $F$, i.e. $x$ satisfies $x'(t) \in -\partial F(x(t))$ for almost every $t$.

If $\lambda > 0$, then the solution to (66) converges exponentially to the unique minimizer of $F$ (which exists by coercivity). So by taking $\tau \to 0$ and $k \to \infty$, we tend towards the minimizer of $F$.
## C.2 Gradient Flows in Metric Spaces

The advantage of the minimizing movement scheme (67) is that it can be adapted to metric spaces. If $(E, d)$ is a metric space, $x_0 \in E$ and $F: E \to \mathbb{R} \cup \{\infty\}$ is lower semi-continuous and bounded from below, we can define a sequence $(x_k^\tau)_k$ with $x_0^\tau = x_0$ and

$$x_{k+1}^\tau \in \arg\min_{x \in E} F(x) + \frac{1}{2\tau} d(x, x_k^\tau)^2 \quad (71)$$

If $E$ is compact then a solution to (71) exists. We study the limit as $\tau$ goes to 0 of the piecewise constant curve $x^\tau$ constructed from these points by (69). In the euclidean setting, we had uniform convergence to a gradient flow of $F$.

A curve $x: [0, T] \to E$ is called a minimizing movement for $F$ and $x_0$ if there exists a sequence of time steps $(\tau_j)_j$ that goes to 0 such that the sequence of piecewise constant curves $(x^{\tau_j})_j$ uniformly converges to $x$ as $j \to \infty$.

Again, if $F(x_0) < \infty$ and $\inf F > -\infty$ then a minimizing movement exists because we can extract a uniformly converging subsequence of curves from the family of piecewise constant curves $(x^\tau)_\tau$ as $\tau \to 0$.

If we add some structure to $(E, d)$ by assuming it is a geodesic space, we can also construct a metric counterpart of the affine interpolation (70) using geodesics

$$\tilde{x}^{\tau}(t) = w_{x_{k-1}^{\tau}, x_k^{\tau}} \left( \frac{t - (k-1)\tau}{\tau} \right) \quad \text{for } t \in [(k-1)\tau, k\tau] \text{ and } k \in \{1, \dots, K\} \quad (72)$$

where $w_{x_{k-1}^{\tau}, x_k^{\tau}}$ is any constant-speed geodesic between $x$ and $y$. The curve $\tilde{x}^\tau$ is continuous and locally Lipschitz (because geodesics are) and coincides with the piecewise constant curve $x^\tau$ at times $k\tau$. If $F(x_0) < \infty$ and $\inf F > -\infty$ then $(\tilde{x}^\tau)_\tau$ is also compact and converges to the same minimizing movement as $x^\tau$ when $\tau \to 0$.

However, if $x$ is a minimizing movement for $F$ we cannot say that it is a gradient flow of $F$, i.e. that it satisfies (66), because $x' = -\nabla F$ does not mean anything in a metric space. We characterize gradient flows in a metric space by finding equivalent conditions to $x' = -\nabla F$ in the euclidean setting that have a metric counterpart. The first such condition is obtained through the Cauchy-Schwartz inequality and is called the energy dissipation equality (EDE)

$$F(x(s)) - F(x(t)) = \int_s^t \frac{1}{2} (\|x'(r)\|^2 + \|\nabla F(x(r))\|^2) dr \quad \forall s < t \quad (73)$$

In the metric setting, we use the metric derivative of a curve $|x'|(t)$ and the descending slope of a real-valued function

$$|\nabla^{-} F|(x) := \limsup_{y \to x} \frac{(F(x) - F(y))_+}{d(x, y)} \quad (74)$$
Under some conditions on the descending slope of F (in particular lower semi-continuity) a minimizing movement for F, i.e. the limit of an interpolation of the minimizing movement scheme, satisfies the EDE condition (73).

If F is λ-convex, another condition which is equivalent to $x'(t) \in -\partial F(x(t))$ in the euclidean setting is the evolution variational inequality (EVI)

$$ \frac{1}{2} \frac{d}{dt} \|x(t) - y\|^2 \leq F(y) - F(x(t)) - \frac{\lambda}{2} \|x(t) - y\|^2 \quad \forall y \in \mathbb{R}^d \quad (75) $$

In the metric setting, we replace the euclidean distance by the distance d and obtain a different definition of gradient flows in metric spaces. The EVI condition leads to uniqueness and stability results for gradient flows and implies the EDE condition, but more assumptions are needed to prove that minimizing movements satisfy it.

### C.3 Gradient Flows in the Wasserstein Space

In the geodesic metric space $W_2(\Omega)$ with $\Omega$ compact and convex, for the functional $F: W_2(\Omega) \to \mathbb{R} \cup \{\infty\}$ lower semi-continuous for the weak convergence of measures in duality with $C(\Omega)$ and $\rho_0^\tau = \rho_0 \in P(\Omega)$, the minimizing movement scheme (71) becomes

$$ \rho_{k+1}^{\tau} \in \arg \min_{\rho \in P_2(\Omega)} F(\rho) + \frac{1}{2\tau} W_2^2(\rho, \rho_k^{\tau}) \quad (76) $$

This problem has a solution because the objective is lower semi-continuous and the minimization is over $P(\Omega)$ which is compact by Banach-Alaoglu.

We can construct a piecewise constant interpolation between the measures $\rho_k^\tau$ as in (69), or a geodesic interpolation as in (72), where we travel along a geodesic between $\rho_{k-1}^\tau$ and $\rho_k^\tau$ in $W_2(\Omega)$ (see Theorem 28 in Appendix B.3 for how to construct geodesics in $W_2(\Omega)$ from optimal transport maps). Again, if $F(x_0) < \infty$ and $\inf F > -\infty$ then both interpolations converge uniformly to a limit curve $(\tilde{\rho}_t)$ as $\tau$ goes to zero. If $F$ is also λ-convex along geodesics, then we might be able to prove the EVI condition for the minimizing movement $\tilde{\rho}$.

However, the space $W_2(\Omega)$ allows us to consider another possible property of the minimizing movement that can also provide uniqueness and stability. We first need to define a notion of differentiation for functionals defined on probability spaces.

**Definition 19.** Given $F: P(\Omega) \to \mathbb{R} \cup \{\infty\}$, we say that $\rho \in P(\Omega)$ is regular for $F$ if $F((1-\varepsilon)\rho + \varepsilon\tilde{\rho}) < \infty$ for all $\varepsilon \in [0, 1]$ and all $\tilde{\rho} \in P(\Omega) \cap L_c^\infty(\Omega)$ (absolutely continuous probability measures with densities in $L_c^\infty(\Omega)$).

**Definition 20.** Given $F: P(\Omega) \to \mathbb{R} \cup \{\infty\}$ and $\rho \in P(\Omega)$ that is regular for $F$ we call the first variation of $F$ at $\rho$ and denote $\frac{\delta F}{\delta \rho}(\rho)$, any measurable function, if it exists, such that

$$ \frac{d}{d\varepsilon} F(\rho + \varepsilon\chi)|_{\varepsilon=0} = \int_{\Omega} \frac{\delta F}{\delta \rho}(\rho) d\chi \quad (77) $$
for every $\chi = \hat{\rho} - \rho$ with $\hat{\rho} \in \mathcal{P}(\Omega) \cap L_c^\infty(\Omega)$.

**Remark 33.** If there is a unique Kantorovitch potential $\varphi$ up to additive constants from $\mu$ to $\nu$ (see Theorem 26 in Appendix B.2), then $\frac{\delta W_p^p(\cdot, \cdot)}{\delta \rho}(\mu) = \varphi$.

We can now say that the minimizing movement $\tilde{\rho}$ may satisfy the following continuity equation:

$$ \partial_t \rho_t(x) - \nabla \cdot (\rho_t(x) \nabla (\frac{\delta F}{\delta \rho}(\rho_t))(x)) = 0 \quad (78) $$

And we can then use the two following results. First from [Santambrogio, 2016]:

**Proposition 9.** If the functional $F: \mathcal{P}(\Omega) \to \mathbb{R} \cup \{\infty\}$ is $\lambda$-geodesically convex and the curves $\rho^0$ and $\rho^1$ satisfy the continuity equation (78) then

$$ \frac{1}{2} \frac{d}{dt} W_2^2(\rho_t^0, \rho_t^1) \leq -\frac{\lambda}{2} W_2^2(\rho_t^0, \rho_t^1) $$

and by Gronwall's inequality

$$ W_2^2(\rho_t^0, \rho_t^1) \leq W_2^2(\rho_0^0, \rho_0^1) e^{-\lambda t} $$

So we have uniqueness of the solution to (78) for the same initial condition $\rho_0$ and stability with exponential convergence as $t \to \infty$ if $\lambda > 0$.

And from [Santambrogio, 2015]:

**Proposition 10.** Given a functional $F: \mathcal{P}(\Omega) \to \mathbb{R} \cup \{\infty\}$ and a probability measure $\rho^*$ that is regular for $F$ and that minimizes it, suppose $\frac{\delta F}{\delta \rho}(\rho^*)$ exists and let $l = \essinf \frac{\delta F}{\delta \rho}(\rho^*)$ be its essential infimum with respect to the Lebesgue measure.

Then if $\frac{\delta F}{\delta \rho}(\rho^*) \in \mathcal{C}(\Omega)$ then $\forall x \in \Omega$, $\frac{\delta F}{\delta \rho}(\rho^*)(x) \geq l$ and $\forall x \in \text{support}(\rho^*)$, $\frac{\delta F}{\delta \rho}(\rho^*)(x) = l$.

And if $\frac{\delta F}{\delta \rho}(\rho^*)$ is only measurable but $\rho^*$ is absolutely continuous, then for almost every $x \in \Omega$, $\frac{\delta F}{\delta \rho}(\rho^*)(x) \geq l$ and for almost every $x \in \{\rho^* > 0\}$, $\frac{\delta F}{\delta \rho}(\rho^*)(x) = l$.

Therefore if $\rho^*$ minimizes $F$ while satisfying the conditions of Proposition 10, then the constant curve $\tilde{\rho}_t = \rho^*$ satisfies the continuity equation (78) in the weak sense or almost everywhere as both terms are equal to zero. If the minimizing movement $\rho$ also satisfies the continuity equation (78) then by Proposition 9, $W_2(\tilde{\rho}_t, \rho^*) \xrightarrow{t \to \infty} 0$ exponentially if $F$ is $\lambda$-geodesically convex for $\lambda > 0$.

Why does it make sense for $\tilde{\rho}$ to satisfy the continuity equation (78)? By Proposition 10 and Remark 33, the solution $\rho_{k+1}^\tau$ to problem (76) satisfies the optimality condition

$$ \frac{\delta F}{\delta \rho}(\rho_{k+1}^\tau) + \frac{\varphi}{\tau} = \text{const} \quad (79) $$
where $\varphi$ is the Kantorovich potential associated with the optimal transport $S$ from $\rho_{k+1}^\tau$ to $\rho_k^\tau$. By Theorem 27 in Appendix B.2 (with $h$ being the Euclidean norm), $S(x) = x - \nabla\varphi(x)$, we consider the velocity $v$ of particles moving from $\rho_k^\tau$ to $\rho_{k+1}^\tau$. According to Appendix B.4, it is given by

$$ -v(x) = \frac{S(x) - x}{\tau} = -\frac{\nabla\varphi(x)}{\tau} = \nabla\left(\frac{\delta F}{\delta\rho}(\rho_{k+1}^\tau)\right)(x) \quad (80) $$

Therefore, as $\tau \to 0$, we can expect to get in $\tilde{\rho}$ a solution to the continuity equation (78), but this has to be proven rigorously for a given $F$.
<figure><img src="image_39.png" /></figure>
# D Background on Module-wise Training

The typical setting of (sequential) module-wise training for minimizing a loss $L$, is, given a dataset $\mathcal{D}$, to solve one after the other, for $1 \le k \le K$, Problems

$$ (T_k, F_k) \in \arg \min_{T,F} \sum_{x \in \mathcal{D}} L(F, T(G_{k-1}(x))) $$

where $G_k = T_k \circ \dots \circ T_1$ for $1 \le k \le K$ and $G_0$=id. Here, $T_k$ is the module (one or many layers) and it receives the output of module $T_{k-1}$, and $F_k$ is an auxiliary classifier that processes the outputs of $T_k$ so the loss can be computed. The final network trained this way is $F_K \circ G_K$, but we can stop at any depth $k$ and use $F_k \circ G_k$ if it performs better.

In fact, module-wise training suffers often from a well-documented stagnation problem observed in [Marquez et al., 2018, Belilovsky et al., 2019, Wang et al., 2021, Pyeon et al., 2021], whereby greedy early modules overfit and learn more discriminative features than end-to-end training, and deeper modules don't improve the test accuracy significantly, or even degrade it.

We present below some recent methods that aim at alleviating this problem and improving the accuracy of module-wise training, either by adding other terms to the loss used to train the modules, or simply through architectural considerations.

## D.1 Additional Loss Terms

**InfoPro [Wang et al., 2021].** InfoPro start by verifying that vanilla module-wise training outperforms end-to-end training in the first modules, but then destroys too much information between the learned features and both the inputs and the labels and stagnates and gets overtaken by end-to-end training in the later modules. See Figure 27 below.

<figure><img src="image_40.png" /><figcaption>Figure 27: Linear separability (left), mutual information with input (center), and mutual information with output (right) of features learned with module-wise training with $K$ modules. $K=1$ is end-to-end training. From [Wang et al., 2021].</figcaption></figure>
They then propose to train each module to maximize the mutual information between the learned features h and both the input x and its label y. This means minimizing

$$-\lambda_1 I(h, x) - \lambda_2 I(h, y)$$

where I is the mutual information given by $I(a, b) = H(b) - H(b|a)$ where H is the entropy [Shannon, 1948]. The mutual information $I(h, y)$ with the label is approximated using an auxiliary classifier with the cross-entropy loss or a contrastive loss inspired by [Chen et al., 2020b, Khosla et al., 2020, He et al., 2020]. The mutual information $I(h, x)$ with the input is approximated using an auxiliary decoder that tries to reconstruct x from h with a reconstruction loss that is binary cross-entropy.

In [Pathak et al., 2022], this method is slightly modified in that the auxiliary decoder now tries to reconstruct the previous feature (i.e the input to the module) and not the initial input x to the network.

**PredSim [Nøkland and Eidnes, 2019].** PredSim adds an auxiliary similarity matching loss inspired by similarity losses used in neuroscience [Patel et al., 2022]. The similarity matching loss is the distance between matrices whose entries are meant to measure the similarity between inputs in a mini-batch. The target for such a matrix is then the similarity matrix constructed from the labels of the inputs, which contains 1 at position (i, j) if inputs i and j have the same label, and 0 otherwise. More precisely, given a matrix A made up of stacked vectors $a_i$, we denote $S(A)$ the adjusted cosine similarity matrix or correlation matrix, whose entries are given by

$$S(A)_{i,j} = S(A)_{j,i} = \frac{\tilde{a}_i^\top \tilde{a}_j}{\|\tilde{a}_i\|_2 \|\tilde{a}_j\|_2}$$

where we mean-center vector $a_i$ into vector $\tilde{a}_i$. Given a mini-batch, we denote Y the one-hot encoded matrix of its labels, and H the matrix of the stacked embeddings $h_i$ of the elements of the mini-batch after the current module. The similarity matching loss is then $\|S(\text{Conv}(H) - S(Y))\|_F^2$, meaning that the embeddings go through a convolutional auxiliary layer to reduce their dimension. This loss is added to the cross-entropy classification loss for the training of every module.

## D.2 Architectural Approaches

**DGL [Belilovsky et al., 2020].** DGL, and previously [Belilovsky et al., 2019], simply consider the best auxiliary architectures for module-wise training. They find that deeper convolutional auxiliary classifiers lead to the best performances. However, these auxiliary classifiers are quite large which reduces the computational savings we desire from module-wise training. We, like InfoPro [Wang et al., 2021], limit ourselves to one convolutional layer. DGL then introduce an auxiliary architecture made up of staged resolution reduction phase before the
linear classification layer, which they find offers a good trade-off between size
and performance.

**Sedona [Pyeon et al., 2021].** Sedona designs light auxiliary classifiers inspired by the architecture of MobileNetV2 [Sandler et al., 2018] that relies on depth-wise and point-wise convolutions and inverted residual blocks. They also introduce an architecture search phase that chooses the best auxiliary classifier architecture from a pool of candidates, and where exactly to insert them in the network. This search phase however leads to a deep early module which reduces the memory savings of module-wise training.
<figure><img src="image_41.png" /></figure>
# E Background on Adversarial Attacks

## E.1 Adversarial Attacks

We present here how some popular adversarial attacks on neural networks that solve an image classification task work. We refer to [Aldahdooh et al., 2022, Li et al., 2022] for more details.

Below, $L$ is the classification loss (usually cross-entropy), $x$ is a clean image, $y$ is its label and $d$ its dimension as a vector, $\tilde{x}$ is the adversarial image found from $x$, $\delta = x - \tilde{x}$ is the perturbation, and $\epsilon > 0$ is the maximal allowed perturbation size in a given $L_p$ norm. For targeted attacks, $t \neq y$ is the target label the attack wants the network to predict for $\tilde{x}$. We denote the pre-softmax outputs of the network $Z$, while $f$ denotes the softmax outputs.

### E.1.1 White-box Attacks

White-box attacks are attacks where the adversary has access to the network's architecture, weights and gradients. We present below three gradient-based attacks (FGM, BIM and PGD), one optimization-based attack (CW), which is one of the strongest and most widely used attacks, and one boundary-based attack (DF).

**FGM [Goodfellow et al., 2015b].** The Fast Gradient Method is a simple early attack. It finds the update direction that maximizes the loss $L$ at each pixel of image $x$ and takes one step in that direction. So

$$ \tilde{x} = x + \epsilon \operatorname{sign}(\nabla_x L(x, y)) \text{ such that } \tilde{x} \in [0, 1]^d $$

This way $\|\tilde{x} - x\|_{\infty} \le \epsilon$.

**BIM [Kurakin et al., 2017].** The Basic Iterative Method iterates the FGM attack $k$ times. So for $0 \le i < k$ and $\tilde{x}_0 = x$,

$$ \tilde{x}_{i+1} = \tilde{x}_i + \alpha \operatorname{sign}(\nabla_x L(\tilde{x}_i, y)) \text{ such that } \tilde{x}_{i+1} \in [0, 1]^d $$

where $\alpha = \epsilon/k$ and $\tilde{x} = \tilde{x}_k$

**PGD [Madry et al., 2018].** Projected Gradient Descent is similar to BIM but starts from a random perturbation in an $L_p$ ball around $x$. Auto-PGD [Croce and Hein, 2020b] adds a momentum to the gradient step and adapts the step size $\alpha$ across iterations.

**CW [Carlini and Wagner, 2017b].** The Carlini-Wagner attack is one of the strongest attacks available. It considers a loss

$$ g(x) = \max(\max_{j \neq t}(Z_j(x)) - Z_t(x), -c) $$

where $c > 0$ is a confidence parameter. Minimizing $g$ finds an image with a high softmax output the target $t$. So the attack then solves the optimization problem

$$ \min_{\delta} \|\delta\|_{p} + c g(x + \delta) \text{ such that } \tilde{x} = x + \delta \in [0, 1]^d $$
The attack also transforms the problem from a box-constrained one to an unconstrained one by rewriting the perturbation as $\delta = \frac{1}{2}(\tanh(w) + 1) - x$.

**DF [Moosavi-Dezfooli et al., 2016].** The DeepFool attack approximates the decision regions of a neural network by polyhedrons, taking inspiration from the case of an affine classifier. The attack then accumulates perturbations that are the shortest projections to the boundary of the polyhedron of the correct label. First the attack computes $w_j = \nabla f_j(x) - \nabla f_y(x)$ and $d_j = f_j(x) - f_y(x)$ for all labels $j$ that are different from the ground truth label $y$. It then finds the closest boundary thus

$$k = \arg \min_{j \neq y} \frac{|d_j|}{\|w_j\|_2}$$

The perturbation added to $x$ is then

$$\delta = \frac{|d_k|}{\|w_k\|_2^2} w_k$$

This process is repeated starting from $x + \delta$ until the predicted label changes.

Among other white-box attack we mention JSMA (Jacobian Saliency Map Attack) [Papernot et al., 2016b], UAP (Universal Adversarial Perturbation) [Moosavi-Dezfooli et al., 2017b], FA (Feature Adversary) [Sabour et al., 2016], ATN (Adversarial Transformation Networks) [Baluja and Fischer, 2017] and FAB (Fast Adaptive Boundary) [Croce and Hein, 2020a].

## E.1.2 Black-box Attacks

Black-box attacks don't have any knowledge about the network and can only query it. We present below the substitute model family of black-box attacks, a random search attack and an optimization-based attack.

**Substitute model attacks [Papernot et al., 2016a, Papernot et al., 2017, Hu and Tan, 2022].** These attacks simply train a network (called a substitute model) that imitates the predictions of the targeted network. The attacker then uses a white-box attack to find an adversarial sample for the substitute model and feeds this sample to the targeted network.

**SA [Andriushchenko et al., 2020].** The Square Attack selects via random search local square-shaped $\epsilon$-bounded perturbations at random positions to solve the optimization problem

$$\min_{\tilde{x}} \left( Z_y(\tilde{x}) - \max_{j \neq y} Z_j(\tilde{x}) \cdot Z_y(\tilde{x}) \right) \text{ such that } \tilde{x} \in [0, 1]^d \text{ and } \|\tilde{x} - x\|_p \le \epsilon$$

The attack has an $L_2$ variant and an $L_\infty$ variant.

**ZOO [Chen et al., 2017b].** The Zeroth Order Optimization attack estimates the gradients of the network $f$ using finite differences. Indeed

$$\frac{\mathrm{d}f(x)}{\mathrm{d}x_i} \approx \frac{f(x + he_i) - f(x - he_i)}{2h}$$
and

$$ \frac{\mathrm{d}^2 f(x)}{\mathrm{d}^2 x_i} \approx \frac{f(x + he_i) - 2f(x) + f(x - he_i)}{h^2} $$

where $h > 0$ is small and $e_i$ is the $i$-th standard basis vector. In one variant of the attack, these estimations are used to perform stochastic coordinate descent optimization methods [Bertsekas, 2016] to minimize a similar target to that of the Carlini-Wagner attack:

$$ g(x) = \max(\max_{j \neq t}(\log f_i(x)) - \log f_t(x), -c) $$

Among other black-box attacks we mention Pixel Attacks [Su et al., 2019, Kotyan and Vargas, 2020], ST (spatial transformation) [Engstrom et al., 2019], HSJ (Hop-Skip-Jump) [Chen et al., 2020a] and BA (Boundary Attack) [Brendel et al., 2018].

## E.2 Adversarial Detectors

We describe here how some recent adversarial detectors operate. Generally, an adversarial detector extracts statistics from the activations of an input going through a network, before training a classifier on these statistics.

### E.2.1 Statistical Approaches

Mahalanobis [Lee et al., 2018]. The Mahalanobis detector uses Mahalanobis distances between the activations and a Gaussian fitted to them during network training as detection features. The means of the Gaussian distributions are fitted for each class $c$ thus

$$ \hat{\mu}_{c,m} = \frac{1}{N_c} \sum_{x \in D_c} h_m(x) $$

where $D_c$ is the set of samples $x$ with class $c$, $N_c = \#D_c$ is its cardinality, and $h_m$ is the embedding function at depth $m$. The covariance matrix is the same for all classes and is fitted thus

$$ \hat{\Sigma}_m = \frac{1}{N} \sum_{(x,y) \in D} (h_m(x) - \hat{\mu}_{y,m})(h_m(x) - \hat{\mu}_{y,m})^T $$

where $D$ is the dataset of sample-label pairs, and $N = \#D$ is the number of samples. The Mahalanobis distances to these Gaussians for an input $x$ are then for every class $c$

$$ \hat{M}_{c,m}(x) = -(h_m(x) - \hat{\mu}_{c,m})^T \hat{\Sigma}^{-1} (h_m(x) - \hat{\mu}_{c,m}) $$

The confidence scores for an input $x$ to be classified as adversarial (or out-of-distribution) or not are then

$$ \hat{S}_m(x) = \max_c - \hat{M}_{c,m}(x) $$
In practice, the detector (a logistic regression) is trained on Mahalanobis scores $\hat{S}_m(x)$ at five different depths $m$.

[Ren et al., 2021] propose an improvement over the Mahalanobis distance for out-of-distribution detection called RMD (Relative Mahalanobis Distance). They fit another Gaussian to all training samples, regardless of class. So the estimated mean is now

$$ \tilde{\mu}_m = \frac{1}{N} \sum_{(x,y) \in \mathcal{D}} h_m(x) $$

And the estimated covariance is now

$$ \tilde{\Sigma}_m = \frac{1}{N} \sum_{(x,y) \in \mathcal{D}} (h_m(x) - \tilde{\mu}_m)(h_m(x) - \tilde{\mu}_m)^{\top} $$

Then they compute the Mahalanobis distance to this Gaussian

$$ \tilde{M}_m(x) = -(h_m(x) - \tilde{\mu}_m)^{\top} \tilde{\Sigma}^{-1} (h_m(x) - \tilde{\mu}_m) $$

And define the relative Mahalanobis distance as follows

$$ \mathrm{RMD}_{c,m}(x) = \hat{M}_{c,m}(x) - \tilde{M}_m(x) $$

The confidence scores used as features for detection are now

$$ \tilde{S}_m(x) = \max_c -\mathrm{RMD}_{c,m}(x) $$

**LID [Ma et al., 2018].** Given a point $x$ sampled from a probability distribution $P$, and a batch $B = \{x_i\}_{i=1}^k$ of points also sampled i.i.d from $P$, the MLE estimator of the *local intrinsic dimensionality* (LID) of $P$ at point $x$ is given by

$$ \hat{L}(x, B) = - \left( \frac{1}{k} \sum_{i=1}^{k} \log \frac{r_i(x)}{r_k(x)} \right)^{-1} $$

where $r_i(x)$ is the distance from $x$ to its $i$-th nearest neighbor in $B$. For adversarial detection, the features used are $\hat{L}(a, B)$ where $a$ is an activation at a certain depth from an input (clean or adversarial) going through the network, and $B$ is always a batch of activations of clean samples.

## E.2.2 Computer Vision Approaches

**NSS [Kherchouche et al., 2020].** The natural scene statistics (NSS) detector uses image quality assessment statistics that detect distortions to detect adversarial samples. It extracts 18 of these statistics from the images themselves and not from their embeddings. We describe here how the first two statistics are extracted.
The mean subtracted contrast normalized (MSCN) coefficients [Mittal et al., 2012] of an image I are given by

$$ \hat{I}(i, j) = \frac{I(i, j) - \mu(i, j)}{\phi(i, j) + c} $$

where $c > 0$ is a small constant, $(i, j)$ is a pixel position, the local mean $\mu$ is

$$ \mu(i, j) = \sum_{k=-3}^{3} \sum_{l=-3}^{3} w(k, l)I(i + k, j + l) $$

and the local standard deviation $\phi$ is

$$ \phi(i, j) = \sqrt{\sum_{k=-3}^{3} \sum_{l=-3}^{3} w(k, l)(I(i + k, j + l) - \mu(i, j))^2} $$

where $w$ is a circularly-symmetric Gaussian weighting function. A generalized Gaussian distribution (GGD) is then fitted to $\hat{I}$ and its shape parameter and variance are used as the first two features for detection.

**Fourier detectors [Harder et al., 2021].** The discrete Fourier transform (DFT) in 2 dimensions of an image $I \in [0, 1]^{d \times d}$ is given by

$$ \mathcal{F}(I)(l, m) = \sum_{j,k=0}^{d} I(j, k) \exp -2\pi i \frac{l j + k m}{d} $$

for $l, k \in \{0, \dots, d-1\}$. The magnitude of the Fourier coefficients is given by

$$ |\mathcal{F}(I)(l, m)| = \sqrt{\mathrm{Re}^2(\mathcal{F}(I)(l, m)) + \mathrm{Im}^2(\mathcal{F}(I)(l, m))} $$

And the phase is given by

$$ \phi(\mathcal{F}(I)(l, m)) = \arctan \frac{\mathrm{Re}(\mathcal{F}(I)(l, m))}{\mathrm{Im}(\mathcal{F}(I)(l, m))} $$

[Harder et al., 2021] propose two detectors. This first is called MFS (magnitude Fourier spectrum) and uses the magnitudes of the Fourier coefficients as detection features. The second is called PFS (phase Fourier spectrum) and uses the phases of the Fourier coefficients as detection features. They find it is better to extract these features from the embeddings of the images at certain layers than from the images themselves. But given the high dimensionality of these embeddings, they have to heuristically select a subset of layers to consider.
<figure><img src="image_42.png" /></figure>
# F Additional Experiments from Section 3

## F.1 Implementation Details

These are further implementation details for the experiments in Section 3.5. Orthogonal initialization with gain 0.01 is used for ResNet models. Kaiming initialization is used for ResNeXt models. SGD is used for training. The momentum is always set to 0.9 and weight decay to 0.0001. For ResNet models, the learning rate is 0.01 and is divided by 5 at epochs 120, 160 and 200 (when the training goes that far). For ResNeXt models, the learning rate is 0.1 and is divided by 10 at epochs 150, 225 and 250. Batch size is 128 for all experiments. Architectures of ResNet [He et al., 2016a] and ResNeXt [Xie et al., 2017] blocks are standard and exactly as in the references. The ResNets used are single representation ResNets (i.e. containing one residual stage) with 9 blocks. The ResNeXt architecture used is the ResNeXt-50-32×4d from [Xie et al., 2017].

## F.2 Transport Visualization

If we pretrain an autoencoder on MNIST, we can use its encoder as the encoder of the ResNet and freeze it during training. This makes it possible to visualize the transport of the data by decoding, via the pretrained decoder, the output of each residual block (see Figure 28 below).

<figure><img src="image_43.png" /><figcaption>Figure 28: Decoding the internal representations of a network.</figcaption></figure>
We show this below on MNIST. In Figure 29, we see the decodings of a vanilla
ResNet trained to achieve 99.4% test accuracy.

<figure><img src="image_44.png" /><figcaption>Figure 29: Decodings of the internal representations (outputs of the ResBlocks) after training a ResNet-9 on MNIST (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>

We add the transport cost with $\lambda_0=5$, $\tau=1$ and $s=5$. The accuracy barely drops (99.3%) and we see in Figure 30 that the decodings change much less.

<figure><img src="image_45.png" /><figcaption>Figure 30: Decodings of the internal representations (outputs of the ResBlocks) after training a LAP-ResNet-9 on MNIST (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>
We do the same thing on CIFAR10, first with a vanilla ResNet-9 in Figure 31.

<figure><img src="image_46.png" /><figcaption>Figure 31: Decodings of the internal representations (outputs of the ResBlocks) after training a ResNet-9 on CIFAR10 (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>

Then with LAP training with $\lambda_0=10$, $\tau=0.1$ and $s=50$ in Figure 32 below.

<figure><img src="image_47.png" /><figcaption>Figure 32: Decodings of the internal representations (outputs of the ResBlocks) after training a LAP-ResNet-9 on CIFAR10 (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>
Finally, on CIFAR100, the decodings without regularization are in Figure 33.

<figure><img src="image_48.png" /><figcaption>Figure 33: Decodings of the internal representations (outputs of the ResBlocks) after training a ResNet-9 on CIFAR100 (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>

With the transport regularization, the decodings are in Figure 34.

<figure><img src="image_49.png" /><figcaption>Figure 34: Decodings of the internal representations (outputs of the ResBlocks) after training a LAP-ResNet-9 on CIFAR100 (og: original image, ae: encoding, b1: output of block 1...)</figcaption></figure>
## F.3 Experiments with Fixed $\lambda$

In this section, we show some additional experimental results with a model where, instead of using an adaptive optimization algorithm, we simply take the transport cost as a regularizer, thus giving us a minimization objective:

$$ \mathcal{L}(\theta) + \lambda \mathcal{C}(\theta) $$

This is an easier and more straightforward method, simply considering a relaxed constraint in the optimization problem (18). Aside from the advantage of simpler implementation, it allows for easier fine-tuning of the regularization hyper-parameter which is useful when the datasets and networks are big. The adaptivity is lost but this still leads to better test accuracy than the non-regularized networks. Results on the same tasks as in Section 3.5 are below.

<table><caption>Table 16: Average highest test accuracy and 95% confidence interval of ResNeXt50 over 10 instances on CIFAR100 with training sets of different sizes (in %)</caption><thead><tr><th>Train size</th><th>ResNeXt</th><th>LAP-ResNeXt50</th><th>ResNeXt50 with λ = 0.01</th></tr></thead><tbody><tr><td>50 000</td><td>72.97 ± 1.18</td><td><b>76.11</b> ± 0.78</td><td>75.96 ± 1.04</td></tr><tr><td>25 000</td><td>62.55 ± 2.37</td><td><b>64.11</b> ± 1.85</td><td>64.10 ± 1.74</td></tr><tr><td>12 500</td><td>45.90 ± 2.74</td><td><b>48.23</b> ± 1.84</td><td>47.77 ± 1.85</td></tr></tbody></table>

<table><caption>Table 17: Average highest test accuracy and 95% confidence interval of ResNet-9 over 10 instances on CIFAR100 with training sets of different sizes (in %)</caption><thead><tr><th>Train size</th><th>ResNet</th><th>LAP-ResNet</th><th>ResNet with λ ∈ {0.05, 0.2}</th></tr></thead><tbody><tr><td>50 000</td><td>72.32 ± 0.24</td><td>72.43 ± 0.18</td><td><b>72.62</b> ± 0.21</td></tr><tr><td>25 000</td><td>64.34 ± 0.23</td><td>64.34 ± 0.23</td><td><b>64.76</b> ± 0.24</td></tr><tr><td>10 000</td><td>49.27 ± 0.42</td><td><b>50.57</b> ± 0.23</td><td>50.46 ± 0.27</td></tr><tr><td>5 000</td><td>34.74 ± 0.84</td><td>37.97 ± 0.30</td><td><b>38.44</b> ± 0.45</td></tr><tr><td>1 000</td><td>15.66 ± 0.43</td><td><b>16.42</b> ± 0.32</td><td>16.03 ± 0.49</td></tr></tbody></table>

<table><caption>Table 18: Average highest test accuracy and 95% confidence interval of ResNet-9 over 20 instances on CIFAR10 with training sets of different sizes (in %)</caption><thead><tr><th>Train size</th><th>ResNet</th><th>LAP-ResNet</th><th>ResNet with λ = 0.2</th></tr></thead><tbody><tr><td>50 000</td><td>91.49 ± 0.09</td><td><b>91.94</b> ± 0.10</td><td>91.36 ± 0.08</td></tr><tr><td>30 000</td><td>88.61 ± 0.14</td><td><b>89.41</b> ± 0.09</td><td>88.50 ± 0.11</td></tr><tr><td>20 000</td><td>85.73 ± 0.14</td><td><b>86.74</b> ± 0.13</td><td>85.82 ± 0.11</td></tr><tr><td>10 000</td><td>79.25 ± 0.25</td><td><b>80.90</b> ± 0.16</td><td>80.15 ± 0.13</td></tr><tr><td>5 000</td><td>70.32 ± 0.32</td><td><b>72.58</b> ± 0.12</td><td>72.03 ± 0.31</td></tr><tr><td>4 000</td><td>67.80 ± 0.26</td><td><b>70.12</b> ± 0.30</td><td>69.64 ± 0.30</td></tr><tr><td>1 000</td><td>49.22 ± 0.52</td><td><b>51.14</b> ± 0.45</td><td>50.38 ± 0.45</td></tr><tr><td>500</td><td>41.55 ± 0.41</td><td><b>42.92</b> ± 0.37</td><td>42.30 ± 0.43</td></tr><tr><td>100</td><td>26.98 ± 0.99</td><td>25.34 ± 0.76</td><td><b>27.53</b> ± 0.94</td></tr></tbody></table>
## F.4 Impact on Weight Distribution of Classifier

LAP training moves the weights of the classification layers away from zero: compare Figure 35 for a vanilla network and Figure 36 for a LAP network.

<figure><img src="image_50.png" /><figcaption>histogram of weights of the second fully connected layer of the classifier</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 35: Weight distribution of the two fully connected layers of the classification head of a vanilla ResNet-9 trained on CIFAR10.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>histogram of weights of the second fully connected layer of the classifier</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>Figure 36: Weight distribution of the two fully connected layers of the classification head of a LAP-ResNet-9 trained on CIFAR10.</figcaption></figure>
# F.5 Additional Experiments in 2 Dimensions

## F.5.1 Comparison to Exact Wasserstein Distances

In 2 dimensions we can compute the exact $W_2$ distance between successive point clouds in the network before and after each block, and compare it, in Figure 37 below, to the displacement that the block causes. This confirms that residual networks are indeed already close to being optimal transport maps, and that LAP training encourages this bias.

<figure><img src="image_54.png" /><figcaption>Figure 37: Average residual function norm and $W_2$ distance between successive point clouds in a ResNet-9 trained on the circles test set. Vanilla training on the left and LAP training on the right.</figcaption></figure>

## F.5.2 Comparison to Batch Normalization

We compare our method with batch normalization. We find that the two cooperate well to improve test accuracy on the same 2-dimensional task as in Section 3.2 when the model is too small (1 block, Table 19), too big (100 blocks, Table 20), poorly initialized ($\mathcal{N}(0, 5)$ initialization, Table 22) and when the data set is small (50 points, Table 21). LAP-ResNets use $\lambda_0 = 0.1$, $\tau = 0.1$ and $s = 5$.

We are interested in cases where adding the transport cost to the loss helps with the generalization of the network. The advantage of Uzawa's algorithm is that fixing $\lambda_0 = 0.1$, $\tau = 0.1$ and $s = 5$ works well in most cases. A fixed $\lambda$ requires more tinkering, but when a good value is found, it often performs
better than Uzawa's algorithm with those values, and finding better values for Uzawa requires a grid search over three parameters instead of just one. Since the results are similar for both the circles and the moons dataset, we report below the results on the circles dataset only.

A first case where adding $C$ to the loss helps is when the network is small. If the network has only two blocks, then it takes on average around 73 epochs to achieve 100% test accuracy and keep that performance for three consecutive epochs (on average over 100 instances of the problem, an instance being a randomly generated circles dataset and a random initialization of the model's weights). If we add batch normalization to the residual blocks as in Figure 1, then the average number of epochs is around 71. But if we add the transport cost $C$ with fixed weight $\lambda = 0.01$, then the average number of epochs to achieve the same performance drops to 54. Adding both batch normalization and the transport cost at the same time does not help, as the average number of epochs goes up to 64.

If the network has only one block (Table 19), then it never reaches 100% test accuracy. The highest test accuracy it reaches after 100 epochs is on average 76.6%. With batch normalization, this figure drops to 75.4%, but with Uzawa's algorithm, it goes up to 82.1%. With both batch normalization and Uzawa, the average best accuracy is even higher at 84.6%.

<table><caption>Table 19: Average test accuracy and 95% confidence interval over 100 instances on the circles 2D dataset with 1000 points and 1 block (in %)</caption><thead><tr><th></th><th>No batch normalization</th><th>Batch normalization</th></tr></thead><tbody><tr><td>ResNet</td><td>76.6 ± 3.5</td><td>75.4 ± 3.1</td></tr><tr><td>ResNet with λ = 0.005</td><td>76.5 ± 3.5</td><td>75.6, ± 3.3</td></tr><tr><td>LAP-ResNet</td><td>82.1 ± 2.6</td><td><b>84.6</b> ± 3.1</td></tr></tbody></table>

A second case is when the network is too big (Table 20). The network with 100 blocks and no batch normalization has an average test performance of only 89.1%. As expected, adding batch normalization helps greatly as the average test accuracy rises to 99.4%. This performance is improved slightly by adding a transport cost via Uzawa's algorithm and reaches 99.8%. Note that adding the transport cost without batch normalization hinders the test performance in this case.

<table><caption>Table 20: Average test accuracy and 95% confidence interval over 100 instances on the circles 2D dataset with 1000 points and 100 blocks (in %)</caption><thead><tr><th></th><th>No batch normalization</th><th>Batch normalization</th></tr></thead><tbody><tr><td>ResNet</td><td>89.1 ± 0.9</td><td>99.4 ± 0.4</td></tr><tr><td>ResNet with λ = 0.09</td><td>69.7 ± 4.0</td><td>99.5 ± 0.6</td></tr><tr><td>LAP-ResNet</td><td>75.7 ± 2.9</td><td><b>99.8</b> ± 0.2</td></tr></tbody></table>
A third instance of the transport cost improving generalization is when we use a big initialization (Table 22). The network has 9 blocks and we initialize the weights and biases of the fully connected layers by drawing from a centered normal distribution with standard deviation equal to 5. Here, batch normalization is essential for the network to recover from this initialization and get close to 100% test accuracy inside 100 epochs. If we further add a transport cost with $\lambda = 0.04$, the average highest test accuracy improves from 97.8% to 99.8%. If we add a transport cost via Uzawa's algorithm instead of a fixed weight, we get an average test accuracy of 99.6%. An example of the behavior we want to avoid with a big initialization is in Figure 39 below. With batch normalization (and the same random seeds), we get the behavior in Figure 40. If we also add $\lambda = 0.04$, we gain further control as seen in Figure 41.

<table><caption>Table 21: Average test accuracy and 95% confidence interval over 100 instances on the circles 2D dataset with 50 points and 9 blocks (in %)</caption><thead><tr><th></th><th>No batch normalization</th><th>Batch normalization</th></tr></thead><tbody><tr><td>ResNet</td><td>88.2 ± 2.3</td><td>92.9 ± 2.0</td></tr><tr><td>ResNet with λ = 0.04</td><td>93.5 ± 2.1</td><td>94.4 ± 2</td></tr><tr><td>LAP-ResNet</td><td>95.8 ± 1.8</td><td><b>96.0 ± 1.3</b></td></tr></tbody></table>

A final case is when the training set is small (Table 21). Here we consider a dataset containing 50 points (40 train and 10 test) and a network with 9 blocks. The test performance of this network is on average 88.2%. We can get to 92.7% average test accuracy with batch normalization but using Uzawa's algorithm and no batch normalization yields the best average performance at 95.7%. Adding batch normalization does not change this.

<table><caption>Table 22: Average test accuracy and 95% confidence interval over 100 instances on the circles 2D dataset with a $\mathcal{N}(0, 5)$ initialization (in %)</caption><thead><tr><th></th><th>No batch normalization</th><th>Batch normalization</th></tr></thead><tbody><tr><td>ResNet</td><td>90.2 ± 1.3</td><td>98.0 ± 0.8</td></tr><tr><td>ResNet with λ = 0.04</td><td>89.7 ± 1.5</td><td><b>99.7 ± 0.2</b></td></tr><tr><td>LAP-ResNet</td><td>79.1 ± 3.9</td><td>99.4 ± 0.4</td></tr></tbody></table>

### F.5.3 Additional Transport Visualizations in 2 Dimensions

We provide here further visualizations of transport dynamics in 2 dimensions inside a ResNet-9, as in Section 3.2.2, but now showing point clouds after every block. In Figure 38, we show the normal behavior of a network trained with transport regularization (to be compared with the first row in Figure 12 in Section 3.2.2). The three following figures show how transport regularization and batch normalization help control the behavior of a badly initialized network.
<figure><img src="image_55.png" /><figcaption>Figure 38: Circles test set after each block after training with transport regularization with weight $\lambda = 0.07$.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>Figure 39: Circles test set after each block after training with $N(0, 5)$ initialization.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>Figure 40: Circles test set after each block after training with $N(0, 5)$ initialization and batch normalization.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>Figure 41: Circles test set after each block after training with $\mathcal{N}(0, 5)$ initialization, batch normalization and transport regularization with fixed weight $\lambda = 0.04$.</figcaption></figure>
<figure><img src="image_59.png" /></figure>
# G Additional Experiments from Section 4

## G.1 Implementation Details

We use standard data augmentation and standard implementations for VGG-19, ResNet-50, ResNet-101, ResNet-110, ResNet-152 and Swin-Tiny Transformer (the same as for the other methods). We use NVIDIA Tesla V100 16GB GPUs for the experiments.

Training a Resnet-152 on TinyImageNet in Table 4 takes about 36 hours. Training a Resnet-152 on CIFAR100 in Table 4 takes about 11 hours. Training a ResNet-110 on STL10 in Table 6 takes about 3 hours. Training a Swin-Tiny Transformer in Table 7 take between 40 minutes and 1 hour.

For sequential and multi-lap sequential training, we use SGD with a learning rate of 0.007. With the exception of the Swin Transformer in Table 7, we use SGD for parallel training with learning rate of 0.003 in all Tables but Table 6 where the learning rate is 0.002. For the Swin Transformer in Table 7, we use the AdamW optimizer with a learning rate of 0.007 and a CosineLR scheduler.

For end-to-end training we use a learning rate of 0.1 that is divided by five at epochs 120, 160 and 200. Momentum is always 0.9. For parallel and end-to-end training, we train for 300 epochs. For sequential and multi-lap sequential training, the number of epochs varies per module (see Section 4.5.5).

For experiments in Section 4.5.1, we use a batch size of 256, orthogonal initialization [Saxe et al., 2014] with a gain of 0.1, label smoothing of 0.1 and weight decay of 0.0002. Only for Table 6, the batch size changes to 64.

For experiments in Section 4.5.5, we use a batch size of 128, orthogonal initialization with a gain of 0.05, no label smoothing and weight decay of 0.0001.

For the experiments on TinyImageNet in Table 4, we use $\tau = 500000$ for the first two modules and then double it for the last two modules for TRGL.

For the experiments on CIFAR100 in Table 5, we use $\lambda_{k,1} = 1, h = 1$ and $s = 50$ for TRGL.

For the experiments on STL10 in Table 6, we use $\tau = 0.5$ and double it at the midpoint, expect for the first line where $\tau = 50$.

## G.2 Additional Experiments

<table><caption>Table 23: Test accuracy of 2-7 ResNet on CIFAR100 with parallel TRGL and VanGL, compared to DDG and FR from [Huo et al., 2018a] that also split their networks in 2 module-wise-parallel-trained modules.</caption><thead><tr><th>Parallel VanGL (ours)</th><th>Parallel TRGL (ours)</th><th>DDG ResNet-101</th><th>FR ResNet-101</th></tr></thead><tbody><tr><td>76.26</td><td><strong>76.96</strong></td><td>75.75</td><td>76.90</td></tr></tbody></table>
<table><caption>Table 24: Test accuracy of parallel TRGL and VanGL with 2 modules, compared to DGL and DDG from [Belilovsky et al., 2020] and InfoPro from [Wang et al., 2021] that also split a ResNet-110 in 2 module-wise-parallel-trained modules, on CIFAR 10.</caption><thead><tr><th>Parallel VanGL (ours)</th><th>Parallel TRGL (ours)</th><th>DGL</th><th>DDG</th><th>InfoPro</th></tr></thead><tbody><tr><td>94.01 ± .17</td><td><b>94.05</b> ± .18</td><td>93.50</td><td>93.41</td><td>93.58</td></tr></tbody></table>

<table><caption>Table 25: Average last block test accuracy and 95% confidence interval of 20-1 ResNet (32 filters, fixed encoder, same classifier) over 20/50 runs on MNIST with different train sizes and parallel (Par) TRGL and VanGL, compared to E2E.</caption><thead><tr><th>Train size</th><th>Par VanGL</th><th>Par TRGL</th><th>E2E</th></tr></thead><tbody><tr><td>60000</td><td>99.07 ± .04</td><td><b>99.08</b> ± .04</td><td>99.30 ± .03</td></tr><tr><td>30000</td><td>98.90 ± .05</td><td><b>98.93</b> ± .06</td><td>99.22 ± .03</td></tr><tr><td>12000</td><td>98.52 ± .06</td><td><b>98.59</b> ± .06</td><td>98.96 ± .06</td></tr><tr><td>6000</td><td>98.05 ± .09</td><td><b>98.16</b> ± .07</td><td>98.62 ± .06</td></tr><tr><td>1500</td><td>96.34 ± .12</td><td><b>96.91</b> ± .07</td><td>97.19 ± .08</td></tr><tr><td>1200</td><td>95.80 ± .12</td><td><b>96.58</b> ± .09</td><td>96.88 ± .09</td></tr><tr><td>600</td><td>91.35 ± .99</td><td><b>95.16</b> ± .15</td><td>95.30 ± .17</td></tr><tr><td>300</td><td>89.81 ± .73</td><td><b>92.86</b> ± .24</td><td>92.87 ± .28</td></tr><tr><td>150</td><td>81.84 ± .122</td><td><b>87.48</b> ± .42</td><td>87.82 ± .59</td></tr></tbody></table>

<table><caption>Table 26: Average last block test accuracy and 95% confidence interval of 20-1 ResNet (100 filters, fixed encoder, same classifier) over 10 runs on CIFAR10 with different train sizes and parallel (Par) TRGL and VanGL, compared to E2E.</caption><thead><tr><th>Train size</th><th>Par VanGL</th><th>Par TRGL</th><th>E2E</th></tr></thead><tbody><tr><td>50000</td><td>85.98 ± .28</td><td><b>86.02</b> ± .26</td><td>93.11 ± .19</td></tr><tr><td>25000</td><td>80.94 ± .25</td><td><b>81.09</b> ± .32</td><td>89.10 ± .29</td></tr><tr><td>10000</td><td>72.49 ± .46</td><td><b>73.01</b> ± .31</td><td>80.52 ± .46</td></tr><tr><td>5000</td><td>62.31 ± .54</td><td><b>64.06</b> ± .57</td><td>69.44 ± .88</td></tr><tr><td>500</td><td>38.61 ± .47</td><td><b>41.44</b> ± .44</td><td>40.40 ± .60</td></tr></tbody></table>

Figures 18 and 42 suggest that our regularization helps most when looking at the accuracy of the last block. We confirm this by including in Table 27 the accuracy achieved by the best block for the same experiment as in Table 11 and we notice a more important improvement from the regularization in the accuracy of the last block than in the accuracy of the best block. We also observe that with the regularization the difference between the accuracy of the last block and that of the best block is smaller than without the regularization. We further confirm this through the following experiment. As the network gets deeper (50 blocks trained for 10 epochs each sequentially), we expect training it block-wise to become more difficult, and indeed the improvement from the regularization is slightly larger than usual when looking at the accuracy of the last block for both sequential training methods (Table 28). We also include in Table 29 results from block-wise training of ResNeXt-50-32×4d [Xie et al., 2017], which turns out to be difficult to train block-wise.