---
title: Bilevel Local Operator Learning (BiLO)
url: https://www.emergentmind.com/topics/bilevel-local-operator-learning-bilo
type: topic
---

# Bilevel Local Operator Learning (BiLO)

Bilevel Local Operator Learning (BiLO) is a neural network based method for solving inverse problems for partial differential equations (PDEs) by formulating the PDE inverse problem as a bilevel optimization problem. At the upper level, it minimizes the data loss with respect to the PDE parameters. At the lower level, it trains a neural network to locally approximate the PDE solution operator in the neighborhood of a given set of PDE parameters, which enables an accurate approximation of the descent direction for the upper level optimization problem. In the second part of the series, the same framework is extended to the Bayesian inference framework: at the lower level, a network approximates the local solution operator by minimizing the local operator loss with respect to the weights of the neural network, and at the upper level, the PDE parameters are sampled from the posterior distribution [2404.17789; 2507.17019].

## 1. Deterministic bilevel formulation

In the deterministic formulation, the spatial, or spatial-temporal, domain is denoted by $\Omega \subset \mathbb{R}^d$, the data-collocation points by $T$, and the observed data by $u_{\mathrm{obs}}(x)$ at $x \in T$. The PDE is assumed to take the form
$$
F(D^k u(x), \ldots, Du(x), u(x), \Theta)=0,
$$
where $\Theta \in \mathbb{R}^p$ are the unknown PDE parameters. BiLO introduces a neural network $u(x,\Theta;W)$ parameterized by weights $W$ to locally approximate the solution operator $u(x;\Theta)$ in a neighborhood of the current $\Theta$ [2404.17789].

The upper-level loss is the data-mismatch
$$
L_{\mathrm{data}}(\Theta,W)=\frac{1}{|T|}\sum_{x\in T}|u(x,\Theta;W)-u_{\mathrm{obs}}(x)|^2.
$$
The lower-level construction begins with the pointwise residual
$$
r(x,\Theta;W)\equiv F(D^k u(x,\Theta;W),\ldots,Du(x,\Theta;W),u(x,\Theta;W),\Theta).
$$
It then penalizes both the squared residual and its sensitivity to $\Theta$:
$$
L_{\mathrm{res}}(\Theta,W)=\frac{1}{|T_r|}\sum_{x\in T_r}|r(x,\Theta;W)|^2,
$$
$$
L_{\mathrm{rg}}(\Theta,W)=\frac{1}{|T_r|}\sum_{x\in T_r}|\partial_\Theta r(x,\Theta;W)|^2.
$$
The local operator loss is
$$
L_{\mathrm{op}}(\Theta,W)=L_{\mathrm{res}}(\Theta,W)+\alpha\cdot L_{\mathrm{rg}}(\Theta,W),
$$
where $\alpha>0$ is a small constant, for example $10^{-3}$, that controls the strength of the residual-gradient term.

The bilevel problem is
$$
\Theta^*=\arg\min_\Theta L_{\mathrm{data}}(\Theta,W^*(\Theta))), \qquad
W^*(\Theta)=\arg\min_W L_{\mathrm{op}}(\Theta,W).
$$
The stated role of the lower level is not to produce a global neural surrogate over the full parameter space, but to provide a local surrogate accurate enough to supply a descent direction for the upper-level optimization.

## 2. Joint gradient descent and the local-operator mechanism

BiLO solves the two levels simultaneously by gradient descent rather than an inner-outer loop because $W^*(\Theta)$ need only be accurate enough to give a descent direction $\partial_\Theta L_{\mathrm{data}}$. The scalar-parameter algorithm takes as input data collocation sets $T_d$ and $T_r$, an initial PDE guess $\Theta^0$, network weights $W^0$, a pretraining solution $u_0(x)$ at $\Theta^0$ as an optional input, and learning rates $lr_\Theta$ and $lr_W$ [2404.17789].

Optional pretraining minimizes
$$
L_{\mathrm{op}}(\Theta^0,W)+L_{\mathrm{pre}}(W),
$$
with
$$
L_{\mathrm{pre}}(W)=\frac{1}{|T_d|}\sum_x |u(x,\Theta^0;W)-u_0(x)|^2.
$$
Fine-tuning then applies simultaneous steps for $k=0,\ldots,K-1$:
1. compute $\partial_\Theta L_{\mathrm{data}}(\Theta^k,W^k)$ via autodiff through $u(\cdot,\Theta;W)$;
2. compute $\partial_W L_{\mathrm{op}}(\Theta^k,W^k)$;
3. update
$$
\Theta^{k+1}=\Theta^k-lr_\Theta\cdot \partial_\Theta L_{\mathrm{data}}(\Theta^k,W^k),
$$
$$
W^{k+1}=W^k-lr_W\cdot \partial_W L_{\mathrm{op}}(\Theta^k,W^k).
$$

In practice, the reported implementation uses ADAM with $lr_\Theta=lr_W\approx 10^{-3}$ and $\alpha\approx 10^{-3}$. The implementation details further specify pretraining $W$ for $10$k steps on $L_{\mathrm{op}}+L_{\mathrm{pre}}$ and fine-tuning for $10$k steps by joint GD. The explicit methodological claim is that BiLO enforces strong PDE constraints and eliminates the need to balance the residual and the data loss, which is inherent to the soft PDE constraints in many existing methods.

## 3. Network architecture and enforcement of boundary or initial conditions

For scalar $\Theta$, the reported architecture begins with an embedding layer
$$
h_0=\sigma(W_x\cdot x+R\cdot \Theta+b),
$$
where $W_x\in\mathbb{R}^{m\times d}$, $R\in\mathbb{R}^{m\times p}$ is a fixed random matrix, $b\in\mathbb{R}^m$, and $\sigma=\tanh$. The hidden portion consists of $2$ fully connected layers of width $128$ with $\tanh$ activations, and the output layer is a single scalar $N(x,\Theta;W)$ [2404.17789].

Boundary-condition enforcement is handled by output transformations. For Dirichlet data $u(\cdot)=0$ on $\partial\Omega$, the construction
$$
u(x,\Theta;W)=\psi(x)\cdot N(x,\Theta;W)+g(x)
$$
is used, where $\psi$ vanishes on $\partial\Omega$ and $g$ enforces the boundary or initial data. The implementation notes specify enforcing BC/IC via output transformations as in Section 3.1 of the paper. For unknown functions constrained by positivity, the activation is $\tanh$ for the $u$-network and softplus when enforcing $f(x)\ge 0$.

The reported collocation choices illustrate the discretization used in the experiments: Fisher–KPP uses $T_r=21\times 21$ on a grid in $x,t$ and $T_d=11$ points at $t=1$; the Poisson problem uses $T_r,T_d=101$ points in $x$; and the heat equation uses $T_r=51\times 51$ in $x,t$ and $T_d=11$ at $t=1$.

## 4. Inference of unknown functions through an auxiliary variable

When the PDE depends on an unknown function $f(x)$, BiLO parametrizes $f(x;V)$ by another network $V$, but avoids treating $V$ as an input to $u$ directly because it would be too high-dimensional. Instead, it introduces an auxiliary field
$$
z(x)=f(x;V).
$$
The residual is then defined with this auxiliary variable:
$$
a(x,z;W)\equiv F(D^k u(x,z;W),\ldots,u(x,z;W),z).
$$
The data loss becomes
$$
L_{\mathrm{data}}(V,W)=\frac{1}{|T_d|}\sum_x |u(x,f(x;V);W)-u_{\mathrm{obs}}(x)|^2,
$$
and regularization on $f$ is
$$
L_{\mathrm{reg}}(V)=\frac{1}{|T_f|}\sum_x |\nabla_x f(x;V)|^2,
$$
or any Tikhonov penalty [2404.17789].

The lower-level loss is
$$
L_{\mathrm{op}}(V,W)=L_{\mathrm{res\_fn}}(V,W)+\alpha\cdot L_{\mathrm{rg\_fn}}(V,W),
$$
with
$$
L_{\mathrm{res\_fn}}=\frac{1}{|T_r|}\sum_x |a(x,f(x;V);W)|^2,
$$
$$
L_{\mathrm{rg\_fn}}=\frac{1}{|T_r|}\sum_x |\partial_z a(x,f(x;V);W)|^2.
$$
The resulting bilevel problem is
$$
V^*=\arg\min_V\bigl[L_{\mathrm{data}}(V,W^*(V))+\lambda L_{\mathrm{reg}}(V)\bigr], \qquad
W^*(V)=\arg\min_W L_{\mathrm{op}}(V,W),
$$
and is again solved by joint gradient descent on $(V,W)$.

The heat-equation example is the principal illustration of this mechanism: BiLO is reported to infer an unknown initial condition from the final state and to recover non-smooth $f(x)$ robustly, while PINN requires extremely large or small $\beta$ to even converge.

## 5. Reported behavior on inverse problems and relation to PINNs

The central comparison in the deterministic paper is between BiLO and PINNs with soft-constraint objective $L_{\mathrm{data}}+\beta L_{\mathrm{res}}$. The stated advantage is strong PDE fidelity: unlike PINNs’ soft-constraint $L_{\mathrm{data}}+\beta L_{\mathrm{res}}$, BiLO enforces zero residual and zero residual-gradient sufficiently tightly in the lower level, eliminating delicate $\beta$-tuning. The other principal claim is robustness to sparse and noisy data [2404.17789].

The Fisher–KPP reaction–diffusion experiment infers $D,\rho$ from single-time data with $1e$–$4$ noise; the reported result is that BiLO recovers $D,\rho$ with $\ll 1\%$ error, while PINN degrades heavily unless $\beta$ is hand-tuned. In the $1$D Poisson problem with variable diffusion $D(x)$, BiLO is reported to capture the “hat” shape and kink, whereas PINN either oversmooths or oscillates depending on $\beta$. In the nonlinear ODE system that infers $a$ in $u'=\ldots + a u_1^3$, BiLO is described as stable under noise $0.05$, while PINN fails for most $\beta$.

These experiments define the framework’s practical positioning. The paper also lists suggested extensions for new PDE inverse problems: increasing network width or depth for higher-dimensional domains, combining adaptive collocation or importance sampling for $T_r$, leveraging multi-fidelity pretraining if a rough solver is available, and extending to time-dependent unknowns $f(x,t)$ by including $t$ in the $z$-embedding. These are presented as suggestions rather than established experimental conclusions.

## 6. Bayesian extension, posterior sampling, and low-rank adaptation

Part II extends BiLO to the Bayesian inference framework for PDE-constrained optimization problems. The finite-dimensional parameter is written as $m\in\Theta\subset\mathbb{R}^p$, the forward PDE is
$$
F(u(\cdot,m),m)=0 \quad \text{on } \Omega,
$$
and the local solution operator is approximated by
$$
u_\theta(x,m)=\mathrm{NN}(x,m;\theta),
$$
where $\theta\in\mathbb{R}^n$ are the network weights. The bilevel formulation is
$$
m\sim \exp(-\Phi(m)), \qquad
\theta^*(m)=\arg\min_\theta \mathcal{L}_{\rm op}(\theta;m),
$$
with lower-level loss
$$
\mathcal{L}_{\rm op}(\theta;m)=\int_\Omega |F(u_\theta(x,m),m)|^2\,dx+\omega\int_\Omega |\partial_m F(u_\theta(x,m),m)|^2\,dx,
$$
and posterior
$$
p(m\mid d)\propto \exp(-\Phi(m)),
$$
where a prior $\pi_0(m)$ and data-misfit
$$
\Phi(m)=\frac{1}{2\sigma^2}\|u_{\theta^*(m)}(\cdot,m)-d\|^2 + \text{(any regularization)}
$$
are assumed [2507.17019].

The upper-level Bayesian inference is solved using a gradient-based MCMC, for example Preconditioned Crank–Nicolson (pCN) or Hamiltonian Monte Carlo (HMC). At each proposed parameter $m$, the lower-level problem in $\theta$ is solved approximately to tolerance $\varepsilon$, after which $\Phi(m)$ and $\nabla_m\Phi(m)$ are computed. The gradient uses $\partial_m u_\theta$, obtained by automatic differentiation through the network and the PDE residual term. The combined algorithm is given as BiLO–LoRA–HMC.

To speed up repeated lower-level solves for varying $m$, Part II introduces LoRA (Low-Rank Adaptation). Let $\theta_0$ be the pre-trained network at some reference $m_0$. Updates are restricted to
$$
\theta(m)=\theta_0+UV^T,\qquad U\in\mathbb{R}^{n\times r},\;V\in\mathbb{R}^{n\times r},\;r\ll n.
$$
At each MCMC step, only $(U,V)$ are optimized to solve the lower-level problem to tolerance $\varepsilon$. The number of trainable parameters per leapfrog step is therefore reduced from $O(n)$ to $O(r\,n)$ with $r\ll n$, and the paper states that this drastically cuts both memory and computation.

## 7. Error analysis and reported uncertainty-quantification results

Part II analyzes the effect of inexact lower-level optimization on both the sampler dynamics and the posterior distribution. Let $W^*(m)$ denote the exact lower-level minimizer with $\mathcal{L}_{\rm op}(\theta^*;m)=0$, and let $W(m)$ be an inexact solution with $\mathcal{L}_{\rm op}(W(m);m)\le \varepsilon$. Under mild assumptions, specifically PDE stability and Lipschitz bounds, the paper states two principal results: a dynamic error in the hypergradient,
$$
\|g_{\rm approx}(m)-g_{\rm true}(m)\| = O(\varepsilon),
$$
and a static error in the posterior,
$$
D_{\mathrm{KL}}(\pi_* \| \bar\pi)=O(\varepsilon).
$$
The stated implication is that solving the lower-level problem to tolerance $\varepsilon$ perturbs both the MCMC dynamics and the target posterior only by $O(\varepsilon)$ [2507.17019].

The numerical experiments are organized around three PDE inverse problems. For nonlinear Poisson, $\lambda u_{xx}+k\tanh(u)=f$, BiLO recovers $k$ and $u$ with accuracy matching a fine-grid MH reference; LoRA reduces per-step cost by factors $>10$ with negligible loss in ESS or accuracy; and pCN and HMC acceptance rates are $\sim 90\%$ with stepsizes unreachable by BPINNs. For the $1$D point-process PDE, $u''(x)+\lambda\delta(x-z)-\mu u=0$, inference of $\lambda,\mu$ from point data yields $< 3\%$ MAP error, and BiLO-LoRA with $r=2\!-\!4$ achieves MOSI comparable to full fine-tuning at $r=512$. For $2$D Darcy flow, $-\nabla\!\cdot\!(D(\mathbf x)\nabla u)=f$, the framework infers $64$ KL-expansion parameters of $D(\mathbf x)$ from noisy observations of $u$; BiLO + LoRA with $r=4$ converges $\sim 4\times$ faster than full fine-tuning for width $\ge 1$k networks, with $< 7\%$ $L^2$ error and accurate $95\%$ confidence intervals.

Across these examples, the paper states that BiLO–LoRA–HMC attains high-fidelity posterior inference, large-stepsize stability unlike standard BPINNs, speedups of $5$–$30\times$, and calibrated UQ. Taken together with Part I, this places BiLO as a bilevel framework whose defining feature is local operator learning at the lower level and parameter inference, deterministic or Bayesian, at the upper level.

Source: https://www.emergentmind.com/topics/bilevel-local-operator-learning-bilo