# Numerical Solution of Mixed-Dimensional PDEs Using a Neural Preconditioner

Nunzio Dimola, Nicola R. Franco, Paolo Zunino

► To cite this version:

Nunzio Dimola, Nicola R. Franco, Paolo Zunino. Numerical Solution of Mixed-Dimensional PDEs Using a Neural Preconditioner. 2025. hal-05056688

HAL Id: hal-05056688

https://hal.science/hal-05056688v1

Preprint submitted on 5 May 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
# Numerical Solution of Mixed-Dimensional PDEs
Using a Neural Preconditioner

Nunzio Dimola¹, Nicola Rares Franco¹, Paolo Zunino¹,\*

¹MOX, Department of Mathematics, Politecnico di Milano, P.zza Leonardo da Vinci 32, Milan, 20133, Italy

## Abstract

Mixed-dimensional partial differential equations (PDEs) are characterized by coupled operators defined on domains of varying dimensions and pose significant computational challenges due to their inherent ill-conditioning. Moreover, the computational workload increases considerably when attempting to accurately capture the behavior of the system under significant variations or uncertainties in the low-dimensional structures such as fractures, fibers, or vascular networks, due to the inevitable necessity of running multiple simulations. In this work, we present a novel preconditioning strategy that leverages neural networks and unsupervised operator learning to design an efficient preconditioner specifically tailored to a class of 3D-1D mixed-dimensional PDEs. The proposed approach is capable of generalizing to varying shapes of the 1D manifold without retraining, making it robust to changes in the 1D graph topology. Moreover, thanks to convolutional neural networks, the neural preconditioner can adapt over a range of increasing mesh resolutions of the discrete problem, enabling us to train it on low resolution problems and deploy it on higher resolutions. Numerical experiments validate the effectiveness of the preconditioner in accelerating convergence in iterative solvers, demonstrating its appeal and limitations over traditional methods. This study lays the groundwork for applying neural network-based preconditioning techniques to a broader range of coupled multi-physics systems.

**Keywords** mixed-dimensional PDEs, finite element approximation, scientific machine learning, preconditioning

## 1 Introduction

Research on mixed-dimensional partial differential equations (PDEs) is advancing to address challenges in complex multiscale applications in fields such as geophysics, biomechanics, and neuroscience. These disciplines often involve thin structures with high aspect ratio, including structural mechanics and geomechanics [30, 24, 14], microcirculation and perfusion modeling [40, 34], axons in neural networks [7], and environmental flows [18]. Modeling these interactions with traditional 3D meshes would require massive computational resources, making direct simulations impractical. Mixed-dimensional PDEs, which enable lower-dimensional representations to be embedded within higher-dimensional domains, present a promising alternative. However, the computational effort associated with these models rises significantly when trying to capture the behavior of the system for varying configurations or interactions within the low-dimensional structures, the latter representing, e.g., fracture networks in geological formations, fiber distributions in structural materials, or vascular and neural networks in biological systems. Similarly, multiple simulations are needed when the structure of the low-dimensional problem is uncertain or unknown. In all such cases, despite the simplifications provided by mixed-dimensional PDEs, the overall computational demand remains high due to the cumulative cost of repeated individual simulations. This growing need underscores the importance of developing efficient solvers and preconditioners capable of handling the repeated solutions essential to representing the statistical variability of intricate networks in real-world applications.

Today, the discretization of mixed-dimensional PDEs is established, providing computationally efficient schemes suitable for applications with intricate networks [5, 16, 22], with rigorous error estimates ensuring

<footnote>\*Corresponding author</footnote>

<footnote>Email addresses: nunzio.dimola@polimi.it, nicolarares.franco@polimi.it, paolo.zunino@polimi.it</footnote>
accuracy. The well-established literature on block and operator preconditioning techniques, see, e.g., [28], has been instrumental in the development of computational solvers for such problems. [21, 23] have paved the way for the efficient resolution of discrete mixed-dimensional problems. Along this direction, we also mention the work of [8] on algebraic multigrid methods for metric-perturbed coupled problems, which has significantly advanced the understanding of block preconditioning techniques in mixed-dimensional PDEs. Parallel to this, other researchers have also explored approximate algorithms that can enhance scalability for large systems: see, e.g. [12].

The goal of this study is to develop a novel preconditioning strategy that addresses the twofold challenges inherent in mixed-dimensional PDEs: the ill-conditioned behavior of these problems and the need for adaptability across diverse geometric configurations of the low-dimensional subproblems. We aim to create a preconditioner that not only mitigates the ill-conditioning of the system, but also maintains this improved conditioning across a variety of configurations for 1D or other low-dimensional structures. To accomplish this, the preconditioner is designed as a nonlinear operator rather than a conventional matrix, enabling it to process the entire solution manifold effectively. Our method takes advantage of the latest advances in the approximation capabilities of neural networks and learning approaches related to machine learning [20, 6, 26]. These advancements have sparked active research in the creation of novel solvers and preconditioners tailored for different types of problem [15, 19, 17, 29, 2], this being a non-exhaustive list of examples supporting the potential of neural network-based preconditioners as a viable alternative for solving large, sparse linear systems across various scientific and engineering domains.

Building upon these studies, we develop a nonlinear preconditioner designed to efficiently manage the numerical challenges associated with the complex structure of mixed-dimensional equations. By incorporating a shape descriptor into the learning framework, the preconditioner can generalize across different configurations of the 1D graph or low-dimensional structures without the need for retraining. In order to demonstrate the efficacy of the proposed approach, we focus our attention on a mixed-dimensional 3D-1D model problem, which serves as a representative template for a broader class of coupled systems with similar characteristics. The results show that the learned preconditioner significantly accelerates the convergence of the iterative Generalized Minimum Residual (GMRES) solver, even as the complexity of the problem increases due to changes in the 1D graph's topology or the size of the 3D domain. Unlike traditional preconditioning techniques, which typically require problem-specific tuning or reconfiguration when the underlying system changes, our method maintains consistent performance across diverse problem settings and mesh resolutions. Moreover, the proposed preconditioning approach can be easily generalized, making it a versatile and robust tool for a wide array of mixed-dimensional problems beyond the specific examples addressed in this study.

This manuscript is organized as follows. First, we introduce the mathematical formulation of the 3D-1D mixed-dimensional PDE and highlight the challenges associated with solving such systems. Then, we present the proposed preconditioner learning framework, discussing both the theoretical foundations and the practical implementation of the nonlinear preconditioner. Next, we apply the preconditioner to the 3D-1D problem, illustrating its construction and providing a detailed analysis of its efficiency and scalability. Finally, we validate the approach through extensive numerical experiments, comparing its performance against conventional preconditioning methods, and conclude with a discussion of the implications and potential future directions for extending this work.

## 2 A mixed-dimensional problem

We consider here a mixed-dimensional elliptic problem as a simple template for a wider class of coupled problems with similar characteristics. It will be our reference for comparisons, analyses, and considerations on preconditioner performances. Originally, we define the problem in a three-dimensional (3D) slender domain embedded in an external one, where we set partial differential equations coupled by suitable interface conditions. The mixed-dimensional model is obtained by a *dimensional reduction* strategy that represents the slender interior domain as a one-dimensional (1D) manifold, namely a metric graph. More precisely, we address a coupled problem defined by an exterior domain $\Omega \subset \mathbb{R}^3$ and an interior domain $\Sigma \subset \Omega$. We assume $\Sigma$ to be a generalized cylinder with a centerline $\Lambda$, the latter being a 1D domain with arc length parameter $s \in (0, S)$. We assume that the cylinder has a constant radius $\epsilon > 0$, so that if $\lambda : (0, S) \to \Lambda$ parameterizes the centerline, then $\Sigma = \{\lambda(s) + \varepsilon : s \in (0, S), \|\varepsilon\| \le \epsilon, \varepsilon \perp \tau(s)\}$ with $\tau(s)$ the tangent vector at $\Lambda$. We assume that the transverse diameter of $\Sigma$, which is equal to $2\epsilon$, is small compared to the diameter of $\Omega$.
The dimensionality reduction process consists of replacing $\Sigma$ with $\Lambda$ while keeping some information about its original 3D structure and its interaction with the external domain $\Omega$. Specifically, let $\Gamma := \partial\Sigma$ be the interface of the coupled problem in the full 3D representation and let $\mathcal{T}_\Lambda$ be the restriction operator from $\Omega$ to $\Lambda$. In [22] it is defined as the composition of the trace operator on $\Gamma$ combined with a projection operator from $\Lambda$ to $\Omega$ based on cross-sectional averages. Then, the continuous 3D-1D coupled problem can be formally formulated as follows,

$$ \left\{ \begin{array}{ll} -\nabla \cdot (k_{\Omega} \nabla u_{\Omega}) + \sigma_{\Omega} u_{\Omega} + 2\pi\epsilon (\mathcal{T}_{\Lambda} u_{\Omega} - u_{\Lambda}) \delta_{\Lambda} = 0, & \text{on } \Omega, \\ -d_{s}(k_{\Lambda} \partial_{s} u_{\Lambda}) + 2\pi\epsilon (u_{\Lambda} - \mathcal{T}_{\Lambda} u_{\Omega}) = 0, & \text{on } \Lambda, \\ -\nabla u_{\Omega} \cdot \mathbf{n} = 0, & \text{on } \partial\Omega, \\ -d_{s} u_{\Lambda} \cdot \mathbf{n} = 0, & \text{on } \partial\Lambda \setminus \partial\Lambda_{D}, \\ u_{\Lambda} = 1, & \text{on } \partial\Lambda_{D}. \end{array} \right. \quad (1) $$

The functions $u_\Omega, u_\Lambda$ are the coupled unknowns on the 3D and the 1D domains, respectively, and $k_\Omega, \sigma_\Omega, k_\Lambda$ and $2\pi\epsilon$ represent respectively the coefficients of the elliptic operators on the 3D domain, the 1D one and the coupling parameter. Here, $\partial\Lambda_D \subset \partial\Lambda$ is a given subset of the 1D domain holding a Dirichlet boundary condition. An analogous formulation holds when the low-dimensional structure is given by the union of multiple substructures $\Sigma = \cup_i \Sigma_i$ and $\Lambda = \cup_i \Lambda_i$, each being a generalized cylinder. In this case, $\Lambda$ becomes a 1D graph embedded within a 3D domain. Figure 1 shows an example of a solution to this mixed-dimensional problem (1D solution on the left, 3D solution in the center).

<figure><img src="image_2.png" /><figcaption>Figure 1: Solution of a 3D-1D coupled problem for a fixed graph geometry $\Lambda$. Left: plot of 1D solution $u_\Lambda$. Center: slices of 3D solution $u_\Omega$. Right: slices plot of the distance function $d(\Lambda)$ associated to the graph.</figcaption></figure>

## 2.1 Parametrization and discretization of the problem

As we mentioned, we are interested in solving multiple instances of the above problem for different configurations of the low-dimensional structure, including, for instance, changes in the topology, arc-length or density of the 1D graph $\Lambda$. To formalize this fact, it is useful to define a suitable parameter space, $\mathcal{P}$, that collects all configurations of interest. To this end, we first note that each 1D graph $\Lambda \subset \Omega$ is uniquely identified by its 3D distance function, $d(\Lambda) : \Omega \to \mathbb{R}$, which maps each point $x \in \Omega$ to its distance from the closest point in $\Lambda$: see the right panel of Figure 1 for a visual representation. This allows us to define the parameter space as a subset of the space of continuous functions on $\Omega$, that is, $\mathcal{P} := \{d(\Lambda)\}_\Lambda \subset \mathcal{C}(\Omega)$. This representation has several advantages. First, it can be easily transferred to the discrete setting by introducing a suitable discretization of the exterior domain $\Omega$ using, for example, Finite Elements (FE). Furthermore, it automatically enriches the parameter space by equipping it with a metrizable topology (the one induced by the supremum norm): this allows us to discuss, for instance, about the continuity of the 3D solution with respect to the 1D graph. Finally, it will be useful when discussing our neural preconditioning approach in Section 4.

We now come to the discretization of (1). To this end, we first address its variational formulation. Let $V := H^1(\Omega) \times H^1_{\partial\Lambda_D}(\Lambda)$ be the mixed-dimensional space where PDE solutions will be sought. Here, $H^1_{\partial\Lambda_D}(\Lambda) = \{w \in H^1(\Lambda) : w|_{\partial\Lambda_d} = 1\}$ accounts for Dirichlet boundary conditions, so that —up to translations— $V$ is
isomorphic to a Hilbert space. In the parametric setting, re-writing Eq. (1) in weak form corresponds to saying that, for every $\mu \in \mathcal{P}$ we aim to find $u^\mu \in V$ such that

$$a^\mu(u^\mu, v) = F^\mu(v), \quad \forall v \in V, \qquad (2)$$

with $a^\mu: V \times V \to \mathbb{R}$ and $F^\mu: V \to \mathbb{R}$ parameter-dependent operators. Specifically,

$$a^\mu(u, v) = a_\Omega(u_\Omega, v_\Omega) + a_\Lambda^\mu(u_\Lambda, u_\Omega) + b_\Lambda^\mu(\mathcal{T}_\Lambda u_\Omega - u_\Lambda, \mathcal{T}_\Lambda v_\Omega - v_\Lambda),$$

where, being $\mathcal{D}(s)$ the cross section of $\Sigma$ at arc-length $s$, and having set

$$(u, v)_\Omega := \int_\Omega u(x)v(x)dx, \quad (u, v)_{\Lambda, |\mathcal{D}|} := \int_\Lambda u(s)v(s)|\mathcal{D}(s)|ds, \quad (u, v)_{\Lambda, |\partial\mathcal{D}|} := \int_\Lambda u(s)v(s)|\partial\mathcal{D}(s)|ds,$$

we have,

$$a_\Omega(u, v) = (k_\Omega \nabla u, \nabla v)_\Omega + (\sigma_\Omega u, v)_\Omega, \qquad (3)$$

$$a_\Lambda^\mu(u, v) = (k_\Lambda d_s u, d_s v)_{\Lambda, |\mathcal{D}|}, \quad b_\Lambda^\mu(u, v) = 2\pi\epsilon(u, v)_{\Lambda, |\partial\mathcal{D}|}. \qquad (4)$$

We note that the adopted notation is slightly redundant because the parametric dependence of $a^\mu$ is only due to the fact that $a_\Lambda^\mu$ and $b_\Lambda^\mu$ have support in $\Lambda$.

In order to discretize (2), we consider its Galerkin projection onto a (broken) FE space $V_h = V_h^\Omega \times V_h^\Lambda$ of dimension $N_h = N_h^\Omega + N_h^\Lambda$, being $N_h^\Omega = \dim(V_h^\Omega)$ and $N_h^\Lambda = \dim(V_h^\Lambda)$, which is suitably chosen depending on the characteristics of the problem at hand. Assuming for simplicity a fully conformal approximation, given $V_h \subset V$ with $V_h^\Omega \subset H^1(\Omega)$, $V_h^\Lambda \subset H^1_{\partial\Lambda_D}(\Lambda)$, we aim to find $u_h^\mu \in V_h$ such that

$$a^\mu(u_h^\mu, v_h) = F^\mu(v_h) \quad \forall v_h \in V_h. \qquad (5)$$

From a discrete point of view, the problem (5) is equivalent to a (large) system of algebraic equations $A_h^\mu \mathbf{u}_h^\mu = \mathbf{F}_h^\mu$. Specifically, we have

$$\left\{ \begin{array}{c} A_h^\mu \\ \left( \begin{array}{cc} K_{h,00} & 0 \\ 0 & K_{h,11}^\mu \end{array} \right) + 2\pi\epsilon \left( \begin{array}{cc} M_{h,00}^\mu & M_{h,01}^\mu \\ M_{h,10}^\mu & M_{h,11}^\mu \end{array} \right) \\ K_h^\mu \end{array} \right\} \left\{ \begin{array}{c} \mathbf{u}_h^\mu \\ \left( \begin{array}{c} \mathbf{u}_{h,0}^\mu \\ \mathbf{u}_{h,1}^\mu \end{array} \right) \\ M_h^\mu \end{array} \right\} = \left\{ \begin{array}{c} \mathbf{F}_h^\mu \\ \left( \begin{array}{c} 0 \\ \mathbf{f}_{h,1}^\mu \end{array} \right) \end{array} \right\}, \qquad (6)$$

where $\mathbf{u}_h^\mu \in \mathbb{R}^{N_h}$ is the vector of degrees of freedom of the FE approximation; $K_{h,00}$ and $K_{h,11}^\mu$ represent the discretized bilinear forms $a_\Omega$ and $a_\Lambda^\mu$; the block matrix $M_h^\mu$ models the coupling enforced by the operator $b_\Lambda^\mu(\mathcal{T}_\Lambda u_\Omega - u_\Lambda, \mathcal{T}_\Lambda v_\Omega - v_\Lambda)$, and the forcing term $\mathbf{f}_{h,1}^\mu$ takes into account the non-homogeneous Dirichlet condition on the 1D graph.

As we discretize the differential problem, we also take the opportunity to discretize the parameter space, which we do by projecting—or interpolating— all 3D distance functions onto $V_h^\Omega$. This allows us to represent each parametric configuration $\mu \in \mathcal{P}$ through a finite-dimensional vector $\mathbf{d}_h^\mu \in \mathbb{R}^{N_h^\Omega}$, obtained by listing the nodal values of $\mu = d(\Lambda)$ at the degrees of freedom. In this way, the discrete/finite-dimensional parameter space becomes $\mathcal{P}_h = \{\mathbf{d}_h^\mu\} \subset \mathbb{R}^{N_h^\Omega}$.

**Remark 1.** As apparent from equations (3)-(4), we notice that problem (2) is symmetric positive definite. However, we will not exploit this property when defining the computational solver. This choice is motivated by two main reasons. The first one is that (2) is just a particular case of a general family of operators in which the bilinear form $a_\Omega(\cdot, \cdot)$ may also describe transport phenomena that violate symmetry. The second reason is that, as it will be clearer later on, learned techniques, such as our neural preconditioner, can hardly preserve symmetry, unless this constraint is explicitly built into the network architecture—a strategy that currently lacks a reliable and efficient implementation.
**Remark 2.** In what follows, we shall make the assumption that the 3D domain has been discretized using a uniform mesh, consisting either of tetrahedrons or bricks. This requirement will be crucial in Section 4 when implementing our neural network-based preconditioner as it will allow us to exploit the relationship between tensor-like structures and convolutional neural networks. Nevertheless, we emphasize that in the context of mixed-dimensional problems, where the primary challenge is representing a complex low-dimensional structure embedded in 3D, this restriction is acceptable. In fact, this requirement does not constrain the low-dimensional domain, which can still exhibit arbitrary geometric complexity.

## 2.2 A preconditioning strategy for the coupled problem

As we mentioned, solving (6) can be computationally intensive due to poor conditioning of the linear system. Here, we address this fact by using a block-preconditioning approach, specifically tailored for the block-structured matrix arising from the discretization of the coupled 3D-1D problem. The goal is to employ a right-preconditioner expressed as:

$$A_h^\mu Q_h^\mu \mathbf{z}_h^\mu = \mathbf{F}_h^\mu, \quad Q_h^\mu \mathbf{z}_h^\mu = \mathbf{u}_h^\mu, \qquad (7)$$

where $Q_h^\mu$ is the inverse of the block upper triangular part of $A_h^\mu$. Employing this right preconditioner to accelerate the convergence of an iterative method, see e.g. [37], involves conditioning the residual from $[\mathbf{z}_{h,0}^\mu, \mathbf{z}_{h,1}^\mu]^\top$ to $[\tilde{\mathbf{z}}_{h,0}^\mu, \tilde{\mathbf{z}}_{h,1}^\mu]^\top$ by solving the following linear system through a two-step back-substitution procedure,

$$\begin{bmatrix} (K_{h,00} + 2\pi\epsilon M_{h,00}^\mu) & 2\pi\epsilon M_{h,01}^\mu \\ 0 & (K_{h,11}^\mu + 2\pi\epsilon M_{h,11}^\mu) \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{z}}_{h,0}^\mu \\ \tilde{\mathbf{z}}_{h,1}^\mu \end{bmatrix} = \begin{bmatrix} \mathbf{z}_{h,0}^\mu \\ \mathbf{z}_{h,1}^\mu \end{bmatrix}. \qquad (8)$$

In general, the first step, which is to solve for $\tilde{\mathbf{z}}_{h,1}^\mu$, does not pose significant challenges. This is because the dimension of the discrete 1D problem is typically smaller than the one of the 3D problem, $N_h^\Lambda \ll N_h^\Omega$, making it more manageable to address using traditional techniques, such as direct solvers. In contrast, the second step, which is to solve for $\tilde{\mathbf{z}}_{h,0}^\mu$, can be computationally demanding and can severely hinder the applicability of the approach when considering many-query scenarios, where the linear system is to be solved for multiple instances of the model parameters. This step, in fact, is the one where the main challenges inherent to mixed-dimensional problems become apparent. Indeed, as noted in [8], the latter involves elliptic operators altered by a low-order semidefinite component, referred to as the *metric term*. This term, whose impact is modulated by a scalable intensity parameter $\epsilon$, can disrupt the structure and conditioning of the elliptic component, thereby complicating the numerical approximation of the coupled system. For these reasons, our objective is to develop a nonlinear preconditioner, $\mathcal{P}$, for the mixed-dimensional system defined on $\Omega$,

$$(K_{h,00} + 2\pi\epsilon M_{h,00}^\mu)\tilde{\mathbf{z}}_{0,h}^\mu = \mathbf{z}_{h,0}^\mu - 2\pi\epsilon M_{01}^\mu \tilde{\mathbf{z}}_{1,h}^\mu. \qquad (9)$$

This nonlinear preconditioner will take two arguments: one that is conformal to the 3D solution (as any classical preconditioner would do), and one that refers to the problem parameters, allowing the preconditioner to adjust depending on the scenario of interest. We shall write $\mathcal{P} = \mathcal{P}(\mathbf{v}, \mathbf{d})$ in the discrete setting and, occasionally, $\mathcal{P} = \mathcal{P}(v, \mu)$ when considering the more abstract continuous formulation. Using this notation, and having set

$$C_h^\mu := (K_{h,00} + 2\pi\epsilon M_{h,00}^\mu), \quad \mathbf{v}_h^\mu := \tilde{\mathbf{z}}_{0,h}^\mu, \quad \mathbf{b}_h^\mu := \mathbf{z}_{h,0}^\mu - 2\pi\epsilon M_{01}^\mu \tilde{\mathbf{z}}_{1,h}^\mu, \qquad (10)$$

the right-preconditioned family of problems reads $C_h^\mu \mathcal{P}(\mathbf{v}_h^\mu, \mathbf{d}_h^\mu) = \mathbf{b}_h^\mu$. Considering the specific structure of $\mathbf{F}_h^\mu$ in (6), we focus on making the preconditioner $\mathcal{P}$ efficient for a limited subset of acceptable right-hand sides within $\mathbb{R}^{N_h^\Omega}$. In particular, we are interested in solving the system for the forcing terms that arise from the solution of the 1D problem, namely all $\mathcal{B} = \{\mathbf{b}_h^\mu\}_{\mu \in \mathcal{P}}$ such that

$$\mathcal{B} := \{-2\pi\epsilon M_{01,h}^\mu \mathbf{x} : \mu \in \mathcal{P} \text{ and } \mathbf{x} \in \mathbb{R}^{N_h^\Omega} \text{ with } (K_{11}^\mu + 2\pi\epsilon M_{11,h}^\mu)\mathbf{x} = \mathbf{f}_h^\mu\}. \qquad (11)$$

As we anticipated, to guarantee generalization across various configurations of the 1D problem, it is essential to develop a suitable nonlinear preconditioner, $\mathcal{P}$. This involves carefully choosing an appropriate hypothesis space and designing a representative training set, as detailed in the following section.
## 3 A learning approach for preconditioning parametrized systems

We shall now discuss the problem of learning preconditioners for parametrized systems. Although this task is inherently related to finite-dimensional spaces, we find it convenient to explore and present the idea within an infinite-dimensional context. This shift to a function space setting, as presented in this section, provides a deeper insight into the mathematical structure of the problem. For these reasons, we start this section by fixing some notation. Given two Banach spaces $(V_1, \|\cdot\|_{V_1})$ and $(V_2, \|\cdot\|_{V_2})$ we denote $\mathcal{L}(V_1, V_2)$ the space of bounded linear operators from $V_1$ to $V_2$, equipped with the operator norm

$$ \|A\|_{\mathcal{L}(V_1, V_2)} := \sup_{v \in V_1 \setminus \{0\}} \frac{\|Av\|_{V_2}}{\|v\|_{V_1}}. $$

It is also useful to introduce the General Linear Group

$$ \mathcal{GL}(V_1, V_2) := \{A \in \mathcal{L}(V_1, V_2) : \exists B \in \mathcal{L}(V_2, V_1), BA = \text{Id}_{V_1} \text{ and } AB = \text{Id}_{V_2}\}, $$

which consists of all linear operators from $V_1 \to V_2$ that are both invertible and continuous. Here, $\text{Id}_{V_i}$ denotes the identity operator of $V_i$ on itself. Given a Banach space $(V, \|\cdot\|_V)$, we also denote by $V'$ its topological dual, that is, $V' := \mathcal{L}(V, \mathbb{R})$. Our purpose is to address the problem of operator preconditioning for parametrized linear problems in arbitrary (complete) normed spaces, i.e. to find a continuous preconditioner operator $\mathcal{P} = \mathcal{P}(v, \mu)$ for a parametrized problem of the form

$$ \text{find } u^{\mu} \in V : \quad A^{\mu}u^{\mu} = b^{\mu}, \quad \forall \mu \in \mathcal{P}, \qquad (12) $$

given a Banach state space $(V, \|\cdot\|_V)$ and a compact metric space $\mathcal{P}$ serving as parameter space. Here, $A^{\mu} \in \mathcal{GL}(V, V')$ and $b^{\mu} \in V'$ are parameter-dependent operators and problem data. In general, we think of the latter as a linear system arising from the discretization of a given PDE, such as (5) or (6). Everything, in fact, can be traced back to our discussion in Section 2: we refer to Remark 3 for additional insight on the matter. Furthermore, we anticipate that the upcoming Section 4 will also contain further clarifications on how these concepts translate when put into practice.

We organize the remainder of this Section into distinct parts. Initially, we explore how operator preconditioning is applied to (12), emphasizing the imperative for a nonlinear preconditioner (Section 3.1). Subsequently, we outline the structure of the learning algorithm, presenting a suitable unsupervised training strategy (Section 3.2).

**Remark 3.** The advantage of working within an abstract environment is that the idea can be easily transferred from the continuous formulation to the discrete one and vice versa. For instance, following our notation in Section 2.2, if we resort to the discrete setting, then $A^{\mu} \cong C_h^{\mu}$ is a matrix, while $u^{\mu} \cong v_h^{\mu}$ and $b^{\mu} \cong \tilde{b}_h^{\mu}$ are finite-dimensional vectors; cf. Eq. (9). Conversely, if we lift the idea to the continuous level, then $A^{\mu}$ is the linear operator acting as $u \mapsto a_{\Omega}(u, \cdot) + b_{\Lambda}^{\mu}(T_{\Lambda}u, \cdot)$, while $u^{\mu}$ and $b^{\mu}$ are functions defined over $\Omega$: cf. Section 2.1.

### 3.1 The need for a nonlinear and matrix-free preconditioner

Classical approaches to operator preconditioning focus on finding a linear preconditioner for a given operator of interest $A \in \mathcal{L}(V', V)$. In principle, one could apply this idea to parametric problems by addressing each parametric scenario separately. Indeed, for a generic, but fixed, parameter $\bar{\mu} \in \mathcal{P}$ one could consider the corresponding operator $A^{\bar{\mu}} \equiv A$, with $A: V \to V'$, and then look for a suitable (right) preconditioner that is, any linear operator $P \in \mathcal{L}(V', V)$ for which $\kappa(AP) \ll \kappa(A)$, being $\kappa$ the operator condition number $\kappa(B) := \|B\|_{\mathcal{L}(V_1, V_2)} \|B^{-1}\|_{\mathcal{L}(V_1, V_2)}$ for $B \in \mathcal{GL}(V_1, V_2)$. Then, once a suitable preconditioner is available, the original problem $Au^{\bar{\mu}} = b^{\bar{\mu}}$ can be replaced with a new one, $APz^{\bar{\mu}} = b^{\bar{\mu}}$, which is easier to solve, so that $u^{\bar{\mu}} = Pz^{\bar{\mu}}$. In general, we note that up to isometries, the "optimal" preconditioner is $P = A^{-1}$, regardless of the right hand side $b^{\bar{\mu}} \in V'$. However, this approach becomes computationally unfeasible if the procedure is repeated multiple times for varying $\mu \in \mathcal{P}$. A first alternative could be to look for a surrogate $\tilde{\mathcal{P}}: \mathcal{P} \to \mathcal{L}(V', V)$ that, for each $\mu \in \mathcal{P}$, returns a suitable preconditioner for $A^{\mu}u^{\mu} = b^{\mu}$. However, we do not find this solution to be very practical as that would entail casting the learning problem within an
extremely high-dimensional space. Notice, in fact, that in the discrete setting one has dim $\mathcal{L}(V', V) = (N_h^\Omega)^2$. Furthermore, even if we restrict our attention to sparse preconditioners, this is still very challenging as the sparsity pattern would likely depend on the structure of the underlying 1D problem, i.e. on $\mu \in \mathcal{P}$. A more practical approach can be to rely on a matrix-free formulation. There, one is not directly interested in the preconditioner by itself, but rather in the way it acts when applied to an input vector. Mathematically speaking, this change of perspective corresponds to looking for a map, potentially nonlinear, of the form $\mathcal{P}: V' \times \mathcal{P} \to V$. This idea is further motivated by the following Lemma.

**Lemma 1.** Let $\mathcal{P} \ni \mu \mapsto A^\mu \in \mathcal{L}(V, V')$ be continuous. For any $\mu \in \mathcal{P}$ and any $v \in V'$ let $x^{\mu,v}$ be the solution of $A^\mu x^{\mu,v} = v$. Then, there exists a continuous (nonlinear) operator $\mathcal{P}: V' \times \mathcal{P} \to V$ such that

$$ \mathcal{P}(v, \mu) = x^{\mu,v} \quad \forall (v, \mu) \in V' \times \mathcal{P}. $$

In particular, $A^\mu \mathcal{P}(v, \mu) = v$ for all $(v, \mu) \in V' \times \mathcal{P}$.

*Proof.* Let $\mathcal{A}: \mathcal{P} \to \mathcal{GL}(V, V')$ be the map $\mathcal{A}(\mu) = A^\mu$. Let $\mathcal{I}: \mathcal{GL}(V, V') \to \mathcal{GL}(V', V)$ be the inversion map, $B \mapsto B^{-1}$. It is well known that $\mathcal{I}$ is continuous — for a rigorous proof, we refer the interested reader to the Appendix, Lemma A.1. Let $\mathcal{P}: V' \times \mathcal{P} \to V$ be defined as $\mathcal{P}(v, \mu) := [(\mathcal{I} \circ \mathcal{A})(\mu)](v)$. Then, $\mathcal{P}$ is continuous and $\mathcal{P}(v, \mu) = (A^\mu)^{-1}v = x^{\mu,v}$, as claimed. $\square$

As seen in Lemma 1, in order to adapt to different parametric scenarios, the preconditioner $\mathcal{P}$ must be nonlinear in its arguments (see also Remark 4). In light of this, hereon, when speaking of a nonlinear preconditioner for a parametric problem such as (12), we intend any continuous map $\mathcal{P}: V' \times \mathcal{P} \to V$ for which $A^\mu \mathcal{P}(v, \mu) \approx v$, or equivalently, $(\mathrm{Id} - A^\mu \mathcal{P}(\cdot, \mu))v \approx 0$. Ideally, this should hold for all $\mu \in \mathcal{P}$ and all right-hand sides of interest $v \in \mathcal{K}^\mu \subset V'$, with the latter subset potentially depending on the parameter instance $\mu$. For example, the simplest choice could be $\mathcal{K}^\mu = \{b^\mu\}$, so that $\mathcal{P}(b^\mu, \mu) \approx u^\mu$. Setting $\mathcal{K}^\mu = V'$ is also possible but not strictly necessary: we recall, in fact, that our ultimate goal is to solve (12). We shall return to this in Section 3.2, as we address the details of our learning strategy, and in Section 4.1, where we discuss the actual implementation of the approach and the definition of $\mathcal{K}^\mu$.

**Remark 4.** In principle, one could restrict oneself to continuous maps $\mathcal{P}: V' \times \mathcal{P} \to V$ that are linear in their first argument, $\mathcal{P} = \mathcal{P}(v, \mu)$, as it happens for the "ideal" preconditioner in Lemma 1. However, allowing for more general maps gives additional freedom in designing the preconditioner, potentially simplifying the learning process. For example, assume that $\mathcal{P}$ is compact and that the map $\mathcal{F}: \mu \to b^\mu$ is both continuous and injective. Then, $\mathcal{F}^{-1}: \mathcal{F}(\mathcal{P}) \to \mathcal{P}$ exists, is continuous and can be extended to $V' \supset \mathcal{F}(\mathcal{P})$ using, e.g., Dugundji's extension of Tieze's Theorem [11]. With this setup, consider the simple scenario $\mathcal{K}^\mu = \{b^\mu\}$ and let $\mathcal{P}$ be the map in Lemma 1. Define $\tilde{\mathcal{P}}: V' \times \mathcal{P} \to V$ as $\tilde{\mathcal{P}}(v, \mu) := \mathcal{P}(v, \mathcal{F}^{-1}(v))$. Then, $\tilde{\mathcal{P}} \neq \mathcal{P}$ is nonlinear in $v$, does not depend explicitly on $\mu$ (only implicitly), but still $A^\mu \tilde{\mathcal{P}}(v, \mu) = v$ for all $\mu \in \mathcal{P}$ and all $v \in \mathcal{K}^\mu \subset V'$.

## 3.2 Learning strategy

Following our previous discussion, our purpose here is to learn a (nonlinear) preconditioner for a parametrized problem given as in Eq. (12). In practice, this involves fixing a suitable hypothesis space $\mathcal{H}$ and finding the optimal preconditioner $\mathcal{P} \in \mathcal{H}$ by minimizing a suitable functional known as loss function. In light of Lemma 1, we let $\mathcal{H} \subset \mathcal{C}(V' \times \mathcal{P}, V) = \{\mathcal{P}: V' \times \mathcal{P} \to V \text{ continuous}\}$ be a subset of the space of continuous functions from $V' \times \mathcal{P}$ to $V$. For instance, if we resort to deep learning techniques, as we shall do in the remaining of the paper, the latter can be thought of as the space of all possible neural network models associated to a given architecture of choice. In this case, we'll refer to $\mathcal{P}$ as to a neural preconditioner. Concerning the loss function, instead, different approaches become available. In principle, one could define the nonlinear preconditioner $\mathcal{P}$ by requiring either $\mathcal{P}(v, \mu) \approx x^{\mu,v}$ or $(\mathrm{Id} - A^\mu \mathcal{P}(\cdot, \mu))v \approx 0$. In fact, we notice that if the original problem admits a unique solution for each $\mu \in \mathcal{P}$, then every $\mathcal{P}$ that fulfills $\mathcal{P}(b^\mu, \mu) = u^\mu$ yields $A^\mu \mathcal{P}(b^\mu, \mu) = b^\mu$ and vice versa. In this sense, the two approaches may seem equivalent. However, our hypothesis space will typically not contain the global minimizer in Lemma 1, reason for which minimizing the error $\|\mathcal{P}(v, \mu) - x^{\mu,v}\|$ or the residual $\|(\mathrm{Id} - A^\mu \mathcal{P}(\cdot, \mu))v\|$ can provide completely different results. Furthermore, the two approaches differ significantly in their implementation: the first one is supervised, as it
requires actually solving (12) in order to compute $x^{\mu,v}$ for multiple $\mu \in \mathcal{P}$; in contrast, the second approach is fully unsupervised.

This work focuses exclusively on the second approach. This is because preliminary investigations by the authors found the supervised approach to be extremely inefficient, with high costs during the training phase, significant challenges in minimizing the loss function, and unsatisfactory results. In contrast, the unsupervised approach is more practical, as it does not require sampling from the solution manifold $\mathcal{S} = \{u^{\mu}\}_{\mu \in \mathcal{P}}$. All of this considered, a first definition of the learning problem can be

$$ \mathcal{P}_* = \arg \min_{\mathcal{P} \in \mathcal{H}} \int_{\mathcal{P}} \|(\mathrm{Id} - A^{\mu} \mathcal{P}(\cdot, \mu))b^{\mu}\|_{V'} \varrho(d\mu), \quad (13) $$

for a given probability measure $\varrho$ defined over the parameter space $\mathcal{P}$. The downside of this approach is that for each $\mu \in \mathcal{P}$, we focus solely on how well $\mathcal{P}$ operates on $b^{\mu}$. As anticipated in Section 3.1, a better approach can be to introduce a subset $\mathcal{K}^{\mu} \subset V'$ of possible right-hand-sides (depending in general on $\mu$), and reformulate (13) as

$$ \mathcal{P}_* = \arg \min_{\mathcal{P} \in \mathcal{H}} \int_{\mathcal{P}} \int_{\mathcal{K}^{\mu}} \|(\mathrm{Id} - A^{\mu} \mathcal{P}(\cdot, \mu))v\|_{V'} \varrho_{\mu}(dv) \varrho(d\mu), \quad (14) $$

where $\varrho_{\mu}$ is a suitable probability measure supported over $\mathcal{K}^{\mu} \subset V'$. Clearly, (13) is a special case of (14), since it can be obtained letting $\mathcal{K}^{\mu} = \{b^{\mu}\}$ be a singleton and $\varrho_{\mu} = \delta_{b^{\mu}}$ be a Dirac delta measure. The question of defining $\mathcal{K}^{\mu}$ and the corresponding probability measure is generally intriguing, yet quite complex. We will defer this discussion to Section 4.1. We anticipate that $\mathcal{K}^{\mu}$ should contain all the components that are detrimental to the well-posedness of $A^{\mu}$ and, in general, to the convergence rate of the algebraic solver considered.

From a practical point of view, the ideal minimization problem in Eq. (14) can be addressed through empirical risk minimization [20, 6]. Specifically, let $\mu_1, \dots, \mu_{N_P}$ be an independent random sample identically distributed (i.i.d.) with $\mu_j \sim \varrho$. For each $j = 1, \dots, N_P$ let $\mathcal{K}^{\mu_j} := \{v_{1,j}, \dots, v_{N_{\mathcal{K}^{\mu_j}}, j}\} \subset \mathcal{K}^{\mu_j}$, be an i.i.d. random sample with $v_{i,j} \sim \varrho_{\mu_j}$. Consider the empirical measures

$$ \varrho_{N_P} = \frac{1}{N_P} \sum_{j=1}^{N_P} \delta_{\mu_j}, \quad \varrho_{N_{\mathcal{K}^{\mu_j}}} = \frac{1}{N_{\mathcal{K}^{\mu_j}}} \sum_{v \in \mathcal{K}_{\mu_j}} \delta_v, \quad \varrho_{N^{\mathcal{K}^{\mu}}} = \sum_{j=1}^{N_P} \mathbf{1}_{\{\mu_j\}}(\mu) \varrho_{N_{\mathcal{K}^{\mu_j}}}, $$

where $\delta_x$ is the Dirac measure centered at $x$. Following classical procedures, we define the empirical risk as

$$ \begin{aligned} \mathfrak{R}(\mathcal{P}) &:= \int_{\mathcal{P}} \int_{\mathcal{K}^{\mu}} \|(\mathrm{Id} - A^{\mu} \mathcal{P}(\cdot, \mu))v\|_{V'} \varrho_{N_{\mathcal{K}^{\mu}}}(dv) \varrho_{N_P}(d\mu) = \\ &= \frac{1}{N_P} \sum_{j=1}^{N} \frac{1}{N_{\mathcal{K}^{\mu_j}}} \sum_{v \in \mathcal{K}_{\mu_j}} \|(\mathrm{Id} - A^{\mu_j} \mathcal{P}(\cdot, \mu_j))v\|_{V'}. \end{aligned} \quad (15) $$

Then, the training phase consists in solving the following minimization problem

$$ \tilde{\mathcal{P}}_* := \arg \min_{\mathcal{P} \in \mathcal{H}} \mathfrak{R}(\mathcal{P}). \quad (16) $$

Note that while this step might be computationally demanding, it has to be performed only once. After training, the preconditioner can be readily applied to any problem in the parametric class, without further optimizations or assembling stages. In particular, if one designs the hypothesis space to consist of neural network models, this can provide significant speed ups, as the computational cost associated with online evaluations of $\tilde{\mathcal{P}}_*$ becomes nearly negligible (on the order of milliseconds).

# 4 Realization of the learning approach

In this section, we return to the finite-dimensional setting to make the preconditioner learning strategy specific to the discretization of the parametrized mixed-dimensional problem (9)-(10). As we address here the 3D problem arising from the elimination of the 1D equations from the coupled system, with a slight abuse
of notation, we write $N_h$ instead of $N_h^\Omega$. Once a suitable preconditioner $\mathcal{P}: \mathbb{R}^{N_h} \times \mathbb{R}^{N_h} \to \mathbb{R}^{N_h}$ has been identified (we recall that the finite-dimensional parameter space is $\mathcal{P}_h \subset \mathbb{R}^{N_h}$), it can be used to accelerate the solution of the discrete problem. In general, this can be achieved in multiple ways. Here, we propose the integration of the nonlinear preconditioner within a GMRES solver,motivated by the fact that we cannot guarantee the symmetry of the neural preconditioner. In addition, considering the nonlinear nature of the preconditioner $\mathcal{P}$, we shall adopt the Flexible GMRES (FGMRES) algorithm [38], as the latter can easily accommodate variations of the preconditioning operator with the input vector.

## 4.1 The training set $\mathcal{K}^\mu$

The selection of the training set $\mathcal{K}^\mu$ is a crucial factor in the success of the operator learning approach for preconditioning. In fact, the training set must match the underlying behavior of the operator to ensure that the learned preconditioner is effective in different instances of the parameter space $\mathcal{P}$. To start, we notice that inside the FGMRES algorithm, the preconditioner will act only on vectors from the unit sphere $\mathbb{S}^{N_h-1} := \{\mathbf{v} \in \mathbb{R}^{N_h} : \|\mathbf{v}\| = 1\}$. Thus, we let $\mathcal{K}^\mu \subset \mathbb{S}^{N_h-1}$. Next, based on the findings presented in [8], we propose that the training set should be guided by the kernel structure of the operator. In fact, the coupling between the 3D and 1D domains introduces significant challenges in the construction of efficient solvers, particularly because the kernel of the operator $M_{h,00}^\mu$ contains high-frequency components that must be taken into account by the computational solvers. In the context of the algebraic multigrid (AMG) approach addressed in [8], it is shown that an essential property for the robustness of the AMG solver is the following kernel decomposition condition:

$$ \ker(M_{h,00}^\mu) = \ker(M_{h,00}^\mu) \cap V_c + \sum_{j=1}^{J} \ker(M_{h,00}^\mu) \cap V_j, $$

where in the case of AMG algorithms $V_c$ refers to the coarse subspace that captures the low-frequency components of the solution and $V_j$ are the fine subspaces that represent localized corrections in specific regions or aggregates. In the context of a multiresolution method such as AMG, this condition suggests that the numerical solver should operate on all components of the kernel of the coupling operator. The richer the kernel subspace, possibly including low- and high-frequency components, the harder it is to develop good solvers. This issue is closely related to the concept of spectral bias in neural network training, which states that lower-frequency components are learned more easily and earlier during training, while higher-frequency components are more challenging to learn [35, 42]. For the approximation of mixed-dimensional problems, this spectral bias has direct implications. Specifically, the high-frequency components in the kernel of the coupling operator $M_{h,00}^\mu$ may not be adequately captured if the training process is not designed to address this inherent difficulty. To mitigate this challenge, we include in $\mathcal{K}^\mu$ samples that span both the low-frequency and high-frequency components of the kernel of $C_h^\mu$, defined in (10). In the context of operator learning, this insight is reflected in the definition of the measure $\varrho_\mu(dv) = \varrho_{N_{\mathcal{K}^\mu}}$. This considered, we propose the following training set structure:

$$ \mathcal{K}^\mu = \bigcup_{j=1}^{N_P} \mathcal{K}^{\mu_j}, \quad \mathcal{K}^{\mu_j} := \left\{ \frac{b_h^{\mu_j}}{\|b_h^{\mu_j}\|} \right\} \cup \mathcal{D}^{\mu_j}, \quad b_h^{\mu_j} \in \mathcal{B}; \qquad (17) $$

where $\mathcal{D}^{\mu_j}$ represents a suitable data-augmentation set, given the parameter $\mu_j$ (see Fig. 2). In this way, the operator $\mathcal{P}(\cdot, \mu)$ can be trained to act as a preconditioner on a suitable subset of the unit sphere with a specified frequency content. In what follows, we examine two distinct strategies for data augmentation, while the next section will present numerical experiments that support this methodology and demonstrate that incorporating both right-hand side vectors and kernel components into $\mathcal{K}^\mu$ improves the convergence rates in FGMRES with the application of the learned preconditioner.

### 4.1.1 Augmented Training Set by Krylov Subspaces

Our initial data-augmentation strategy involves incorporating high-frequency elements derived from the construction of Krylov subspaces related to the matrix. As a starting point, the Krylov subspace associated with a general operator $O \in \mathcal{L}(X, X)$ and a vector $x \in X$ is defined as:

$$ \mathcal{K}_q(O, x) := \operatorname{span}(x, Ox, O^2x, \dots, O^{q-1}x) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Pictorial summary for the construction of the augmented training set $\mathcal{K}^{\mu}$. The data augmentation set $\mathcal{D}^{\mu_j}$ contains vectors in the unit sphere with the desired frequency content. Here, $\mathcal{B}/\|\mathcal{B}\| := \{b/\|b\| : b \in \mathcal{B}\} \subset \mathbb{S}^{N_h-1}$ is the set of normalized right-hand side vectors.</figcaption></figure>

where $q$ is the order of the Krylov subspace. This subspace contains progressively higher-order components of $x$ under the action of the operator $O$, capturing information about the behavior of the system across multiple scales. By including vectors from these Krylov subspaces in our training set, we effectively augment the representation of both low- and high-frequency components in the solution space.

Specifically, for a generic parameter $\mu_j$, let $O = C_h^{\mu_j}$ be the discrete operator associated with the mixed-dimensional PDEs described in previous sections. We construct the data-augmentation set $\mathcal{D}^{\mu_j}$ as

$$ \mathcal{D}^{\mu_j} := \{ q_{p'}^{\mu_j} \cdots q_{(p+p')}^{\mu_j} \}, $$

where $0 \le p'$ is a suitable shift index and $q_i^{\mu_j}$, $i = 0, \dots, p'+p$ is a suitable orthonormal basis of $\mathcal{K}_{(p+p')}(C_h^{\mu_j}, b_h^{\mu_j})$ with $q_0^{\mu_j} = b_h^{\mu_j}/\|b_h^{\mu_j}\|$. These additional basis functions aim to capture a wider spectrum of the operator's behavior, thus enhancing the preconditioner's ability to address different scales and frequencies present in the solution. In this context, it is convenient to keep $p$ and $p'$ small and of the same magnitude (i.e. $p \simeq p'$). This lies in the balance between representing the essential high-frequency characteristics of the operator $C_h^{\mu_j}$ and maintaining computational efficiency. The choice of small $p, p'$ allows us to efficiently construct the augmented training set without significantly increasing the computational burden during the training phase. As will be shown later (6), even with a small value of $p+p'$, the Krylov subspace captures important high-order effects that improve the robustness of the preconditioner.

### 4.1.2 Augmented Training Set by Random Vectors

A second approach involves adding to the physics-based training subset $\mathcal{B}$ a new set of completely unrelated random functions. Specifically, we propose augmenting the training set by generating random unity vectors $\mathbf{r}_h \in \mathbb{S}^{N_h-1} \subset \mathbb{R}^{N_h}$, uniformly sampled from the unit hypersphere. In practice, each $\mathbf{r}_h$ can be generated by drawing and normalizing a random vector $\mathbf{v} \sim \mathcal{N}(0, \mathbf{I}_{N_h})$ with standard multivariate normal distribution, that is, $\mathbf{r}_h = \mathbf{v}/\|\mathbf{v}\|$. Such vectors are designed to introduce high-frequency content that is independent of the physics of the problem, thereby complementing the functions based on physics $b_h^{\mu} \in \mathcal{B}$. In fact, given a finite element discretization based on a FEM space $V_h(\Omega)$, where $\Omega$ represents the three-dimensional computational domain and $h$ denotes the characteristic element size, the vectors $\mathbf{r}_h$ can be interpreted as the highest frequency modes that can be resolved by the mesh. The frequency associated with these random vectors is proportional to $|\Omega|/h$, which corresponds to the smallest possible wavelength that can be captured in the discretized domain. Mathematically, we can define the data augmentation set as $\mathcal{D}^{\mu_j} := \{r_{h,1}, r_{h,2}, \dots\}$. Note that these vectors are independent of $\mu$ but a different instance of $\mathcal{D}$ is defined for each $\mu_j$. We expect this strategy to improve robustness with respect to high-frequency vectors that are inside $\ker(M_{h,00}^{\mu_j})$.
## 4.2 Hypothesis space for $\mathcal{P}$

At this stage, it remains to define the hypothesis space $\mathcal{H}$, over which the empirical risk minimization will be performed (16). Following a deep-learning approach, the hypothesis space is specified once a particular neural architecture $\mathcal{N}_H$ is selected. This architecture is characterized by a composition of linear and nonlinear functions, which in turn are defined by a set of hyper-parameters $H$. The resulting hypothesis space $\mathcal{H}$ can thus be expressed as $\mathcal{H} = \{\mathcal{N}_H(\cdot; \theta)\}_{\theta \in \mathbb{R}^t}$ where $\theta$ is the vector of trainable parameters associated with $\mathcal{N}_H$, over which the optimization is performed. Drawing inspiration from cutting-edge deep learning techniques in image processing, as proposed in [2], which employ varying resolutions to effectively manage the diverse scales of images through pooling and upscaling operations, the neural architecture chosen to represent the preconditioner will be a U-net, see for example [36] for a landmark paper and [41] for recent developments. This architecture is highly capable of addressing multi-scale issues as it effectively captures both global and local features across varying resolutions. This makes it a strong candidate function for representing solutions of mixed-dimensional problems, such as the ones represented in Figure 1. According to this choice, we define $\mathcal{N}_H(\cdot; \theta) \equiv \mathcal{U}_L(\cdot; \theta)$, where $\mathcal{U}_L$ represents a U-net architecture with $L$ levels. For the sake of simplicity, we will refer only to the hyper-parameter $L$ to indicate the depth of the U-net, although additional parameters such as bottleneck size, number of channels, and convolutional kernel sizes are also involved. Consequently, the desired nonlinear preconditioner operator can be represented as $\mathcal{P} \equiv \mathcal{U}_L(\cdot; \theta^*)$, where $\theta^*$ denotes the set of optimized parameters. The U-net architecture is structured as a combination of two components, an encoder function $\Phi$ and a decoder function $\Psi$:

$$ \Phi := \Phi_0 \circ \Phi_1 \circ \dots \circ \Phi_{L-2} \circ \Phi_{L-1} : \mathbb{R}^{c_{in} \times n} \to \mathbb{R}^{c_{b} \times m}; $$

$$ \Psi := \Psi_{L-1} \circ \Psi_{L-2} \circ \dots \circ \Psi_1 \circ \Psi_0 : \mathbb{R}^{c_b \times m} \to \mathbb{R}^{c_{out} \times n}. $$

Then, given an input tensor $\mathbf{X} = [\mathbf{X}_1|\mathbf{X}_2|\dots|\mathbf{X}_{c_{in}}]$, where $\mathbf{X}_i \in \mathbb{R}^{n_1 \times n_2 \times \dots \times n_d}$, and a prescribed number of output channels $c_{out}$, a U-net with $j$ levels can be defined recursively as follows:

$$ \mathcal{U}_j(\mathbf{X}) := \Psi_j([\mathcal{U}_{j-1} \circ \Phi_j|\Phi_j]) (\mathbf{X}), \quad \mathcal{U}_0(\mathbf{X}) := \Psi_0 \circ \Phi_0(\mathbf{X}), $$

where the notation $[a|b]$ represents a channel stacking operation that concatenates tensors $a$ and $b$ along the channel dimension. Specifically, if $a \in \mathbb{R}^{c_a \times m}$ and $b \in \mathbb{R}^{c_b \times m}$, then we have $[a|b] \in \mathbb{R}^{(c_a+c_b) \times m}$. Operation $[\mathcal{U}_{j-1} \circ \Phi_j(\mathbf{X})|\Phi_j(\mathbf{X})]$ is called the $j$-th skip connection, which transmits information from the $j$-th encoding layer to the corresponding $j$-th decoding layer. Skip connections are critical in retaining and merging high-resolution features, ensuring that the decoder can effectively recover fine-scale information. In the literature, the term $c$ refers to the number of channels at each layer. Specifically, $c_b$ is the number of channels in the bottleneck layer of the U-net, while $c_{in}$ and $c_{out}$ denote the number of input and output channels, respectively. The convolutional layer operates locally by sliding a kernel across the input tensor, extracting features from small spatial neighborhoods. Given a $d$-dimensional input tensor $\mathbf{X} \in \mathbb{R}^{c_{in} \times n_1 \times n_2 \times \dots \times n_d}$, the convolution operation is defined as:

$$ \mathbf{Y}_i = \sum_{j=1}^{c_{in}} \mathbf{k}_{i,j} * \mathbf{X}_j \text{ for } i = 1, 2, \dots, c_{out}, \text{ with } (\mathbf{k}_{i,j} * \mathbf{X}_j)(\mathbf{p}) = \sum_{\mathbf{q} \in \mathbb{Z}^d} \mathbf{k}_{i,j}(\mathbf{q}) \mathbf{X}_j(\mathbf{p} - \mathbf{q}), $$

$\mathbf{k}_{i,j} \in \mathbb{R}^{s_1 \times s_2 \times \dots \times s_d}$ being the kernel for the $(i,j)$-th input-output channel pair, $\mathbf{p}$ indexing the output tensor, and $\mathbf{q}$ spanning the kernel positions. This operation enables convolutional layers to act as local feature extractors, capturing patterns from the input tensor. Input data, often provided as ordered vectors, must be reshaped into tensor format to preserve spatial correlations. Once trained, the kernel can operate on input tensors of arbitrary dimensions, provided the structure is compatible, such as tensor-product meshes. As a result, the U-Net architecture can handle inputs of any size, which may increase its application beyond the dimensions of the training data. However, reliance on structured grids restricts the application of convolutional layers to simple geometries. Future work will address this limitation by integrating mesh-informed neural networks, as proposed in [13], to generalize U-net architectures while retaining their feature extraction capabilities.
<figure><img src="image_4.png" /><figcaption>Figure 3: Schematic representation of the U-Net architecture $\mathcal{U}_3$. Tensor data are represented by blocks, where the number of channels corresponds to the depth of the blocks. The tensor shape, $c \times n_1 \times n_2 \times n_3$, is reported on the top of the blocks. Solid line arrows represent layer action, while dashed arrows represent channel-staking skip connection.</figcaption></figure>

<table><caption>Table 1: Summary of the $\mathcal{U}_3$ architecture. Input data flows from the top row to the bottom row; see also the scheme in Fig.3. The total number of trainable parameters is $\theta \approx 10^6$.</caption><thead><tr><th>Layer type</th><th>input size</th><th>output size</th><th>input channels</th><th>output channels</th><th># parameters θ</th></tr></thead><tbody><tr><td>Conv.</td><td>21<sup>3</sup></td><td>21<sup>3</sup></td><td>2</td><td>32</td><td>13,840</td></tr><tr><td>Conv.+Max P.</td><td>21<sup>3</sup></td><td>10<sup>3</sup></td><td>32</td><td>64</td><td>55,296</td></tr><tr><td>Conv.+Max P.</td><td>10<sup>3</sup></td><td>5<sup>3</sup></td><td>64</td><td>128</td><td>221,184</td></tr><tr><td>Conv.</td><td>5<sup>3</sup></td><td>5<sup>3</sup></td><td>128</td><td>128</td><td>442,368</td></tr><tr><td>Trasp. conv. + Skip</td><td>5<sup>3</sup></td><td>10<sup>3</sup></td><td>128</td><td>64</td><td>286,784</td></tr><tr><td>Trasp. conv. + Skip</td><td>10<sup>3</sup></td><td>21<sup>3</sup></td><td>64</td><td>32</td><td>71,912</td></tr><tr><td>Conv.</td><td>21<sup>3</sup></td><td>21<sup>3</sup></td><td>32</td><td>1</td><td>545</td></tr></tbody></table>

# 5 Numerical solution of mixed-dimensional PDEs using the neural preconditioner

We analyze how the neural preconditioner accelerates the convergence of iterative solvers, in the case of mixed-dimensional PDEs. We discuss the computational setup and analyze the preconditioner's impact on FGMRES convergence rate, mesh size scalability, and computational time. In doing so, we shall take the opportunity to assess the importance of the several design choices entailed by the implementation of the neural preconditioner, such as, e.g., the data-augmentation strategy and the parametrization of the 1D geometry. For the sake of this preliminary exploration, in this Section we shall restrict ourselves to numerical solution of problem (10), that is, to the 3D block of the mixed-dimensional system, thus ignoring the 3D-1D coupling. The latter will be addressed right after, in Section 6.

We anticipate that the main message emerging from this numerical exploration is that the neural preconditioner proves competitive when compared to other state-of-the-art algorithms —such as AMG or ILU— thanks to a reduction in the computational cost of each preconditioning step. However, the iteration count and the scalability of the neural preconditioner with respect to the mesh size are not optimal. In fact, as shown later in Table 9, a trade-off between computational cost and algorithmic efficiency appears. However, our findings indicate that, in general, the neural preconditioner offers a competitive approach to solving mixed-dimensional PDEs, which motivates its use in Section 6 when addressing the coupled problem.

## 5.1 General setup of the numerical tests

The domain $\Omega$ is the unit cube and the problem is discretized using a structured mesh with 9261 nodes (21³ grid points). The physical constants are set as $k_\Omega = 10^{-3}$, $\sigma_\Omega = 10^{-3}$ and $\epsilon = 10^{-3}$. In addition, we performed tests using $k_\Omega = 10^{-2}$, $\sigma_\Omega = 10^{-2}$, confirming the adaptability of the method. For the training
<figure><img src="image_5.png" /><figcaption>Figure 4: Examples of the graphs $\Lambda$ considered in the numerical tests, with increasing geometrical complexity.</figcaption></figure>

set $\mathcal{K}^{\mu}$, different 1D graph geometries are considered (see Fig. 4), ranging from $O(1)$ to $O(10^2)$ number of branches.

We considered a three-level U-net architecture, $\mathcal{U}_3$, with two input channels: one relative to the vector in the training set $\mathcal{K}^{\mu}$ and one carrying the parameter information $\mu$, that is, taking the discrete distance function $\mathbf{d}_{h}^{\mu} \in \mathbb{R}^{N_h}$ associated with the graph $\Lambda$. More technical details can be found in Table 1 and Fig. 3. The actual training set $\tilde{\mathcal{B}} \approx \mathcal{B}$ is obtained by solving approximately the system in (11), with a relative tolerance of $10^{-4}$, thus ensuring computational efficiency in training set generation. The data-augmentations sets $\mathcal{D}^{\mu j}$ vary depending on the scenario analyzed, resulting in three different U-Net preconditioners:

$$ 
\begin{align*}
\mathcal{D}^{\mu j} &= \emptyset &\Rightarrow &\mathcal{U}_3^{\emptyset}(\cdot; \theta), \\
\mathcal{D}^{\mu j} &= \{q_5^{\mu j}, \dots, q_8^{\mu j}\} &\Rightarrow &\mathcal{U}_3^{\text{kry}}(\cdot; \theta), \\
\mathcal{D}^{\mu j} &= \{r_{h,1}, \dots, r_{h,4}\} &\Rightarrow &\mathcal{U}_3^{\text{hf}}(\cdot; \theta).
\end{align*}
$$

Here, $\{q_5^{\mu j}, \dots, q_8^{\mu j}\}$ are orthonormal bases of the Krylov subspace $\mathcal{K}_8(C_h^{\mu}, b_h^{\mu})$, and $\{r_{h,1}, \dots, r_{h,4}\}$ are high-frequency vectors sampled from the unit sphere $\mathbb{S}^{N_h-1}$. The networks $\mathcal{U}_3$ are trained through the minimization of the empirical risk:

$$ \mathcal{U}_3^* = \mathcal{U}_3(\theta^*) \text{ with } \theta^* = \arg\min_{\theta} = \frac{1}{N_P} \sum_{j=1}^{N_P} \frac{1}{N_{\mathcal{K}_{\mu j}}} \sum_{\mathbf{v} \in \mathcal{K}_{\mu j}} \|\mathbf{v} - \alpha^{-1} C_h^{\mu j} \mathcal{U}_3(\mathbf{v}, \mathbf{d}_h^{\mu j}; \theta)\|^2, $$

where $\alpha := \left( \frac{1}{N_P} \sum_{j=1}^{N_P} \frac{1}{\sqrt{N_h}} \| \mathrm{diag}(C_h^{\mu j}) \|_2 \right)$ is a fixed normalizing constant representing the mean matrices diagonal root mean square, which helps stabilize the training process. For comparison purposes, an equal number of graph samples and training epochs is used for $\mathcal{U}_3^{\emptyset}$, $\mathcal{U}_3^{\text{kry}}$, and $\mathcal{U}_3^{\text{hf}}$. More data on the training process are reported in Table 2.

<table><caption>Table 2: Details on the training dataset and parameters. Due to data augmentation strategies, the total number of training samples is five times the one of the 1D graphs.</caption><thead><tr><th>Training</th><th>1D graphs</th><th>Tot. samples</th><th>Validation</th><th>1D graphs</th><th>Tot. samples</th></tr></thead><tbody><tr><td>Samples</td><td>480</td><td>2400</td><td>Samples</td><td>120</td><td>600</td></tr><tr><td>Batch Size</td><td></td><td>5</td><td>Epochs</td><td></td><td>250</td></tr><tr><td>Learn. rate</td><td></td><td>10<sup>-3</sup></td><td>Learn. strategy</td><td></td><td>decaying</td></tr></tbody></table>

Figure 5 displays the convergence history of training and validation for the U-Net architecture when the physical constants $k_{\Omega} = 10^{-3}$, $\sigma_{\Omega} = 10^{-3}$ are considered. The graphs illustrate that the training and validation error for $\mathcal{U}_3^{\emptyset}$, $\mathcal{U}_3^{\text{kry}}$ and $\mathcal{U}_3^{\text{hf}}$ possesses similar trends but quite different values; nevertheless, as will be clarified later on, the training and validation error alone are not discriminant for the performance of the U-Net preconditioner. What is fundamental is the data set over which the given error is obtained; e.g., even if $\mathcal{U}_3^{\emptyset}$ and $\mathcal{U}_3^{\text{kry}}$ possess similar error behavior in the validation set, their performance as a preconditioner is quite different (cf. Table 3).
<figure><img src="image_6.png" /><figcaption>Figure 5: Convergence of train and test relative errors of the U-Net $\mathcal{U}_3^*$ when trained on different datasets, for the physical constants $k_\Omega = 10^{-3}$, $\sigma_\Omega = 10^{-3}$.</figcaption></figure>

The performance of the proposed neural preconditioner is evaluated by comparing it with two well-established preconditioning strategies: the Incomplete LU (ILU) factorization and the Algebraic Multigrid (AMG) method [37]. The ILU preconditioner approximates the LU factorization of a matrix by allowing a controlled level of fill-in, which provides a balance between sparsity and approximation accuracy, making it effective for moderately ill-conditioned systems. AMG is a hierarchical method that constructs a sequence of progressively coarser spaces to efficiently capture both low- and high-frequency error components, thereby accelerating convergence for a wide range of linear systems. Both ILU and AMG preconditioners have been implemented using the PETSc library (https://petsc.org/release/) [3], which offers efficient and scalable linear algebra routines for large-scale computational problems, together with the PETSc python wrapper `petsc4py` (https://pypi.org/project/petsc4py) [10].

All these preconditioners are tested on a flexible GMRES algorithm with right preconditioning, designated as FGMRES($k$), restarted every $k = 20$ steps. Iterations are stopped when the relative residual satisfies $||r||/||r_0|| \le 10^{-6}$. In all tables, we show the average number of iterations required by the algorithm to converge (denoted as **Mean Iter.**), for an ensemble of linear systems corresponding to $10^2$ unseen graphs $\Lambda$, ranging from $O(1)$ to $O(10^2)$ branches.

## 5.2 Effect of the data-augmentation strategy

We analyze the performance of the proposed preconditioners under three augmentation strategies: (i) no augmentation ($\mathcal{U}_3^\theta$), (ii) Krylov subspace augmentation ($\mathcal{U}_3^{kry}$) and (iii) random high-frequency vector augmentation ($\mathcal{U}_3^{hf}$). To explore the robustness of the approach, the tests are repeated for two values of the physical parameters. The results are reported in Table 3.

We discuss here the results for the most challenging case, $k_\Omega = \sigma_\Omega = 10^{-3}$. The unaugmented U-Net preconditioner ($\mathcal{U}_3^\theta$) is less effective, requiring 92.41 iterations on average, indicating its limited ability to capture the full spectrum of the solution space. Incorporating Krylov subspace augmentation ($\mathcal{U}_3^{kry}$) reduces iterations significantly (24.54), as this strategy enhances the preconditioner's ability to handle higher-order components of the operator. Random vector augmentation ($\mathcal{U}_3^{hf}$) achieves the lowest mean iteration count (17.82), demonstrating superior performance over other versions of the U-net preconditioner.

## 5.3 Impact of the 1D graph parameterization

We also investigated to what extent including the distance function $d_h^\mu$ as input enhances the effectiveness of the neural preconditioner $\mathcal{U}(\cdot, d_h^\mu; \theta)$ for different instances of 1D domains. As shown in Table 4, including the distance function $d_h^\mu$ as an input significantly improves its performance, particularly in scenarios where the
<table><caption>Table 3: Mean FGMRES iterations for different preconditioners. The variability of the data corresponds to 100 test configurations.</caption><thead><tr><th rowspan="2">Preconditioner</th><th rowspan="2">Type</th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-3</sup></th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-2</sup></th></tr><tr><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th></tr></thead><tbody><tr><td>None</td><td>-</td><td>147.72</td><td>+83.28</td><td>-52.72</td><td>187.04</td><td>+121.96</td><td>-70.04</td></tr><tr><td>U<sub>3</sub><sup>0</sup></td><td>No Augmentation</td><td>92.41</td><td>+17.59</td><td>-16.41</td><td>42.13</td><td>+10.87</td><td>-6.13</td></tr><tr><td>U<sub>3</sub><sup>kry</sup></td><td>Krylov Subspace</td><td>24.54</td><td>+3.89</td><td>-4.11</td><td>15.54</td><td>+1.46</td><td>-2.54</td></tr><tr><td>U<sub>3</sub><sup>hf</sup></td><td>Random Vectors</td><td>17.82</td><td>+5.18</td><td>-2.82</td><td>16.63</td><td>+2.37</td><td>-1.63</td></tr></tbody></table>

μ-dependent component $M_{h,00}^μ$ dominates. In fact, we recall that the structure of the linear system is defined as $C_h^μ := K_{h,00} + 2πε M_{h,00}^μ$, where $K_{h,00}$ is the fixed component, and $M_{h,00}^μ$ introduces the dependence of parameters. Without providing $d_h^μ$, the preconditioner mainly targets $K_{h,00}$, leading to lower performance on the whole system $C_h^μ$.

<table><caption>Table 4: Mean FGMRES iterations with and without influence of the distance function.</caption><thead><tr><th rowspan="2">Preconditioner</th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-3</sup></th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-2</sup></th></tr><tr><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th></tr></thead><tbody><tr><td>U<sub>3</sub><sup>0</sup>(·)</td><td>108.65</td><td>+30.35</td><td>-29.65</td><td>41.50</td><td>+8.50</td><td>-5.50</td></tr><tr><td>U<sub>3</sub><sup>0</sup>(·, d<sub>h</sub><sup>μ</sup>)</td><td>92.41</td><td>+17.59</td><td>-16.41</td><td>42.13</td><td>+10.87</td><td>-6.13</td></tr><tr><td>U<sub>3</sub><sup>kry</sup>(·)</td><td>44.91</td><td>+8.09</td><td>-7.91</td><td>20.69</td><td>+2.31</td><td>-3.69</td></tr><tr><td>U<sub>3</sub><sup>kry</sup>(·, d<sub>h</sub><sup>μ</sup>)</td><td>24.54</td><td>+3.89</td><td>-4.11</td><td>15.54</td><td>+1.46</td><td>-2.54</td></tr><tr><td>U<sub>3</sub><sup>hf</sup>(·)</td><td>60.31</td><td>+26.69</td><td>-19.31</td><td>24.98</td><td>+6.02</td><td>-4.98</td></tr><tr><td>U<sub>3</sub><sup>hf</sup>(·, d<sub>h</sub><sup>μ</sup>)</td><td>17.82</td><td>+5.18</td><td>-2.82</td><td>16.63</td><td>+2.37</td><td>-1.63</td></tr></tbody></table>

## 5.4 Effect of pre- and post-smoothing

To improve the efficacy of the preconditioning method, similar to other studies [2, 19], we can employ techniques for pre- and post-smoothing. The smoothing procedure helps mitigate high-frequency error components, particularly for preconditioners that lack high-frequency training data. For a given residual $r$ and initial guess $x_0$, we consider a smoothing algorithm with Jacobi relaxation $S_J(A, r, x_0, \text{maxit})$ [39]. The impact of pre-post smoothing is summarized in Table 5. The results show that preconditioning strategies with low-frequency training data ($U_3^0$) benefit the most from smoothing, reducing iterations by 60%.

This test sheds light on an important property of the data-augmented training strategy. Essentially, in Table 5 we observe that pre- and post-smoothing has an effect comparable to that of data augmentation. This reveals that the data-augmentation strategy is equivalent to adding a pre- and post-smoother to the preconditioner, helping in learning how to process high frequencies and contrast the spectral bias. This property may turn out to be particularly effective in those cases where an efficient smoother is not known for the problem at hand. This conclusion is further supported by the fact that the high-frequency data augmentation strategy $U_3^{hf}$ is less sensitive to smoother application.

## 5.5 Spectral analysis of training data and operator kernel

As previously observed, the introduction of high-frequency vectors in the training set is beneficial, with the augmentation strategy based on random vectors slightly outperforming that using the Kylov subspace. In
<table><caption>Table 5: Effect of pre- and post-smoothing on mean FGMRES iterations.</caption><thead><tr><th rowspan="2">Precond.</th><th rowspan="2">Relax</th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-3</sup></th><th colspan="3">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-2</sup></th></tr><tr><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th><th>Rel. Gain.</th><th>Mean Iter.</th><th>Δ<sup>+</sup></th><th>Δ<sup>-</sup></th><th>Rel. Gain</th></tr></thead><tbody><tr><td>U<sub>3</sub><sup>0</sup></td><td>None</td><td>92.41</td><td>+17.59</td><td>-16.41</td><td rowspan="2">0.64</td><td>42.13</td><td>+10.87</td><td>-6.13</td><td rowspan="2">0.47</td></tr><tr><td>U<sub>3</sub><sup>0</sup></td><td>Jacobi</td><td>33.72</td><td>+13.28</td><td>-7.72</td><td>22.47</td><td>+3.53</td><td>-4.47</td></tr><tr><td>U<sub>3</sub><sup>kry</sup></td><td>None</td><td>24.54</td><td>+3.89</td><td>-4.11</td><td rowspan="2">0.49</td><td>15.54</td><td>+1.46</td><td>-2.54</td><td rowspan="2">0.26</td></tr><tr><td>U<sub>3</sub><sup>kry</sup></td><td>Jacobi</td><td>12.6</td><td>+3.40</td><td>-2.60</td><td>11.46</td><td>+0.54</td><td>-1.46</td></tr><tr><td>U<sub>3</sub><sup>hf</sup></td><td>None</td><td>17.83</td><td>+5.17</td><td>-2.83</td><td rowspan="2">0.26</td><td>16.63</td><td>+2.37</td><td>-1.63</td><td rowspan="2">0.22</td></tr><tr><td>U<sub>3</sub><sup>hf</sup></td><td>Jacobi</td><td>13.12</td><td>+2.88</td><td>-2.12</td><td>12.96</td><td>+2.04</td><td>-0.96</td></tr></tbody></table>

this section, we seek to explain the observed behavior by contrasting the two training-set enrichment techniques with the operator's intrinsic kernel structure. In detail, to quantify the operator spectral content, we analyzed the null spaces of the problem matrices $K_{h,00}$, $2\pi\epsilon M_{00,h}^{\mu}$, and $C_h^{\mu}$ using Singular Value Decomposition (SVD). Given a matrix $A \in \mathbb{R}^{m \times n}$ with rank $r$, the associated SVD is given by $A = U\Sigma V^{\top}$, where $U = [\mathbf{u}_1 \ \mathbf{u}_2 \ \dots \ \mathbf{u}_r]$ and $V = [\mathbf{v}_1 \ \mathbf{v}_2 \ \dots \ \mathbf{v}_r]$ contain the left and right singular vectors, respectively. The right singular vectors $\mathbf{v}_i$ corresponding to zero singular values ($\sigma = 0$) span the null space, capturing kernel components that influence preconditioning performance. Using sparse SVD algorithms (scipy sparse svds), we extracted some of these vectors and analyzed them in frequency space via the multidimensional discrete Fourier transform. Let us denote by $\mathbf{V}$ the $d$-dimensional tensor corresponding to a right singular vector $\mathbf{v}$. Then we denote the DFT as $\mathcal{F}$, defined as follows,

$$ \mathcal{F}: \mathbf{V} \in \mathbb{R}_{n_1} \times \dots \times \mathbb{R}_{n_d} \rightarrow \mathcal{F}(\mathbf{V}) \in \mathbb{C}_{k_1} \times \dots \times \mathbb{C}_{k_d}, \quad \mathcal{F}(\mathbf{V})|_{\mathbf{k}} = \sum_{\mathbf{n}=0}^{N-1} e^{-i2\pi\mathbf{k} \cdot (\mathbf{n}/N)} \mathbf{V}|_{\mathbf{n}}, $$

where $\mathbf{n} = (n_1, n_2, \dots, n_d) \in \mathbb{N}^d$, $\mathbf{k} = (k_1, k_2, \dots, k_d) \in \mathbb{N}^d$ and $\mathbf{N} = (N, N, \dots, N) \in \mathbb{N}^d$ are multi-integers and is $\mathbf{N}$ the number of edge subdivision for the considered cubic structured mesh.

In Figure 6 (columns I-III), two-dimensional slices of the spectrum magnitudes $||\mathcal{F}(\mathbf{V}_i)||$ associated with the kernel vectors $\{\mathbf{v}_i\}_{i=1,2,3}$ of matrices $K_{h,00}$, $2\pi\epsilon M_{h,00}^{\mu}$ and $C_h^{\mu}$ are shown. For better visualization, the so-called centered spectrum is considered, where a suitable phase shift $e^{i\pi\mathbf{k}(N-1)}$ in the frequency space is applied to $\mathcal{F}(\mathbf{V})$; in this way, the point $\mathcal{F}(\mathbf{V})|_{(0,...,0)}$ is centered on the plot. We note that the frequency content in the spectra of matrices $K_{h,00}$ and $2\pi\epsilon M_{h,00}^{\mu}$ is remarkably different. As expected, the kernel of the metric term $2\pi\epsilon M_{h,00}^{\mu}$ introduces very noisy high-frequency components that need to be preconditioned. Figure 6 (columns IV-VI) presents a comparison between the different data enhancement strategies, visualized in the frequency domain. Column IV shows the spectrum magnitude of the normalized forcing term. Column V illustrates the spectrum magnitude of the orthonormal Krylov basis. Finally, the last column presents the spectrum magnitude of the random high-frequency data augmentation set. In Figure 6 we observe that Krylov subspace augmentation captures higher-order components by repeated application of the operator to a basis vector, producing structured high-frequency content (column V). In contrast, random vector augmentation introduces unstructured high-frequency content by randomly drawing components in an independent manner. This results in a broader and more uniform spectral coverage (Fig. 6, column VI), enabling the preconditioner to handle fine-scale variations in the solution, critical for mixed-dimensional problems like the 3D-1D coupled system. From a comparison between the spectral content of different training sets and the spectrum of the kernel component of the metric term $2\pi\epsilon M_{h,00}^{\mu}$ (Fig. 6, column II), one can conclude that the Krylov subspace augmentation introduces structured high-frequency components, while random vector augmentation spans a broader frequency range, closer to the true kernel spectrum of the operator. This broader frequency coverage can explain the better convergence behavior observed in numerical experiments, where random vector-based augmentation results are more effective than Krylov-based strategies in handling ill-conditioning and diverse parameter variations.
<figure><img src="image_7.png" /><figcaption>Figure 6: Logarithmic plot in the frequency domain ($k_1, k_2, k_3$) of the $k_2k_3$-cross section for the spectrum magnitude $\|\mathcal{F}(\mathbf{V}_i)\|$; $k_1 = 0$ is considered and a value for $\mu$ is fixed given $k_\Omega = \sigma_\Omega = 10^{-3}$. Frequency components increase with the distance from the center (0, 0, 0). **Column I:** magnitude of first three right-singular vectors associated to $K_{h,00}$; **Column II:** magnitude of first three right-singular vectors associated to $2\pi\epsilon M_{h,00}^{\mu}$; **Column III:** magnitude of first three right-singular vectors associated to $C_h^\mu = K_{h,00} + 2\pi\epsilon M_{h,00}^{\mu}$; **Column IV:** spectrum magnitude of the normalized forcing term $b_h^\mu$; **Column V:** spectrum magnitude of the orthonormal Krylov bases $\mathcal{D}^{kry,\mu} := \{q_5^\mu, \dots, q_8^\mu\}$ (note that $q_0^\mu \equiv b_h^\mu / \|b_h^\mu\|$); **Column VI:** spectrum magnitude of the random high frequency data augmentation set $\mathcal{D}^{hf,\mu} := \{r_{h,1} \dots r_{h,4}\}$.</figcaption></figure>

## 5.6 Scalability with respect to the problem size

Due to the intrinsic properties of the convolutional layers in the U-Net architecture, although trained on a low-resolution input tensor $\mathbf{X} \in \mathbb{R}^{21 \times 21 \times 21 \times 2}$, the proposed neural preconditioner $U_3^{hf}$ can be applied to inputs of arbitrary dimensions, corresponding to higher resolutions of the problem. This adaptability makes it suitable for evaluating its scalability with respect to increasing problem size, quantified by the number of unknowns $N_h$ of the discrete problem.

We test here the convergence rate of the preconditioned FGMRES for the solution of the 3D problem, for three different levels of mesh resolution $N$, namely [21, 41, 81] points per side of the cubic domain $\Omega$ corresponding respectively to a number of degrees of freedom (d.o.f.), $N_h$, for a $\mathbb{P}^1$ FEM discretization equal to [21³; 41³; 81³]. We measure the average number of iterations across 100 different test cases, all referring to 1D graphs that were not used during the training phase. Then, we analyze how the average iteration count scales with the number of d.o.f., $N_h$.

To assess the capacity of the neural preconditioner to generalize on the mesh resolution, we compare the neural preconditioner $U_3^{hf}$ trained on two different resolutions, $N = 21$ and $N = 41$, respectively, and analyze how these two operators behave in the case of discretization based on 81 points. The comparison also includes incomplete LU factorization (ILU) and Algebraic Multigrid (AMG) with 3 and 10 V-cycles, respectively.
<table><caption>Table 6: Mean FGMRES iterations and execution time for different preconditioning strategies of the problem with $k_{\Omega} = \sigma_{\Omega} = 10^{-3}$. All times are in milliseconds. (*) Denotes that FGMRES did not converge in the maximum number of iterations.</caption><thead><tr><th>Precond.</th><th>N<sub>h</sub></th><th>N</th><th>Mean Iter.</th><th>Rate (N<sub>h</sub>)</th><th>Time</th><th>Rate (N<sub>h</sub>)</th><th>Time/Iter</th><th>Rate (N<sub>h</sub>)</th></tr></thead><tbody><tr><td rowspan="3">None</td><td>9,261</td><td>21</td><td>147.72</td><td>-</td><td>365.33</td><td>-</td><td>-</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>431.91</td><td>0.53</td><td>7293.49</td><td>1.49</td><td>-</td><td>-</td></tr><tr><td>531,441</td><td>81</td><td>1000*</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">U<sub>3</sub><sup>hf</sup>(21<sup>3</sup>)</td><td>9,261</td><td>21</td><td>17.38</td><td>-</td><td>45.227</td><td>-</td><td>2.54</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>37.24</td><td>0.37</td><td>497.86</td><td>1.2</td><td>13.37</td><td>1.21</td></tr><tr><td>531,441</td><td>81</td><td>100.73</td><td>0.49</td><td>9212.77</td><td>1.43</td><td>91.46</td><td>1.06</td></tr><tr><td rowspan="3">U<sub>3</sub><sup>hf</sup>(41<sup>3</sup>)</td><td>9,261</td><td>21</td><td>28.81</td><td>-</td><td>60.54</td><td>-</td><td>2.10</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>45.54</td><td>0.27</td><td>519.95</td><td>1.07</td><td>11.42</td><td>1.19</td></tr><tr><td>531,441</td><td>81</td><td>119.27</td><td>0.47</td><td>11985.8</td><td>1.54</td><td>100.49</td><td>0.94</td></tr><tr><td rowspan="3">ILU</td><td>9,261</td><td>21</td><td>24.73</td><td>-</td><td>37.7</td><td>-</td><td>1.52</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>46.2</td><td>0.31</td><td>486.04</td><td>1.27</td><td>10.52</td><td>1.04</td></tr><tr><td>531,441</td><td>81</td><td>102.55</td><td>0.39</td><td>9422.78</td><td>1.45</td><td>91.88</td><td>0.94</td></tr><tr><td rowspan="3">AMG(3)</td><td>9,261</td><td>21</td><td>9.56</td><td>-</td><td>62.35</td><td>-</td><td>6.52</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>17.13</td><td>0.29</td><td>696.14</td><td>1.20</td><td>40.64</td><td>1.10</td></tr><tr><td>531,441</td><td>81</td><td>33.66</td><td>0.33</td><td>9631.04</td><td>1.29</td><td>286.13</td><td>1.05</td></tr><tr><td rowspan="3">AMG(10)</td><td>9,261</td><td>21</td><td>3.64</td><td>-</td><td>55.67</td><td>-</td><td>15.29</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>4.01</td><td>0.05</td><td>329.32</td><td>1.02</td><td>82.12</td><td>1.19</td></tr><tr><td>531,441</td><td>81</td><td>4.9</td><td>0.10</td><td>3921.46</td><td>1.07</td><td>800.30</td><td>0.90</td></tr></tbody></table>

In Table 6 and 7, we see that the scaling behavior of the neural preconditioner $U_3^{hf}$ closely matches that of the ILU preconditioner. These methods show a dependence on $N_h$ of order $\simeq 0.4 - 0.5$. The AMG(3) preconditioner shows a lower iteration count but features a similar scaling with respect to the size of the problem. None of these methods achieves optimal scaling as the AMG(10) preconditioner, the iteration rates of which are almost independent on $N_h$.

Finally, we note that the performance of the neural preconditioners $U_3^{hf}(21^3)$ and $U_3^{hf}(41^3)$ is almost equivalent, indicating that the size and resolution of the training data do not affect the overall performance of the neural operator. This means that it is possible to train the neural preconditioner on low-resolution data and deploy it on higher resolutions, with considerable savings on the computational cost of the offline training phase.

Previous analyses have focused solely on the number of FGMRES iterations as a performance metric. Although the neural preconditioners demonstrated sub-optimal results in terms of convergence rates, it is also essential to consider execution time, as the forward pass through a neural network is computationally lighter than, e.g., a single evaluation of the AMG algorithm. Although the implementation has not been optimized for performance (e.g., no batch-optimized algorithm or GPU acceleration techniques [31]), we also present in Tables 6-7 the results for execution time. Looking at the scaling with $N_h$ of the mean time required to solve the 3D problem, the advantage of AMG(10) is significantly reduced with respect to the neural preconditioner, ILU, and AMG(3), which provide aligned performances. This can be interpreted by observing that for FGMRES with sparse matrices, the overall computational cost is $O(m \times N_h)$, where m is the number of iterations required for convergence. On the one hand, m is sensitive to the choice of the preconditioner; precisely as shown in Table 6-7 column Mean Iter. we have $m = (N_h)^{p_i}$ which being $p_i \simeq 0.5$ for $U_3^{hf}$, $p_i \simeq 0.4$ for ILU, $p_i \simeq 0.4$ for AMG(3) and $p_i \simeq 0$ for AMG (10), this being the only case where m is independent of $N_h$. On the other hand, the time per iteration scales linearly with $N_h$, that is, $T = c_i N_h$ for all methods, but the proportionality constant may be radically different, namely $c_i = 1.72310^{-4}$, $1.88310^{-4}$ for neural preconditioners trained on points 21³ and 41³, $c_i = 1.72110^{-4}$ for ILU, $c_i = 5.38310^{-4}$ for AMG(3) and $c_i = 14.9810^{-4}$ for AMG(10), with neural and ILU preconditioners showing advantage over AMG. As a result, the following relation between the total time and the size of the problem holds, $\log T = c_i + (p_i + 1) \log N_h$,
<table><caption>Table 7: Mean FGMRES iterations and execution time for different preconditioning strategies of the problem with $k_\Omega = \sigma_\Omega = 10^{-2}$. All times are in milliseconds. (*) Denotes that FGMRES did not converge in the maximum number of iterations.</caption><thead><tr><th>Precond.</th><th>N<sub>h</sub></th><th>N</th><th>Mean Iter.</th><th>Rate (N<sub>h</sub>)</th><th>Time</th><th>Rate (N<sub>h</sub>)</th><th>Time/Iter</th><th>Rate (N<sub>h</sub>)</th></tr></thead><tbody><tr><td rowspan="3">None</td><td>9,261</td><td>21</td><td>187.94</td><td>-</td><td>216.17</td><td>-</td><td>-</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>511.34</td><td>0.50</td><td>3,329.16</td><td>1.36</td><td>-</td><td>-</td></tr><tr><td>531,441</td><td>81</td><td>1000*</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">U<sub>3</sub><sup>hf</sup>(21<sup>3</sup>)</td><td>9,261</td><td>21</td><td>16.63</td><td>-</td><td>41.18</td><td>-</td><td>2.48</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>37.98</td><td>0.41</td><td>461.18</td><td>1.20</td><td>12.14</td><td>1.26</td></tr><tr><td>531,441</td><td>81</td><td>113.98</td><td>0.54</td><td>9,908.77</td><td>1.50</td><td>86.93</td><td>1.04</td></tr><tr><td rowspan="3">U<sub>3</sub><sup>hf</sup>(41<sup>3</sup>)</td><td>9,261</td><td>21</td><td>24.27</td><td>-</td><td>55.73</td><td>-</td><td>2.30</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>41.39</td><td>0.27</td><td>458.34</td><td>1.05</td><td>11.07</td><td>1.28</td></tr><tr><td>531,441</td><td>81</td><td>104.99</td><td>0.46</td><td>9,586.00</td><td>1.49</td><td>91.30</td><td>0.97</td></tr><tr><td rowspan="3">ILU</td><td>9,261</td><td>21</td><td>31.46</td><td>-</td><td>42.09</td><td>-</td><td>1.34</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>62.18</td><td>0.34</td><td>548.30</td><td>1.28</td><td>8.82</td><td>1.06</td></tr><tr><td>531,441</td><td>81</td><td>147.30</td><td>0.42</td><td>12,123.32</td><td>1.52</td><td>82.30</td><td>0.91</td></tr><tr><td rowspan="3">AMG(3)</td><td>9,261</td><td>21</td><td>11.50</td><td>-</td><td>63.95</td><td>-</td><td>5.56</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>20.69</td><td>0.29</td><td>727.50</td><td>1.21</td><td>35.16</td><td>1.09</td></tr><tr><td>531,441</td><td>81</td><td>42.14</td><td>0.35</td><td>10,839.18</td><td>1.32</td><td>257.22</td><td>1.03</td></tr><tr><td rowspan="3">AMG(10)</td><td>9,261</td><td>21</td><td>3.92</td><td>-</td><td>53.25</td><td>-</td><td>13.58</td><td>-</td></tr><tr><td>68,921</td><td>41</td><td>4.00</td><td>0.01</td><td>404.91</td><td>1.01</td><td>101.23</td><td>0.999</td></tr><tr><td>531,441</td><td>81</td><td>4.00</td><td>0.00</td><td>3,450.91</td><td>1.05</td><td>862.73</td><td>0.95</td></tr></tbody></table>

showing that methods such as neural preconditioners and ILU provide a potential advantage for systems of small size.

We finally observe that the neural preconditioner struggles to address large systems, for example, the one arising from a discretization of the 3D domain with 101 points per side, corresponding to $N_h \simeq 10^6$, for which $U_3^{hf}(21^3)$ requires on average 188.65 iterations to solve the 3D problem. The behavior does not improve with increasing the dimension of the training set, as this is true for both $U_3^{hf}(21^3)$ and $U_3^{hf}(41^3)$. We hypothesize that this is due to the effect of high frequencies that become dominant components when the discrete solution belongs to a high-dimensional space. For this reason, we tested the preconditioner $U_3^{hf}(21^3)$ with pre- and post-smoothing, as described in Section 5.4, with very encouraging results given by an average number of iterations equal to [13.12; 27.60; 62.32; 91.55] for $N_h = [21^3; 41^3; 81^3; 101^3]$, respectively. This corresponds to the parameters $p_i = 0.31$ for the dependence of the number of iterations on $N_h$ and $c_i = 2.30210^{-4}$ for the proportionality constant between time and iterations. This performance is in line with the ILU preconditioner, actually slightly better.

## 5.7 Coupling the neural preconditioner with ensembled calculations

In the context of parametric simulations of mixed-dimensional PDEs, solving multiple instances of the underlying linear system for varying configurations of the low-dimensional structure can become computationally prohibitive. To mitigate this, the so-called *ensembled* framework has recently been developed to accelerate the solution of families of linear systems that share a common structure [33, 25]. Ensembled methods consist of stacking multiple coefficient matrices and right-hand side vectors, each corresponding to a different parameter instance, and solving them simultaneously as a batched linear system. This paradigm is especially well-suited for settings such as uncertainty quantification or parametric sweeps, where many problems must be solved independently but share structural similarities. The main advantage of this approach lies in the optimized memory access.

We note that a key feature of the proposed neural preconditioner lies in its ability to generalize across the parameter space and that this property is well suited to ensemble methods, as the neural operator can process a group of inputs at the same time, each representing various parameter settings.

We have implemented a prototype of the *ensembled* calculation using the algorithmic infrastructure devel-
<figure><img src="image_8.png" /><figcaption>Figure 7: Left: plot of the speed-up associated to the batched execution of the neural preconditioning pass, using $U_3^{hf}(\mathbf{X}_e)$ with $\mathbf{X}_e \in \mathbb{R}^{N \times 2 \times 21 \times 21 \times 21}$, and batch size $N \in [1, 100]$. Right: Plot of the execution time, expressed in millisecond, associated respectively to the serial (non-batched) and batched preconditioner pass of $U_3^{hf}$.</figcaption></figure>

oped for the neural preconditioner, based on the PyTorch library ([32]) together with the CUDA kernel ([27]). In fact, we recall that the input of the neural preconditioner consists of tensors $\mathbf{X}_i \in \mathbb{R}^{c \times n}$ with $c$ representing the number of channels and $n = n_1 \times n_2 \times \dots \times n_d$ a multi-index encoding the physical dimension of the data (see Section 4.2). Given a set of tensors $\{\mathbf{X}_i\}_{i=1..N}$, it is possible to stack them along an additional batch dimension to obtain a new tensor $\mathbf{X}_e = [\mathbf{X}_1|\dots|\mathbf{X}_N] \in \mathbb{R}^{N \times c \times n}$ that encodes all input tensors to be fed to the neural preconditioner.

We evaluated the speed-up factor achieved by implementing the neural preconditioner $U_3^{hf}$ on $N \in [1, 100]$ input tensors $\mathbf{X}_i \in \mathbb{R}^{2 \times 21 \times 21 \times 21}$, which correspond to $N$ different 1D graphs in eq. 9, i.e.

$$s = \frac{\sum_{i=1}^{N} t(U_3^{hf}(\mathbf{X}_i))}{t(U_3^{hf}(\mathbf{X}_e))}$$

with $t(U_3^{hf}(\cdot))$ representing the execution time of the neural preconditioner pass. The numerical tests are performed on an AMD EPYC 7301 CPU and a NVIDIA Tesla V100 GPU with 64 MB of dedicated RAM.

From the results shown in Table 8 and Fig. 7, we observe a speed-up around 4.5 when up to 100 input tensors are considered. Saturation occurs as the batched input approaches the memory limit of the hardware. The beneficial acceleration is associated with the contiguous memory structure, coupled with CUDA streams and execution via the PyTorch library. We prospect that the synergy between the ensembled approach and the neural preconditioner, along with the use of GPU based architectures, can enable significant gains in solver throughput, especially in multi-query scenarios. A comprehensive study encompassing computational architecture, algorithmic design, memory management, and preconditioner implementation is deferred to future research.

<table><caption>Table 8: Execution times [ms] and speed-up of the forward pass through the neural preconditioner as a function of the batch size.</caption><thead><tr><th>Batch size (N)</th><th>1</th><th>5</th><th>10</th><th>15</th><th>20</th><th>30</th><th>50</th><th>70</th><th>90</th><th>100</th></tr></thead><tbody><tr><th>time non-batched [ms]</th><td>1.5</td><td>6.1</td><td>11.8</td><td>17.5</td><td>23.2</td><td>34.7</td><td>58.2</td><td>81.3</td><td>104.7</td><td>114.9</td></tr><tr><th>time batched [ms]</th><td>1.5</td><td>2.6</td><td>4.1</td><td>5.5</td><td>6.6</td><td>9.7</td><td>13.5</td><td>18.2</td><td>22.9</td><td>25.5</td></tr><tr><th>Speed-up</th><td>1.0</td><td>2.3</td><td>2.9</td><td>3.5</td><td>3.7</td><td>3.6</td><td>4.3</td><td>4.5</td><td>4.6</td><td>4.5</td></tr></tbody></table>
# 6 Solution of the 3D-1D coupled problem

We now focus on the fully coupled 3D-1D problem, assessing how well the suggested block-preconditioning strategy (7)-(8) performs when applied to the linear system (6). We use a right-preconditioned FGMRES solver with a restart parameter $k = 20$ and a stopping criterion $||r||/||r_0|| \le 1 \cdot 10^{-15}$. This iterative solver is preconditioned by an approximation $\tilde{Q}_h^\mu$ of the block preconditioner defined in (8). More precisely, the block preconditioner $\tilde{Q}_h^\mu$ is defined by an approximate solution to the 1D problem, obtained with ILU factorization, returning $\tilde{z}_{h,1}^\mu$ (a direct solver has also been tested, leading to no remarkable differences); then, 3 steps of the preconditioned FGMRES algorithm (precisely P-FGMRES with various preconditioners P) return $\tilde{z}_{h,0}^\mu$, an approximate solution of the 3D system.

The performance of this approach for the 3D-1D coupled problem (6) with a comparison of monolithic solvers and different choices of preconditioners for the 3D problem can be found in Table 9. The latter presents a comparative analysis of the different strategies employed to solve the coupled 3D-1D system using the FGMRES iterative solver, with a focus on the impact of preconditioning.

<table><caption>Table 9: Mean FGMRES iterations and execution times with associated rates for the preconditioning strategies considered, in different mesh configurations. Times are reported in seconds. Rates concerning iteration counts and execution times are computed with respect to $N_h$. (*) Denotes that the maximum number of iterations set to 2000 was exceeded. (n.r.) means that the data was not recorded.</caption><thead><tr><th rowspan="2">Precond.</th><th rowspan="2">N<sub>h</sub></th><th colspan="4">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-3</sup></th><th colspan="4">k<sub>Ω</sub> = σ<sub>Ω</sub> = 10<sup>-2</sup></th></tr><tr><th>Iter</th><th>Rate</th><th>Time</th><th>Rate</th><th>Iter</th><th>Rate</th><th>Time</th><th>Rate</th></tr></thead><tbody><tr><td rowspan="3">Non prec.<br>3D-1D problem</td><td>9,261</td><td>1973.55*</td><td></td><td></td><td></td><td>1768.65*</td><td></td><td></td><td></td></tr><tr><td>68,921</td><td>2000*</td><td></td><td>n.r.</td><td></td><td>2000*</td><td></td><td>n.r.</td><td></td></tr><tr><td>531441</td><td>2000*</td><td></td><td></td><td></td><td>2000*</td><td></td><td></td><td></td></tr><tr><td rowspan="3">AMG(10)<br>monolithic</td><td>9,261</td><td>1181.48*</td><td></td><td></td><td></td><td>1220.83*</td><td></td><td></td><td></td></tr><tr><td>68,921</td><td>1366.02*</td><td></td><td>n.r.</td><td></td><td>1423.44*</td><td></td><td>n.r.</td><td></td></tr><tr><td>531441</td><td>1862.71*</td><td></td><td></td><td></td><td>1918.29*</td><td></td><td></td><td></td></tr><tr><td rowspan="3">Haznics-AMG(10)<br>3D-1D problem</td><td>9,261</td><td>38.71</td><td>-</td><td>1.47</td><td>-</td><td>35.36</td><td>-</td><td>1.49</td><td>-</td></tr><tr><td>68,921</td><td>48.04</td><td>0.11</td><td>2.14</td><td>0.19</td><td>41.39</td><td>0.08</td><td>2.26</td><td>0.21</td></tr><tr><td>531441</td><td>725.44*</td><td>1.35</td><td>84.232</td><td>1.80</td><td>497.57*</td><td>1.22</td><td>64.95</td><td>1.64</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = none</td><td>9,261</td><td>153.82</td><td>-</td><td>0.50</td><td>-</td><td>123.25</td><td>-</td><td>0.49</td><td>-</td></tr><tr><td>68,921</td><td>417.82</td><td>0.50</td><td>8.51</td><td>1.42</td><td>316.88</td><td>0.47</td><td>4.48</td><td>1.10</td></tr><tr><td>531441</td><td>1564.84</td><td>0.65</td><td>260.41</td><td>1.32</td><td>916.19</td><td>0.52</td><td>152.88</td><td>1.73</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = U<sub>3</sub><sup>hf</sup>(21<sup>3</sup>)</td><td>9,261</td><td>60.74</td><td>-</td><td>0.47</td><td>-</td><td>22.28</td><td>-</td><td>0.27</td><td>-</td></tr><tr><td>68,921</td><td>123.30</td><td>0.35</td><td>3.12</td><td>0.94</td><td>54.33</td><td>0.44</td><td>1.54</td><td>0.87</td></tr><tr><td>531441</td><td>438.69</td><td>0.62</td><td>93.49</td><td>1.66</td><td>157.69</td><td>0.52</td><td>32.98</td><td>1.50</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = U<sub>3</sub><sup>hf</sup>(41<sup>3</sup>)</td><td>9,261</td><td>69.15</td><td>-</td><td>0.51</td><td>-</td><td>26.95</td><td>-</td><td>0.29</td><td>-</td></tr><tr><td>68,921</td><td>104.80</td><td>0.21</td><td>2.68</td><td>0.83</td><td>36.55</td><td>0.15</td><td>1.06</td><td>0.64</td></tr><tr><td>531441</td><td>257.91</td><td>0.44</td><td>54.62</td><td>1.48</td><td>59.99</td><td>0.23</td><td>12.78</td><td>1.22</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = AMG(10)</td><td>9,261</td><td>50.19</td><td>-</td><td>0.72</td><td>-</td><td>16.85</td><td>-</td><td>0.35</td><td>-</td></tr><tr><td>68,921</td><td>52.57</td><td>0.02</td><td>5.40</td><td>1.00</td><td>17.84</td><td>0.002</td><td>1.95</td><td>0.61</td></tr><tr><td>531441</td><td>56.13</td><td>0.03</td><td>49.17</td><td>1.08</td><td>18.85</td><td>0.02</td><td>17.49</td><td>0.81</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = AMG(3)</td><td>9,261</td><td>51.53</td><td>-</td><td>0.63</td><td>-</td><td>18.23</td><td>-</td><td>0.33</td><td>-</td></tr><tr><td>68,921</td><td>56.61</td><td>0.05</td><td>4.64</td><td>0.99</td><td>23.00</td><td>0.12</td><td>2.01</td><td>0.90</td></tr><tr><td>531441</td><td>68.64</td><td>0.08</td><td>45.58</td><td>1.12</td><td>37.1</td><td>0.23</td><td>24.75</td><td>1.23</td></tr><tr><td rowspan="3">1D: ILU<br>3D: P-FGMRES<br>P = ILU</td><td>9,261</td><td>52.14</td><td>-</td><td>0.4</td><td>-</td><td>18.72</td><td>-</td><td>0.26</td><td>-</td></tr><tr><td>68,921</td><td>58.08</td><td>0.05</td><td>1.97</td><td>0.79</td><td>26.52</td><td>0.17</td><td>1.05</td><td>0.70</td></tr><tr><td>531441</td><td>76.33</td><td>0.13</td><td>20.32</td><td>1.14</td><td>47.83</td><td>0.29</td><td>13.25</td><td>1.24</td></tr></tbody></table>
As seen in the table, the results clearly highlight the importance of preconditioning for the efficient resolution of mixed-dimensional PDEs. Indeed, the first key observation is that the absence of a preconditioner results in a solver failing to converge within the budget of 2000 iterations. This emphasizes the severity of the ill-conditioning that characterizes the 3D-1D coupled problem and establishes the necessity of a robust preconditioning strategy. The block-diagonal preconditioner $\tilde{Q}_h^\mu$ emerges as a good option, underlining the effectiveness of a decoupled block-structured approach that aligns with the intrinsic structure of the mixed-dimensional system. Interestingly, the full monolithic preconditioning strategy based on algebraic multigrid with 10 smoothing steps (denoted AMG(10) monolithic) fails to deliver competitive performance. This suggests that generic monolithic approaches are suboptimal for problems with strong 3D-1D coupling, as they are not designed to exploit the specific features of the system. In contrast, the Haznics preconditioner, which involves a combination of AMG in the entire 3D domain and Schwarz smoother on the interface between the 3D and 1D domain [9], shows the best performance in terms of iterations, validating its role as the current state-of-the-art method for mixed-dimensional problems. A degradation of the mean iteration count can be observed for the highest value of $N_h$, especially due to lack of convergence when very complex graph geometries are considered.

The table also reports the results for strategies that combine the block-diagonal preconditioner $\tilde{Q}_h^\mu$ with various solvers applied to its diagonal blocks. In particular, the first block, corresponding to the 3D subproblem, is preconditioned with ILU, AMG (3), AMG (10) or different instances of the neural preconditioner, namely $U_3^{hf}(21^3)$ and $U_3^{hf}(41^3)$. Although the iteration counts vary distinctively between the different strategies, we notice that their overall computational times are remarkably similar. This is a result of a trade-off between the cost per iteration and the number of iterations required. For instance, classical preconditioners like AMG generally result in lower iteration counts but at a higher per-iteration cost. Neural preconditioners, on the other hand, often yield higher iteration counts but with significantly lower cost per iteration, because of their efficient matrix-free representation. This balance ensures that all these strategies achieve similar total run-times, with performance comparable to or superior to monolithic alternatives.

These findings support the potential of the neural preconditioners as a viable and efficient strategy for solving complex mixed-dimensional problems of small-medium scale. The flexibility in trading iteration count for speed suggests that these learned preconditioners can serve as lightweight generalizable components in the broader preconditioning pipeline.

# 7 Conclusions and Future Perspectives

In this work, we have developed a novel unsupervised neural network-based preconditioning strategy tailored to efficiently solve mixed-dimensional PDEs, specifically targeting the computational challenges arising from the coupling between 3D domains and embedded 1D structures. Using operator learning through convolutional neural networks, our approach successfully generalizes across various topological configurations of the 1D graphs and adapts robustly to different mesh resolutions without requiring retraining.

The neural preconditioner was extensively validated against classical preconditioning methods, including incomplete LU factorization and algebraic multigrid methods. Numerical experiments demonstrated that the proposed neural network preconditioner significantly accelerates convergence of iterative solvers, such as the Flexible GMRES, effectively handling both low- and high-frequency error components. Among the data augmentation strategies investigated, enriching the training set with random high-frequency vectors emerged as the most effective, enabling the neural network to overcome the inherent spectral bias common to most machine learning-based approaches.

Furthermore, our analysis underscored the critical role of including explicit parametric information through the distance function associated with the 1D graph, which improved the effectiveness of the preconditioner in varying geometric configurations. We also explored the complementary effects of pre- and post-smoothing strategies, illustrating their capacity to enhance preconditioning performance, particularly when the training set lacks high-frequency data.

However, there are limitations and areas for future improvement. A significant current limitation is that the proposed approach is applicable only to tensor-product domains, which restricts its use in realistic applications. To address this, we are actively exploring the use of graph- and mesh-informed neural networks [4, 13], as suitable alternatives to convolutional layers, with the aim of extending the applicability of this method to more general computational domains. Additionally, developing a solid theoretical framework for
selecting data augmentation subsets and, more broadly, for learning-based preconditioning methods is crucial
to consolidate this approach and broaden the range of application areas. Although the scalability of the neural
preconditioner shows promise, more work is required to enhance its performance in large-scale problems.

The promising results obtained from this study not only establish neural preconditioning as a viable
and efficient alternative to traditional preconditioners for mixed-dimensional PDEs, but also lay a robust
foundation for further exploration into broader classes of coupled multiphysics problems.

Acknowledgments

The authors thank Prof. Kent-Andre Mardal and Dr. Miroslav Kuchta of the Simula Research Laboratory for useful
discussions. PZ acknowledges the support of the MUR PRIN 2022 grant No. 2022WKWZA8 *Immersed methods for
multiscale and multiphysics problems* (IMMEDIATE) part of the Next Generation EU program Mission 4, Comp.
2, CUP D53D23006010006. NF has received funding under the project *Reduced Order Modeling and Deep Learning
for the real-time approximation of PDEs (DREAM)*, grant no. FIS00003154, funded by the Italian Science Fund
(FIS) and by Ministero dell'Università e della Ricerca (MUR). The present research is part of the activities of the
*Dipartimento di Eccellenza* 2023-2027, Department of Mathematics, Politecnico di Milano. All authors are members
of the Gruppo Nazionale per il Calcolo Scientifico (GNCS), Istituto Nazionale di Alta Matematica (INdAM).

References

[1] Ben Adcock, Simone Brugiapaglia, and Clayton G. Webster. *Sparse Polynomial Approximation of High-Dimensional Functions*. SIAM, Philadelphia, PA, 2022.

[2] Yael Azulay and Eran Treister. Multigrid-augmented deep learning preconditioners for the helmholtz equation. *SIAM Journal on Scientific Computing*, 45(3):S127-S151, 2023.

[3] S Balay, S Abhyankar, Mark F Adams, S Benson, J Brown, P Brune, K Buschelman, EM Constantinescu, L Dalcin, A Dener, et al. Petsc/tao users manual (rev. 3.20). Technical report, Argonne National Laboratory (ANL), Argonne, IL (United States), 2023.

[4] Peter Battaglia, Jessica Blake Chandler Hamrick, Victor Bapst, Alvaro Sanchez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Caglar Gulcehre, Francis Song, Andy Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey Allen, Charles Nash, Victoria Jayne Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet Kohli, Matt Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu. Relational inductive biases, deep learning, and graph networks. *arXiv*, 2018.

[5] Wietse M. Boon, Jan M. Nordbotten, and Jon E. Vatne. Functional analysis and exterior calculus on mixed-dimensional geometries. *Ann. Mat. Pura Appl.*, 200(2):757-789, 2021.

[6] Nicolas Boullé and Alex Townsend. Chapter 3 - a mathematical guide to operator learning. In Siddhartha Mishra and Alex Townsend, editors, *Numerical Analysis Meets Machine Learning*, volume 25 of *Handbook of Numerical Analysis*, pages 83–125. Elsevier, 2024.

[7] Alessio Paolo Buccino, Miroslav Kuchta, Jakob Schreiner, and Kent-André Mardal. Improving neural simulations with the emi model. In Aslak Tveito, Kent-Andre Mardal, and Marie E. Rognes, editors, *Modeling Excitable Tissue: The EMI Framework*, pages 87–98. Springer International Publishing, Cham, 2021.

[8] Ana Budiša, Xiaozhe Hu, Miroslav Kuchta, Kent-Andre Mardal, and Ludmil Zikatanov. Algebraic multigrid methods for metric-perturbed coupled problems. *SIAM Journal on Scientific Computing*, 46(3):A1461–A1486, 2024.

[9] Ana Budiša, Xiaozhe Hu, Miroslav Kuchta, Kent-André Mardal, and Ludmil T. Zikatanov. Haznics – software components for multiphysics problems. *ACM Trans. Math. Softw.*, 49(4), December 2023.
[10] Lisandro D Dalcin, Rodrigo R Paz, Pablo A Kler, and Alejandro Cosimo. Parallel distributed computing using python. *Advances in Water Resources*, 34(9):1124–1139, 2011.

[11] J. Dugundji. An extension of Tietze's theorem. *Pacific Journal of Mathematics*, 1(3):353–367, 1951.

[12] Max Firmbach, Ivo Steinbrecher, Alexander Popp, and Matthias Mayr. An approximate block factorization preconditioner for mixed-dimensional beam-solid interaction. *Computer Methods in Applied Mechanics and Engineering*, 431:117256, 2024.

[13] Nicola Rares Franco, Andrea Manzoni, and Paolo Zunino. Mesh-informed neural networks for operator learning in finite element spaces. *Journal of Scientific Computing*, 97(35), 2023.

[14] Alessio Fumagalli and Anna Scotti. A numerical method for two-phase flow in fractured porous media with non-matching grids. *Advances in Water Resources*, 62:454 – 464, 2013.

[15] Moritz Geist, Philipp Petersen, Mones Raslan, Reinhold Schneider, and Gitta Kutyniok. Numerical solution of the parametric diffusion equation by deep neural networks. *Journal of Scientific Computing*, 88(1):22, 2021.

[16] Gjerde, Ingeborg G., Kumar, Kundan, Nordbotten, Jan M., and Wohlmuth, Barbara. Splitting method for elliptic equations with line sources. *ESAIM: M2AN*, 53(5):1715–1739, 2019.

[17] Viktor Grimm, Alexander Heinlein, and Axel Klawonn. Learning the solution operator of two-dimensional incompressible navier-stokes equations using physics-aware convolutional neural networks. *arXiv*, math.NA(2308.02137), 2023.

[18] Nora Hagmeyer, Matthias Mayr, Ivo Steinbrecher, and Alexander Popp. One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa. *Advanced Modeling and Simulation in Engineering Sciences*, 9(1), 2022.

[19] Alena Kopaničáková and George Em Karniadakis. Deeponet based preconditioning strategies for solving parametric linear systems of equations. *SIAM Journal on Scientific Computing*, 47(1):C151–C181, 2025.

[20] Nikola B. Kovachki, Samuel Lanthaler, and Andrew M. Stuart. Chapter 9 - operator learning: Algorithms and analysis. In Siddhartha Mishra and Alex Townsend, editors, *Numerical Analysis Meets Machine Learning*, volume 25 of *Handbook of Numerical Analysis*, pages 419–467. Elsevier, 2024.

[21] M. Kuchta, M. Nordaas, J.C.G. Verschaeve, M. Mortensen, and K.-A. Mardal. Preconditioners for saddle point systems with trace constraints coupling 2d and 1d domains. *SIAM Journal on Scientific Computing*, 38(6):B962–B987, 2016.

[22] Miroslav Kuchta, Federica Laurino, Kent-Andre Mardal, and Paolo Zunino. Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers. *SIAM Journal on Numerical Analysis*, 59(1):558 – 582, 2021.

[23] Miroslav Kuchta, Kent-Andre Mardal, and Mikael Mortensen. Preconditioning trace coupled 3d-1d systems using fractional laplacian. *Numerical Methods for Partial Differential Equations*, 35(1):375–393, 2019.

[24] Fabien Lespagnol, Céline Grandmont, Paolo Zunino, and Miguel A. Fernández. A mixed-dimensional formulation for the simulation of slender structures immersed in an incompressible flow. *Computer Methods in Applied Mechanics and Engineering*, 432, 2024.

[25] Kim Liegeois, Romain Boman, Eric T. Phipps, Tobias A. Wiesner, and Maarten Arnst. Gmres with embedded ensemble propagation for the efficient solution of parametric linear systems in uncertainty quantification of computational models. *Computer Methods in Applied Mechanics and Engineering*, 369:113188, 2020.

[26] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. *Nature Machine Intelligence*, 3(3):218–229, 2021.
[27] David Luebke. Cuda: Scalable parallel programming for high-performance scientific computing. In *2008 5th IEEE international symposium on biomedical imaging: from nano to macro*, pages 836–838. IEEE, 2008.

[28] Kent-Andre Mardal and Ragnar Winther. Preconditioning discretizations of systems of partial differential equations. *Numerical Linear Algebra with Applications*, 18(1):1–40, 2011.

[29] Nils Margenberg, Dirk Hartmann, Christian Lessig, and Thomas Richter. A neural network multigrid solver for the navier-stokes equations. *J. Comput. Phys.*, 460:110983, 2022.

[30] Christoph Meier, Maximilian J. Grill, Wolfgang A. Wall, and Alexander Popp. Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures. *International Journal of Solids and Structures*, 154:124 – 146, 2018.

[31] Sparsh Mittal. A survey on optimized implementation of deep learning models on the nvidia jetson platform. *Journal of Systems Architecture*, 97:428–442, 2019.

[32] A. Paszke. Pytorch: An imperative style, high-performance deep learning library. *arXiv preprint arXiv:1912.01703*, 2019.

[33] Eric T. Phipps, Marta D'Elia, H. Carter Edwards, Mark Hoemmen, Jonathan Hu, and Sivasankaran Rajamanickam. Embedded ensemble propagation for improving performance, portability, and scalability of uncertainty quantification on emerging computational architectures. *SIAM Journal on Scientific Computing*, 39(2):C162–C193, 2017.

[34] L. Possenti, A. Cicchetti, R. Rosati, D. Cerroni, M. L. Costantino, T. Rancati, and P. Zunino. A mesoscale computational model for microvascular oxygen transfer. *Annals of Biomedical Engineering*, 29(12):3356–3373, 2021.

[35] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred A. Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. *arXiv*, stat.ML(1806.08734), 2019.

[36] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In *Medical image computing and computer-assisted intervention–MICCAI 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part III 18*, pages 234–241. Springer, 2015.

[37] Y. Saad. *Iterative Methods for Sparse Linear Systems*. Applied Mathematics. Society for Industrial and Applied Mathematics, 2003.

[38] Youcef Saad. A flexible inner-outer preconditioned gmres algorithm. *SIAM Journal on Scientific Computing*, 14(2):461–469, 1993.

[39] Ulrich Trottenberg, Cornelius W Oosterlee, and Anton Schuller. *Multigrid*. Elsevier, 2000.

[40] Ettore Vidotto, Timo Koch, Tobias Köppl, Rainer Helmig, and Barbara Wohlmuth. Hybrid models for simulating blood flow in microvascular networks. *Multiscale Modeling and Simulation*, 17(3):1076 – 1102, 2019.

[41] Christopher Williams, Fabian Falck, George Deligiannidis, Chris Holmes, Arnaud Doucet, and Saifuddin Syed. A unified framework for u-net design and analysis. *arXiv*, stat.ML(2305.19638), 2024.

[42] Zhi-Qin John Xu, Yaoyu Zhang, and Tao Luo. Overview frequency principle/spectral bias in deep learning. *Communications on Applied Mathematics and Computation*, 2024.
# A Auxiliary results

**Lemma A.1.** Let $(V, \|\cdot\|_V)$ and $(W, \|\cdot\|_W)$ be two Banach spaces. Let $\mathcal{GL}(V, W) \subset \mathcal{L}(V, W)$ be the set of invertible bounded linear operators from $V$ to $W$, equipped with the operator norm

$$ \|A\|_{\mathcal{L}(V,W)} := \sup_{v \in V \setminus \{0\}} \frac{\|Av\|_W}{\|v\|_V}. $$

Then,

1. $\mathcal{GL}(V, W)$ is open in the operator norm topology;

2. for every $A \in \mathcal{GL}(V, W)$ there exists a unique inverse $A^{-1} \in \mathcal{L}(W, V)$. Furthermore, $A^{-1} \in \mathcal{GL}(W, V)$;

3. the inversion map $\mathcal{I}: A \mapsto A^{-1}$ is $\mathcal{GL}(V, W) \mapsto \mathcal{GL}(W, V)$ continuous.

*Proof.* We mention that a stronger result holds, as $\mathcal{I}$ is in fact holomorphic: see, e.g., the proof in [1, Proposition 4.8]. Nonetheless, we can easily prove the statements in the Lemma as follows.

1. Given $A \in \mathcal{GL}(V, W)$, let $B \in \mathcal{L}(V, W)$ be such that $\|A - B\|_{\mathcal{L}(V,W)} < 1/\|A^{-1}\|_{\mathcal{L}(W,V)}$. We notice that the operator $C = A^{-1}(A - B) = \mathrm{Id}_V - A^{-1}B$ is an endomorphism from $V$ onto itself. Furthermore,

$$ \|C\|_{\mathcal{L}(V,V)} \le \|A^{-1}\|_{\mathcal{L}(W,V)} \cdot \|A - B\|_{\mathcal{L}(V,W)} < 1. $$

It follows that, for every $k$, the $k$th power of $C$, defined via composition, satisfies $\|C^k\|_{\mathcal{L}(V,V)} \le \|C\|_{\mathcal{L}(V,V)}^k < 1$. As a consequence, the Nuemann series $\sum_{k=0}^{+\infty} C^k$ is shown to be absolutely convergent in $\mathcal{L}(V, V)$. Indeed,

$$ \sum_{k=0}^{+\infty} \|C^k\|_{\mathcal{L}(V,V)} \le \sum_{k=0}^{+\infty} \|C\|_{\mathcal{L}(V,V)}^k = \frac{1}{1 - \|C\|_{\mathcal{L}(V,V)}}. $$

Since $\mathcal{L}(V, V)$ is a Banach space, this ensures that $\sum_{k=0}^{+\infty} C^k$ converges to some element of $\mathcal{L}(V, V)$. In particular, by construction of the Nuemann series, $\sum_{k=0}^{+\infty} C^k = (\mathrm{Id}_V - C)^{-1} = (A^{-1}B)^{-1}$. Let

$$ B^{\dagger} := (A^{-1}B)^{-1}A^{-1}. \qquad (18) $$

Clearly, $B^{\dagger} \in \mathcal{L}(W, V)$ by composition. Furthermore, we have $A^{-1}BB^{\dagger} = A^{-1}$. Since $A^{-1}$ is invertible, and thus of full rank, this shows that $BB^{\dagger} = \mathrm{Id}_W$. Similarly, composing (18) with $B$ yields

$$ B^{\dagger}B = (A^{-1}B)(A^{-1}B) = \mathrm{Id}_V. $$

It follows that $B \in \mathcal{GL}(V, W)$ and $B^{\dagger} = B^{-1}$.

2. This follows immediately from the fact that invertible operators have full rank. To see this, let $A \in \mathcal{GL}(V, W)$ and assume $B_1, B_2 \in \mathcal{L}(W, V)$ are two -possibly different- inverses of $A$. For all $w \in W$ we have

$$ B_1 w = B_2 A(B_1 w) = B_2 (A B_1 w) = B_2 w $$

since $B_1 w \in V$ and $AB_1 w \in W$. Then $B_1 \equiv B_2 =: A^{-1}$. The fact that $A^{-1} \in \mathcal{GL}(W, V)$ is obvious.

3. Let $A_n \to A$ in $\mathcal{GL}(V, W)$. Define the error $E_n = A - A_n$. There exists some $n_0$ such that $\|E_n\|_{\mathcal{L}(V,W)} < 1/\|A\|_{\mathcal{L}(V,W)}$ for all $n \ge n_0$. As for our proof of (1), we notice that this ensures the convergence of the following Neumann series

$$ \sum_{k=0}^{+\infty} (A^{-1}E_n)^k = (\mathrm{Id}_V - A^{-1}E_n)^{-1}. $$
whenever $n \ge n_0$. Since

$$ \|A^{-1} - A_n^{-1}\|_{\mathcal{L}(W,V)} \le \|\mathrm{Id}_V - A_n^{-1}A\|_{\mathcal{L}(V,V)} \cdot \|A^{-1}\|_{\mathcal{L}(W,V)}, $$

it suffices to prove that $\|\mathrm{Id}_V - A_n^{-1}A\|_{\mathcal{L}(V,V)} \to 0$ as $A_n \to A$. To this end, we note that for all $n \ge n_0$

$$ \begin{aligned} \| \mathrm{Id}_V - A_n^{-1} A \|_{\mathcal{L}(V,V)} &= \| \mathrm{Id}_V - (A - E_n)^{-1} A \|_{\mathcal{L}(V,V)} \\ &= \| \mathrm{Id}_V - (\mathrm{Id}_V - A^{-1}E_n)^{-1} A^{-1} A \|_{\mathcal{L}(V,V)} \\ &= \| \mathrm{Id}_V - (\mathrm{Id}_V - A^{-1}E_n)^{-1} \|_{\mathcal{L}(V,V)}. \end{aligned} $$

In particular, $\|\mathrm{Id}_V - A_n^{-1}A\|_{\mathcal{L}(V,V)} = \left\| \sum_{k=1}^{+\infty} (A^{-1}E_n)^k \right\|_{\mathcal{L}(V,V)}$, and thus

$$ \| \mathrm{Id}_V - A_n^{-1} A \|_{\mathcal{L}(V,V)} \le \|A^{-1}\|_{\mathcal{L}(W,V)} \sum_{k=1}^{+\infty} \|E_n\|_{\mathcal{L}(V,W)}^k = \frac{\|A^{-1}\|_{\mathcal{L}(W,V)} \|E_n\|_{\mathcal{L}(V,W)}}{1 - \|E_n\|_{\mathcal{L}(V,W)}}, $$

which is vanishing for $n \to +\infty$, as $E_n \to 0$.

□