PaperlayerView exact revision on arXiv

Source-linked AI summary

Running error bounds in finite element kernels

Michal Habera, Paul T. Kühner, Matteo Croci, Andreas Zilian

arXiv:2609.37844v1math.NAcs.MScs.NA

TL;DR

Rounding errors can undermine finite element accuracy, while their accumulation in automatically generated kernels is difficult to analyze in advance. This paper implements running error analysis to estimate errors alongside assembled values; tests detect instability and mesh-related error growth, with modest overhead in most cases.

  • Problem

    Rounding errors can limit finite element accuracy and solver convergence, but their accumulation in automatically generated and compiled kernels is difficult to analyze a priori.

    Source: §1 · p. 1

  • Method

    The framework applies running error analysis to generated finite element kernels, computing a posteriori rounding-error estimates alongside assembled values.

    Source: §1 · p. 1

  • Results

    In most tested cases, the two error modes show a slowdown of up to 3.5× versus plain assembly, while detecting Neo-Hooke instability and near-degenerate-mesh error growth.

    Source: §5 · p. 14

  • Takeaways & Limitations

    The estimates can diagnose rounding-error sensitivity in the tested finite element assemblies, including instability and mesh-conditioning effects.

    Source: §5 · p. 14

    The “exact” mode does not fully account for subnormal values and underflow in multiplication and division, which can invalidate residual exactness.

    Source: §3 · p. 6

Abstract

from arXiv · show

Rounding errors in finite element computations can lead to a complete loss of accuracy, stalled convergence, and incorrect results. Moreover, the effects of rounding errors accumulated within automatically generated and compiled kernels are difficult to analyze a priori. We present the first software framework for automated rounding error estimation within finite element kernels. The proposed methodology is based on an a posteriori technique called Running Error Analysis (REA), where a forward error estimate is automatically computed concurrently with the value. An open-source implementation is provided for the FEniCS Form Compiler (FFCx), based on a C++ backend for generating type-generic templated kernels over a custom arithmetic type that tracks both the value and its error estimate. We demonstrate REA on two examples. First, we use it to detect catastrophic cancellation in the assembly of a Neo-Hooke hyperelastic model in the small deformation regime. A series expansion circumvents the cancellation problem and the computed error estimates show this. Second, we study the assembly of the Laplace operator on a near-degenerate mesh. We demonstrate that our REA implementation typically incurs only 2-4x performance overhead. Applications of this work include robust reduced-precision computations in embedded systems, numerical debugging of new, possibly ill-conditioned or unstable PDE formulations, and guiding the design of mixed-precision kernels.

1 Introduction

The paper introduces automated a posteriori rounding-error estimation for finite element kernels, motivated by the difficulty of analyzing errors in generated code and the risks of reduced precision. Its examples show the estimates detecting instability in a hyperelastic model and error growth on near-degenerate meshes.

1 Introduction

  • Rounding-error analysis matters for reduced-precision computation, while finite-strain models near the undeformed state and ill-shaped elements can heighten sensitivity to roundoff.The paper notes that such errors can affect constitutive evaluations, solver convergence, and decisions near contact or yield thresholds.Source: §1 · p. 2
  • Prior FE-kernel work by Croci and Wells (2024) provides theoretical a priori bounds, whereas this work focuses on practical a posteriori estimates.The cited prior work covers basis and geometry evaluation, gradient transformations, quadrature accumulation, and local-to-global assembly.Source: §1 · p. 3
  • The paper presents the first software framework for automated a posteriori rounding-error estimation of finite element computations.It combines code generation with running error analysis to compute bounds alongside assembled values.Source: §1 · pp. 1, 3
  • The implementation integrates a custom C++ arithmetic type with FFCx-generated templated kernels and FEniCSx assemblers, providing “worst” bounds and “exact” signed estimates.In most tested Neo-Hooke cases, the slowdown is at most 3.5× relative to plain assembly.Source: §1 · p. 3
  • Both error modes detect instability in naïve Neo-Hooke energy evaluation despite bounded conditioning, and capture improved accuracy from a series-expansion implementation.The instability occurs in the small-deformation regime.Source: §1 · p. 3
  • For near-degenerate meshes, “worst” and “exact” estimates match conditioning-predicted asymptotic error growth, with “exact” capturing milder growth from roundoff cancellations.The “worst” mode captures worst-case growth.Source: §1 · p. 3

2 A posteriori rounding error analysis

The analysis derives a posteriori error estimates by separating local rounding error from propagated input error under a floating-point model. It covers both nonnegative first-order running bounds and signed estimates, which retain error cancellation information but are not certified bounds.

2.1 Standard floating-point model and notation

  • The floating-point model bounds each basic arithmetic operation’s relative rounding error by machine epsilon, excluding underflow and overflow.The model also applies to number representation when the value lies in the floating-point range.Source: §2.1 · p. 4

2.2 Running error analysis

  • Running Error Analysis (REA) estimates forward error by separating the current operation’s local rounding error from error propagated through its inputs.The resulting nonnegative upper bound is approximate and first-order in machine epsilon and the input-error bounds.Source: §2.2 · pp. 4–5
  • The computable running bound evaluates derivatives at the computed inputs and stores the estimate alongside the value for propagation into later operations.For twice continuously differentiable operations, this derivative substitution changes the estimate only by terms of second order in the input errors.Source: §2.2 · p. 5

2.3 Signed error estimates

  • Signed estimates combine a higher-precision reference value with signed first-order error propagation, preserving cancellation that absolute-value bounds discard.Error-free transformations can compute local rounding residuals exactly for supported operations; multiplication options include TwoProductFMA when FMA is available or TwoProduct at greater floating-point cost.Source: §2.3 · p. 5
  • Neither the signed estimates nor the running bounds are certified bounds, and a large estimated error may indicate ill-conditioning or numerical instability.The general best-case lower bound is zero, so automatically obtaining more informative lower bounds is nontrivial.Source: §2.3 · p. 5
  • Both error modes are implemented in a common arithmetic type and integrated into generated FEniCSx kernels and finite element assembly for comparison of accuracy and computational cost.Source: §2.3 · p. 5

3 Implementation for automatically generated FE kernels

The implementation generates templated C++ finite-element kernels so they can run with a custom arithmetic type that tracks values and errors. The type supports two estimation modes and is carried through FEniCSx assembly.

Custom backend

  • A custom FFCx backend generates templated C++ kernels, allowing compilation with a custom type that overloads arithmetic to track running error.The generated kernel evaluates basis functions and geometry, builds the integrand evaluation tree, and accumulates into tensor A.Source: §3 · p. 6

Custom C++ type

  • running_error_t stores an N-bit value followed by an N-bit error bound or signed estimate, for a total size of 2N bits.WORST selects the REA bound, while EXACT selects the higher-precision error estimate.Source: §3 · p. 6
  • The exact mode uses error-free transformations for addition, subtraction, and multiplication, but subnormal values and underflow can invalidate residual exactness in multiplication and division.The implementation notes that scaling and squeezing techniques could mitigate these range limitations.Source: §3 · p. 6
  • The custom type addresses implementation needs including control of binary layout for C++/Python interoperability and support for 16-bit floating-point values.The implementation uses a NumPy scalar representation with views of the value and error because structured conversion is unsupported in the described interface.Source: §3 · p. 7

Implementation workflow

  • The implementation compiles generated kernels and assembly loops for running_error_t<T, Mode>, then assembles scalars, vectors, or matrices in that type.The workflow uses FFCx, dolfiny, and dolfinx; assembly returns values and corresponding error bounds or estimates.Source: §3 · p. 7

4 Examples

Two examples test running error estimates: Neo-Hooke strain-energy evaluation at small deformation and Laplace assembly on a near-degenerate mesh. The estimates detect cancellation in the former and reflect ill-conditioning in the latter.

4.1 Hyperelastic energy of the Neo-Hooke model

  • For small deformations, the unstable expression’s relative error scales as Π^-2, while the stable expansion’s error remains nearly flat across deformation scales.The error estimates are consistent with a bounded condition number for the stable expression.Source: Figure 7 · p. 10, §4.1 · p. 10
  • The strain-energy evaluation is bounded in condition number as deformation decreases under the stated assumptions, despite the naïve expression’s numerical instability.The bound is independent of deformation scale on a fixed mesh.Source: §4.1 · pp. 10–11
  • Error estimation slows assembly by up to 3.5× in most tested cases, but the float64 unstable-expression exact mode reaches 22× because its reference logarithm is unoptimized.The per-cell overhead varies mildly with cell count and does not change linear assembly scaling.Source: §4.1 · p. 11, Table 2 · p. 12

4.2 Laplace operator on a near-degenerate mesh

  • For the float32 near-degenerate mesh, the exact mode gives much tighter row-wise error estimates than the worst-case bound.The mesh visualization shows the matrix entries and row-wise relative errors.Source: §4.2 · p. 13, Figure 9 · p. 13
  • For the Laplace operator, relative rounding errors scale as ∆^-1 in both estimation modes across float16, float32, and float64; the exact-mode error is around machine epsilon when ∆>0.1.For float16, data for ∆<5 × 10^-3 are excluded because of underflow or overflow.Source: Figure 10 · p. 13, §4.2 · p. 14
  • The measured ∆^-1 growth in the assembled matrix’s relative forward error indicates loss of conditioning as the mesh approaches degeneracy, consistent with the cited theoretical bounds.The reported bounds relate geometry-tensor error to the Jacobian condition number.Source: §4.2 · p. 14

5 Conclusion

The conclusion presents a running-error framework for generated finite-element kernels and evaluates it on hyperelastic and near-degenerate Laplace examples. The examples show how the estimates capture rounding-error behavior, while the implementation remains limited in its precision control.

5 Conclusion

  • The framework uses running_error_t to track error in “worst” bound and “exact” signed-estimate modes, passing that type to generated FFCx kernels linked into FEniCS assemblers.The “exact” estimate uses a higher-precision format and error-free transformations.Source: §5 · p. 14
  • In most cases, both error modes slow assembly by up to 3.5× compared with plain assembly.The authors describe this overhead as small enough to make the methods usable in practice.Source: §5 · p. 14
  • For Neo-Hooke energy at small deformations, both modes detect instability in the naïve implementation, while a series expansion improves accuracy.The underlying problem has a bounded condition number, and both modes indicate the series expansion’s improved accuracy.Source: §5 · p. 14
  • For the near-degenerate Laplace mesh, estimated rounding errors grow as distance to degeneracy decreases, with both modes matching the conditioning-based asymptotic scaling.The “worst” mode ignores error cancellations; the “exact” mode estimates smaller errors associated with stochastic roundoffs.Source: §5 · p. 14
  • FFCx’s type controls prevent choosing precision independently for geometry operations or different parts of a bilinear form.For example, geometry evaluation cannot use low-precision inputs with high-precision accumulation, and separate form components cannot be assigned different precisions.Source: §5 · p. 14
  • Future work could address low-precision range effects through scaling and squeezing, or dynamically choose operation precision from running-error estimates.Source: §5 · p. 14

A Proof of Lemma 1

The proof estimates perturbations in the Neo-Hooke strain energy using small-deformation expansions, an inverse inequality, and a non-rigid-body-modes assumption. It bounds the energy below by a quadratic term for sufficiently small deformation.

A Proof of Lemma 1

  • The proof first derives perturbations of the trace and determinant and then uses them to compute the strain-energy perturbation.Source: A Proof of Lemma 1 · p. 15
  • On a fixed cell, the proof applies the finite-element inverse inequality to P1 vector fields and sets H, S, and M from the displacement gradient and mesh parameters.It then invokes the non-rigid-body-modes assumption.Source: A Proof of Lemma 1 · p. 15
  • For sufficiently small ∇u, the deformation-gradient difference is bounded by the symmetric gradient plus a quadratic remainder.The displayed expansion is F − F^-T = ∇u + ∇u^T + O(∥∇u∥^2), bounded by C∥sym(∇u)∥ + O(∥∇u∥^2).Source: A Proof of Lemma 1 · p. 15
  • The quadratic energy expansion yields W ≥ c∥S∥^2 for sufficiently small ∥H∥, after its remainder is absorbed into the quadratic term.Source: A Proof of Lemma 1 · p. 15
Loading 2609.37844v1…