# Mixed precision and local error in ordinary differential equations

Arsène Marzorati, Mouhamad Al-Sayed Ali, Samuel Bernard, Jonathan Rouzaud-Cornabas

► To cite this version:

Arsène Marzorati, Mouhamad Al-Sayed Ali, Samuel Bernard, Jonathan Rouzaud-Cornabas. Mixed precision and local error in ordinary differential equations. Congrès National d'Analyse Numérique 2024, May 2024, Ile de Ré, France. hal-04953537

HAL Id: hal-04953537

https://inria.hal.science/hal-04953537v1

Submitted on 18 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Mixed precision and local error for ODE solvers
CANUM - 2024

M. Al-Sayed Ali

S. Bernard

A. Marzorati

J. Rouzaud-Cornabas

May 30, 2024

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

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

Institut
Camille
Jordan

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

INSTITUT NATIONAL
DES SCIENCES
APPLIQUÉES
CYM
■ Computational challenge:

Solving high dimensional and complex ODE systems for modelling biological systems.

■ General model:

$$
\dot{X}_i = F_i(X) = H_i(X_i) + \frac{1}{N} \sum_{j=1}^{N} G_{ij}(X_i, X_j) \quad i \in \{1, \dots, N\}, \ X_i \in \mathbb{R}^d.
$$
Computational limits (Benchmark: Kuramoto model)

<figure><img src="image_5.png" /></figure>
# Mixed precision: a trade-off ?

## Definition:

Use several numerical arithmetic precisions inside one computational tool.

## Benefits:

Computational acceleration, less memory needed, better error control

## Tools:

Linear Algebra, Machine Learning

## Applications:

Physics, Meteorology
Mixed-Precision and modern chips

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

<figure><img src="image_7.png" /></figure>
Due to the lack of transparency in Matlab computations, we choose to force the arithmetic precision before each operation to ensure the corresponding accuracy. We use the following theoretical assumption. For an arithmetic operation:

$$T_D \sim 2T_S.$$

$$T_D \sim \alpha T_{MP}, \alpha \in [1, 2].$$

We count the number of evaluations performed in one specific arithmetic precision.
An ODE problem: $\dot{X} = f(t, X), t \in [0, T]$.

1) **Step time:** $h = \frac{T}{N}$ and setting $X_i = X(t_i) = X(ih)$.

2) **Get $X_{n+1}$:** ($X_n$ known).

The stage(s): $k_i = f(t_n + c_i h, X_n + h \sum_{j=1}^{i-1} a_{i,j} k_j), \forall i \in \{1, ..., p\}$.

$$X_{n+1} = X_n + h \sum_{i=1}^{p} b_i k_i.$$

3) **Order of the scheme:** truncation error.
# Solver ODE23: Adaptive scheme

* **ODE23:** Combination of RK3 with 3 stages and RK2 with 4 stages (first three stage are shared).

* **Adaptive scheme:** Error evaluation (RK2 vs RK3) and step validation (with the relative tolerance) modifying the step size ($h \to h_n$).

* **First Same As Last:** Last stage at step $n$ is used as First stage at step $n+1$.

* **Bonus:** Already implemented in Matlab.
■ For function evaluation, 3 possibilities can be chosen :

$$ \dot{X}_i = H_i(X_i) + \frac{1}{N} \sum_{j=1}^{N} G_{ij}(X_i, X_j). $$

■ For each stage ($p$-stages, 3 in our case) a different *cocktail* can be chosen.

<table><thead><tr><th rowspan="2">RK-Stage</th><th colspan="3">Precision</th></tr><tr><th>H<sub>i</sub></th><th>&Sigma;<sup>N</sup><sub>j=1</sub></th><th>G<sub>j</sub></th></tr></thead><tbody><tr><td>K<sub>2</sub></td><td>S</td><td>D</td><td>S</td></tr><tr><td>K<sub>3</sub></td><td>S</td><td>D</td><td>S</td></tr><tr><td>K<sub>4</sub></td><td>D</td><td>D</td><td>S</td></tr></tbody></table>
Model equations

$$ \left\{ \begin{array}{l} \frac{dx_i}{dt} = y_i + \frac{1}{N} \sum_{j=1}^{N} (x_j - x_i) , \quad \forall i \in \{1, \dots, N\} \\ \\ \frac{dy_i}{dt} = -x_i \end{array} \right. $$

Why ?

* Appropriate structure

* Analytic solution

* Biological application: biological rhythm
Model equations

$$ \frac{dx_i}{dt} = \omega_i + \frac{1}{N} \sum_{j=1}^{N} K \sin(x_j - x_i), \quad \forall i \in \{1, \dots, N\}. $$

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

Why ?

* Dense literature

* Non-linear interaction term (sine)

* Biological application: neuroscience
Benchmark 3: Circadian clock

Model equations

$$
\begin{cases}
\dot{x}_i = \left( \frac{k_0 \theta^{\frac{h}{a}}}{\theta^{\frac{h}{a}} + y_i^{\frac{h}{a}}} x_i - k_1 \right) x_i + \frac{1}{N} \sum_{j=1}^{N} \frac{k_0 \theta^{\frac{h}{a}}}{\theta^{\frac{h}{a}} + y_j^{\frac{h}{a}}} K \arctan(x_j - x_i), \\
\dot{y}_i = k_2 (x_i - y_i), \\
\dot{u}_i = u_i \left(1 - \frac{u_i^2}{3}\right) - v_i + I_0 \left(1 - \frac{k_3^2}{k_3^2 + x_i^2}\right), \\
\dot{v}_i = \epsilon (u_i + b - C v_i).
\end{cases}
$$

With $X_i = (x_i, y_i, v_i, u_i)^T \in \mathbb{R}^4$ and $X \in \mathbb{R}^{4 \times N}$

■ Why ?

● Real (simplified) model <fnref n="1" />

● Non-linear interaction term (arctan)

<footnote n="1">¹El Cheikh R., Bernard S., & El Khatib N. (2017). A multiscale modelling approach for the regulation of the cell cycle by the circadian clock. J Theor Biol, 426, 117-125.</footnote>
## Error computation

Choose one solution as reference: $X_{ref}(t)$.

At final time, compute:

$$ ||X_{ref}(t_f) - X(t_f)||_\alpha = \frac{||X_{ref}(t_f) - X(t_f)||_2}{\sqrt{N}} $$

## In the next slides:

* The reference solution is computed with ODE45 of Matlab with a relative tolerance of $10^{-9}$.

* The values are averaged over all the tests completed by all the solvers.
# Why not switch completely to low precision ?

<figure><img src="image_9.png" /><figcaption>Circadian clock ($N = 1500$)</figcaption></figure>

<figure><img src="image_10.png" /></figure>
Relative Tolerance: $10^{-3}$

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

Relative Tolerance: $10^{-6}$

<figure><img src="image_12.png" /></figure>
Relative Tolerance: $10^{-7}$

<figure><img src="image_13.png" /></figure>
Relative Tolerance: $10^{-3}$

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

Relative Tolerance: $10^{-6}$

<figure><img src="image_15.png" /></figure>
Relative Tolerance: $10^{-7}$

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

Relative Tolerance: $10^{-8}$

<figure><img src="image_17.png" /></figure>
## Conclusion

* Mixed precision is necessary to benefit from the recent development of chips.

* A good trade-off for solving high-dimensional systems with sufficient accuracy with lower numerical precision.

## Outlook

* Develop numerical methods that select the appropriate numerical precision from the many possibilities.

* Performance measurements on different architectures.

**Thank you for your listening!**
1 Ackmann J., Dueben P.D., Palmer T.N. and Smolarkiewicz P.K., "Mixed-Precision for Linear Solvers in Global Geophysical Flows", J Adv Model Earth Syst (2022).

2 Bogacki P. and Shampine L.F., "A 3 (2) pair of Runge-Kutta formulas." Applied Mathematics Letters 2.4 (1989): 321-325.

3 Croci M. and de Souza G.R., "Mixed-precision explicit stabilized Runge-Kutta methods for single- and multi-scale differential equations." J. Comput. Phys. (2022)

4 Higham N.J. and Mary T. "Mixed precision algorithms in numerical linear algebra." Acta Numerica 31 (2022): 347-414.

5 Haidar A., Tomov S., Dongarra J. and Higham N.J., "Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers." SC18: Int. Conf. High Perform. Comput. Netw. Storage Anal. IEEE (2018).

6 Choquette J. and Wish G., "Nvidia a100 gpu: Performance & innovation for gpu computing." IEEE Hot Chips 32 Symposium (HCS). IEEE Computer Society, (2020).
7. Burnett B., Gottlieb S., Grant Z.J., and Heryudono A., "Performance evaluation of mixed-precision Runge-Kutta methods. IEEE High Performance Extreme Computing Conference (HPEC) (2021).

8. Blanchard P., Higham N.J. and Mary T., "A class of fast and accurate summation algorithms, MIMS EPrint 2019.6." Manchester Institute for Mathematical Sciences, The University of Manchester, UK (2019).

9. Higham, N.J and Pranesh S., "Simulating low precision floating-point arithmetic", SIAM J Sci Comput (2019).

10. Hayford J., Goldman-Wetzler J., Wang E. and Lu L., "Speeding up and reducing memory usage for scientific machine learning via mixed precision", arXiv preprint (2024).

11. El Cheikh R., Bernard S., and El Khatib N., "A multiscale modelling approach for the regulation of the cell cycle by the circadian clock". J Theor Biol (2017), 426, 117-125.

12. Palmer T.N., "More reliable forecasts with less precise computations: a fast-track route to cloud-resolved weather and climate simulators?" Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372.2018 (2014).
Linear oscillators (Relative Tolerance: $10^{-3}$, N=1000)

<figure><img src="image_18.png" /><figcaption>Local error:</figcaption></figure>

Time

<figure><img src="image_19.png" /><figcaption>Step size (power 4):</figcaption></figure>

Time
Linear oscillators (Relative Tolerance: $10^{-6}$, N=1000)

Local error:

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

Time

Step size (power 4):

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

Time
Tolerance and step validation

At each step $X_1^{(n)}$ and $X_2^{(n)}$ are computed, $X_2^{(n-1)}$ is the solution at the previous step.

At final time ($t_f$), compute:

$$err := \left\| (X_1^{(n)} - X_2^{(n)}) / \max \left( |X_2^{(n)}|, |X_2^{(n-1)}|, \frac{Ab}{Rel} \right) \right\|_{\infty}.$$

'./' division term by term.
<table><thead><tr><th>Parameters</th><th>Value(s)</th></tr></thead><tbody><tr><td>Tests number</td><td>1000</td></tr><tr><td>N (number of oscillators)</td><td>{100, 500, 1000, 2000, 5000}</td></tr><tr><td>Tol - Ab</td><td>{10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-3</sup>}</td></tr><tr><td>Tol - Rel</td><td>{10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-3</sup>}</td></tr><tr><td>T<sub>f</sub> (Final time)</td><td>10&pi;</td></tr><tr><td>Initial conditions</td><td>[0, 2&pi;]</td></tr></tbody></table>
<table><thead><tr><th>Parameter</th><th>Value(s)</th></tr></thead><tbody><tr><td>Number of tests</td><td>5000</td></tr><tr><td>N (Number of oscillators)</td><td>{100, 500, 1000, 2000, 3000, 5000}</td></tr><tr><td>Tol - Ab</td><td>{10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-3</sup>}</td></tr><tr><td>Tol - Rel</td><td>{10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-3</sup>}</td></tr><tr><td>σ</td><td>[0, 1]</td></tr><tr><td>ω<sub>i</sub></td><td>[-σ, σ]</td></tr><tr><td>K (coupling coefficient)</td><td>[0, 3]</td></tr><tr><td>T<sub>f</sub> (Final time)</td><td><sup>4π</sup> / med(|W|)K + 0.001</td></tr><tr><td>Initial conditions</td><td>[0, 2π]</td></tr></tbody></table>
<table><thead><tr><th>Parameter</th><th>Value(s)</th></tr></thead><tbody><tr><td>Number of tests</td><td>1000</td></tr><tr><td><i>N</i> (Number of oscillators)</td><td>{100, 300, 700, 1000, 1500}</td></tr><tr><td><i>T</i><sub>ol</sub> &ndash; <i>Ab</i></td><td>{10<sup>&minus;8</sup>, 10<sup>&minus;7</sup>, 10<sup>&minus;6</sup>, 10<sup>&minus;3</sup>}</td></tr><tr><td><i>T</i><sub>ol</sub> &ndash; <i>Rel</i></td><td>{10<sup>&minus;8</sup>, 10<sup>&minus;7</sup>, 10<sup>&minus;6</sup>, 10<sup>&minus;3</sup>}</td></tr><tr><td><i>K</i> (coupling coefficient)</td><td>{0.001, 0.1, 1, 10}</td></tr><tr><td><i>l</i><sub>0</sub></td><td>{0.228249, 1.5, 10}</td></tr><tr><td><i>T</i><sub>f</sub> (Final time)</td><td>3 &times; <i>T</i><sub>cycle</sub> = 72</td></tr></tbody></table>
Accuracy: Circadian clock (Relative Tolerance: $10^{-8}$)

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