First-Order System Least-Squares (FOSLS)

• FOSLS is a formulation that rewrites partial differential equations as first-order systems, enabling the minimization of the squared residual in an L² norm.
• The method transforms second-order problems into symmetric and coercive variational forms, which naturally provide accurate error estimates and support adaptive refinement.
• FOSLS is applied to a wide range of problems—from elliptic and parabolic equations to fluid dynamics and neural-network discretizations—offering robust, parameter-independent solutions.

First-Order System Least-Squares (FOSLS) denotes a class of formulations in which a partial differential equation is rewritten as an equivalent first-order system and the solution is characterized as the minimizer of the squared residual of that system in an $L^2$-type data norm. In an abstract setting, one writes $\mathcal A U = F$ with boundedly invertible $\mathcal A:V\to L$ and solves

$U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$

This residual functional is then equal or equivalent to the error in a natural graph or energy norm, which explains the central FOSLS features repeatedly emphasized across the literature: symmetric least-squares variational forms, quasi-best approximation from conforming subspaces, and residual-based a posteriori control (Opschoor et al., 2024). The framework has been developed for elliptic, parabolic, Helmholtz, Stokes, Darcy/Brinkman, Navier–Stokes, sea-ice, two-phase flow, obstacle, and transmission problems, and it now spans conforming finite elements, adaptive wavelet solvers, simultaneous space-time formulations, and meshless neural-network discretizations (Gantner et al., 2022).

1. Residual minimization as a first-order variational principle

The defining step in FOSLS is not merely least-squares minimization, but least-squares minimization after an equivalent first-order reformulation. For a first-order operator $\mathcal A$, the functional

$J(v)=\|F-\mathcal A v\|_L^2$

induces the Euler–Lagrange condition

$$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$

so the discrete method is a Galerkin method in the residual inner product. In the abstract framework for neural and finite-dimensional trial spaces, the bounded invertibility assumption

$$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$

implies the residual-error equivalence

$$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$

which is the basic stability mechanism behind FOSLS (Opschoor et al., 2024).

This formulation is repeatedly exploited because the least-squares bilinear form is symmetric, and in the cited finite-element settings it is coercive or positive definite in the relevant norm. For singularly perturbed reaction-diffusion, the weighted FOSLS weak form is explicitly shown to be symmetric and positive definite (Adler et al., 2019). For time-stepping parabolic reaction–convection–diffusion, the per-step least-squares bilinear form is coercive, so each time step has a unique solution (Führer et al., 2018). For high-frequency Helmholtz, the first-order least-squares discretization yields a Hermitian positive definite algebraic system despite the indefiniteness of the underlying PDE (Chen et al., 2014).

A recurrent consequence is that the least-squares functional is not only a loss or objective, but also an error measure. Standard FOSLS formulations therefore come with a “free” residual estimator, whereas related $LL^\ast$-type variants do not always retain this property. The div FOSLL* formulation for elliptic problems was introduced precisely as a related method that preserves efficiency and generality, but it requires an explicit residual estimator instead of inheriting the direct estimator structure typical of FOSLS (Cai et al., 2014).

2. Canonical first-order reformulations and function spaces

The first-order variables introduced in FOSLS are problem dependent, but the pattern is stable: second-order operators are replaced by flux, stress, gradient, pseudostress, or space-time flux variables so that only first derivatives enter the residual.

PDE class	Introduced variable(s)	Representative first-order structure
Second-order elliptic	$\mathcal A U = F$0	$\mathcal A U = F$1, $\mathcal A U = F$2 (Bersetche et al., 2022)
Poisson / reaction–diffusion	$\mathcal A U = F$3 or $\mathcal A U = F$4	$\mathcal A U = F$5, $\mathcal A U = F$6; weighted divergence/curl systems (Führer, 2022, Adler et al., 2019)
Helmholtz	$\mathcal A U = F$7	$\mathcal A U = F$8, $\mathcal A U = F$9 (Chen et al., 2014)
Heat / parabolic	$\mathcal A:V\to L$0 or $\mathcal A:V\to L$1	$\mathcal A:V\to L$2, $\mathcal A:V\to L$3 (Gantner et al., 2022, Gantner et al., 2022)
Instationary Stokes / Darcy–Brinkman	$\mathcal A:V\to L$4, $\mathcal A:V\to L$5	stress or pseudostress reformulations with pressure recovery or mean-zero constraints (Gantner et al., 2022, Führer et al., 2022)
Nonlinear flow models	stress tensor, $\mathcal A:V\to L$6, $\mathcal A:V\to L$7	constitutive, momentum, incompressibility, and phase-field residuals (Bertrand, 2018, Adler et al., 2010)

For scalar elliptic problems, the natural trial space is typically $\mathcal A:V\to L$8, or a constrained version encoding Dirichlet and Neumann data (Bersetche et al., 2022). For time-harmonic Maxwell-type systems, the corresponding spaces involve $\mathcal A:V\to L$9 and $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$0 (Opschoor et al., 2024). For simultaneous space-time formulations of parabolic problems, the natural space is a graph space combining $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$1, $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$2, and a divergence term in $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$3 (Gantner et al., 2022). For instationary Stokes, the solution space couples $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$4, stress variables in $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$5, a space-time momentum residual, and an $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$6-in-time divergence term (Gantner et al., 2022).

The importance of these choices is methodological as well as analytic. First-order reformulations move the problem into spaces where bounded invertibility or norm equivalence can be proved without demanding $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$7-regularity of the primal field. This is one reason FOSLS is repeatedly used for interface problems, high-frequency problems, and singularly perturbed regimes where direct second-order residual formulations are less natural (Cai et al., 2019).

3. Discretization architectures

The classical FOSLS discretization is conforming and mixed in the first-order variables. Across the finite-element literature, scalar or velocity unknowns are typically placed in $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$8-conforming Lagrange spaces, while fluxes or stresses are approximated in $U=\arg\min_{v\in V}\|F-\mathcal A v\|_L^2.$9-conforming Raviart–Thomas spaces. This combination appears in Darcy/Brinkman, sea-ice dynamics, Poisson-type FOSLS, parabolic time-stepping, and tensor-product space-time heat discretizations (Führer et al., 2022, Bertrand, 2018, Führer et al., 2018, Führer et al., 2024).

Several works adjust the discretization to the structure of the PDE rather than relying on a single generic mesh design. For singularly perturbed reaction-diffusion, a weighted FOSLS method is paired with a tensor-product Shishkin mesh, with a transition point

$\mathcal A$0

so that thin boundary and corner layers are resolved while the formulation remains symmetric (Adler et al., 2019). For space-time heat and Stokes problems, prismatic partitions are preferred because they support commuting or near-commuting interpolation operators and avoid continuity requirements that would be awkward on general simplicial space-time meshes (Gantner et al., 2022, Gantner et al., 2022).

Adaptive discretization is built directly into many FOSLS formulations because the residual is local and additive. In the abstract neural/finite-element framework, the cellwise decomposition

$\mathcal A$1

is used explicitly as a local error indicator (Opschoor et al., 2024). For fully space-time parabolic PDEs, adaptive refinement on conforming simplicial meshes is driven by local indicators $\mathcal A$2, and plain convergence of the adaptive loop is proved (Gantner et al., 2020). For singularly perturbed Darcy/Brinkman flow, the least-squares functional evaluated at the discrete solution is both reliable and efficient, which makes it particularly natural for adaptive refinement (Führer et al., 2022).

Beyond finite elements, FOSLS has also been coupled with adaptive wavelet Galerkin methods. In that setting, any well-posed second-order PDE is first converted into a well-posed first-order least-squares system, and the normal equations are then solved by an adaptive wavelet solver with optimal computational complexity and best-possible convergence rates (Rekatsinas et al., 2017).

4. Stability, norm equivalence, and higher-order structural properties

The theoretical core of FOSLS is the identification of a norm in which the least-squares residual is equivalent to the true error. In parameter-dependent Darcy/Brinkman flow, the central theorem states that the reduced least-squares functionals $\mathcal A$3 and $\mathcal A$4 are uniformly equivalent to a parameter-dependent norm $\mathcal A$5, with constants independent of the perturbation parameter $\mathcal A$6. From this follow unique solvability, quasi-optimality,

$\mathcal A$7

and a residual estimator that is efficient and reliable uniformly in the singular limit (Führer et al., 2022).

Analogous norm-equivalence statements underlie several space-time formulations. For instationary incompressible Stokes with slip boundary conditions, the operator $\mathcal A$8 satisfies

$\mathcal A$9

which yields well-posedness, quasi-best approximation for any conforming finite-element triple, and an a posteriori estimator

$J(v)=\|F-\mathcal A v\|_L^2$0

that is reliable and efficient up to the continuity and inf-sup constants $J(v)=\|F-\mathcal A v\|_L^2$1 (Gantner et al., 2022). For general second-order parabolic PDEs, the space-time first-order operator $J(v)=\|F-\mathcal A v\|_L^2$2 is shown to be a linear isomorphism between the graph space $J(v)=\|F-\mathcal A v\|_L^2$3 and the data space $J(v)=\|F-\mathcal A v\|_L^2$4, which again makes the least-squares residual equivalent to the $J(v)=\|F-\mathcal A v\|_L^2$5-error (Gantner et al., 2022).

The theory also extends beyond baseline quasi-optimality. For parabolic time-stepping FOSLS, a non-symmetric elliptic projection associated with the elliptic part of the problem is used to derive optimal error estimates in the natural norm and, under additional regularity assumptions, in $J(v)=\|F-\mathcal A v\|_L^2$6 for the scalar variable (Führer et al., 2018). For the space-time heat equation, Aubin–Nitsche-type estimates show that the scalar field error converges faster in $J(v)=\|F-\mathcal A v\|_L^2$7 than in the full graph norm, and the divergence residual enjoys a higher-order conservation property (Führer et al., 2024). For time-dependent heat and Stokes problems discretized by FOSLS, the discrete energy-law residual is shown theoretically and numerically to converge like $J(v)=\|F-\mathcal A v\|_L^2$8, which is one order better in squared-error scaling than the usual $J(v)=\|F-\mathcal A v\|_L^2$9 convergence of the field approximation (Adler et al., 2017).

A recurrent methodological distinction is that these “free” estimator and norm-equivalence properties belong to standard FOSLS but not automatically to all least-squares relatives. The div FOSLL* formulation retains quasi-optimal a priori estimates under the same regularity assumptions as the standard Galerkin method, but it requires an explicit residual estimator with separately proved reliability and efficiency bounds (Cai et al., 2014).

5. Specialized regimes, robustness issues, and methodological caveats

A major branch of FOSLS research concerns singular perturbations, high contrast, or high frequency, where naive discretizations can lose robustness. For singularly perturbed reaction-diffusion, weighted FOSLS is designed so that its natural norm acts as a balanced norm, meaning that thin layer components and smooth components are seen at comparable scale; continuity and coercivity are proved under an explicit bound on the gradient of the weight $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$0, and numerical experiments report parameter robustness across perturbation levels $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$1 (Adler et al., 2019). For singularly perturbed Darcy/Brinkman flow, standard discrete spaces can lock as $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$2; the remedy is an enriched pseudostress space

$$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$3

which removes the locking phenomenon (Führer et al., 2022).

FOSLS has also been adapted to indefinite wave propagation. For the Helmholtz equation at high wave number $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$4, the first-order least-squares formulation remains Hermitian positive definite, and the $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$5-version is shown to be quasi-optimal in $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$6 under the explicit resolution conditions that $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$7 is sufficiently small and $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$8 is at least $$\langle \mathcal A u,\mathcal A v\rangle_L=\langle F,\mathcal A v\rangle_L \qquad \forall v\in V,$$9 (Chen et al., 2014). This is a concrete counterexample to the common assumption that least-squares formulations are restricted to coercive or low-frequency problems.

Another specialized direction concerns singular data. For Poisson with $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$0, standard lowest-order FOSLS is not directly defined because $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$1 need not lie in $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$2. The modified method replaces $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$3 by a bounded projector $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$4, yielding a regularized FOSLS that remains quasi-optimal in weaker norms and restores an $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$5-estimate for the scalar variable; on convex domains with $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$6, the paper proves

$$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$7

for the regularized scheme (Führer, 2022).

Nonlinearity is not excluded from the FOSLS paradigm. In the viscous-plastic sea-ice model, stress is treated as an independent variable, the full residual functional includes momentum, constitutive, concentration, and thickness equations, and the nonlinear least-squares system is solved by Gauss–Newton iteration (Bertrand, 2018). In two-phase flow based on the energetic variational approach, the Navier–Stokes/Allen–Cahn system is rewritten with gradient variables $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$8 and $$c_{\mathcal A}\|v\|_V \le \|\mathcal A v\|_L \le C_{\mathcal A}\|v\|_V$$9; the resulting Newton-FOSLS discretization is used together with nested iteration, algebraic multigrid, and adaptive local refinement, while numerically retaining the global energy law of the model (Adler et al., 2010).

The cited literature also records several caveats. Space-time FOSLS norms can be stronger than the natural solution norm, and this can lead to disappointing rates for rough parabolic solutions when standard simplicial spaces are used; the move to prismatic partitions and near-commuting quasi-interpolants is presented precisely as a remedy (Gantner et al., 2022, Gantner et al., 2022). For instationary Stokes, the well-posedness theory relies on domain assumptions such as convexity or $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$0 boundary, and the authors note numerical deterioration of stability constants on an L-shaped domain, as well as loss of well-posedness when the divergence term is measured only in $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$1 or when no-slip conditions are used in the analogous formulation (Gantner et al., 2022).

6. Neural-network, meshless, and reduced-order extensions

Recent work has extended FOSLS from conforming finite-element and wavelet settings to neural-network trial spaces. An early deep least-squares study focused on second-order elliptic problems and emphasized that the FOSLS residual only requires first derivatives of the neural approximation, unlike direct second-order least-squares losses that effectively demand $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$2-type regularity. In that formulation, two neural branches approximate the scalar solution and flux simultaneously, and the discrete loss is a quadrature approximation of the FOSLS functional with boundary penalties (Cai et al., 2019).

A later deep FOSLS method for elliptic PDEs makes the first-order structure explicit by introducing $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$3 and minimizing

$$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$4

over neural-network trial spaces. Dirichlet and normal-flux boundary conditions are imposed strongly by construction rather than by penalty, and the meshless loss is approximated by Monte Carlo sampling,

$$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$5

A central result is almost-sure $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$6-convergence of the sampled loss to the continuous FOSLS functional, together with minimizer convergence under equi-coercivity and quasi-minimality hypotheses (Bersetche et al., 2022).

The neural extension is not limited to unconstrained elliptic equations. For the obstacle problem, a deep FOSLS variant simultaneously approximates the displacement $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$7, the flux $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$8, and multiplier information through $$c_{\mathcal A}\|U-v\|_V \le \|F-\mathcal A v\|_L \le C_{\mathcal A}\|U-v\|_V,$$9 or $LL^\ast$0. The admissible set is enforced structurally by writing

$LL^\ast$1

so that $LL^\ast$2 and $LL^\ast$3 hold by construction. The paper proves coercivity, local Lipschitz continuity, and almost-sure $LL^\ast$4-convergence for the Monte Carlo loss, and reports numerical experiments in dimensions $LL^\ast$5 and $LL^\ast$6, as well as on a non-Lipschitz slit domain (Acosta et al., 26 Aug 2025).

A more structural neural formulation views FOSLS as a bridge between adaptive finite elements and deep networks. In that framework, neural architectures exactly emulate finite-element spaces in De Rham-compatible complexes, the FOSLS residual is additive over localized subnetworks, and adaptive neural-network growth is guided by AFEM-style residual marking. Assuming exact minimization, the resulting adaptive neural sequences converge rate-optimally to the exact FOSLS solution (Opschoor et al., 2024).

The most recent robust deep variant in the supplied literature addresses transmission problems with discontinuous coefficients. There the continuous weighted loss

$LL^\ast$7

is shown to be equivalent to a weighted energy norm with constants

$LL^\ast$8

independent of $LL^\ast$9 and the material coefficient $\mathcal A U = F$00. The same work proves a passive variance reduction property for stochastic gradient estimates and adopts a reduced-order strategy in which a neural network defines a low-dimensional trial space, the optimal coefficients are obtained by a least-squares solve, and the space is then improved by gradient-based updates. ReQU activations are used to mitigate quasi-Gibbs oscillations near interfaces (Duque et al., 19 Apr 2026).

Taken together, these developments suggest a broad but technically coherent view of FOSLS: it is not a single discretization recipe, but a residual-based first-order reformulation principle that can be specialized to weighted balanced norms, parameter-robust graph norms, simultaneous space-time operators, explicit conservation or energy-law analyses, adaptive wavelet solvers, and meshless neural approximations, provided that the first-order operator remains well posed in the chosen function space.