---
title: First-Order System Least-Squares (FOSLS)
url: https://www.emergentmind.com/topics/first-order-system-least-squares-fosls
type: topic
---

# First-Order System Least-Squares (FOSLS)

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 [2409.20264]. 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 [2208.09616].

## 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 [2409.20264].

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 [1909.08598]. For time-stepping parabolic reaction–convection–diffusion, the per-step least-squares bilinear form is coercive, so each time step has a unique solution [1805.04147]. For high-frequency Helmholtz, the first-order least-squares discretization yields a Hermitian positive definite algebraic system despite the indefiniteness of the underlying PDE [1409.3362].

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 [1407.4558].

## 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 | \(\sigma=A\nabla u\) | \(\sigma-A\nabla u=0\), \(-\operatorname{div}\sigma+Bu-f=0\) [2204.07227] |
| Poisson / reaction–diffusion | \(\sigma=\nabla u\) or \(\vec w\) | \(\operatorname{div}\sigma=-f\), \(-\nabla u+\sigma=0\); weighted divergence/curl systems [2210.14063], [1909.08598] |
| Helmholtz | \(\boldsymbol\phi=ik^{-1}\nabla u\) | \(ik\boldsymbol\phi+\nabla u=0\), \(iku+\nabla\cdot\boldsymbol\phi=-ik^{-1}f\) [1409.3362] |
| Heat / parabolic | \((u_1,\mathbf u_2)=(u,-\nabla_x u)\) or \((u,-{\bf A}\nabla_x u)\) | \(\partial_t u_1+\operatorname{div}_x\mathbf u_2=\cdots\), \(-\mathbf u_2-\nabla_x u_1=\cdots\) [2208.10824], [2208.09616] |
| Instationary Stokes / Darcy–Brinkman | \(\mathbf w=-\mathbf T(\nu\mathbf u,p)\), \(M=t u-\frac{p}{t}I\) | stress or pseudostress reformulations with pressure recovery or mean-zero constraints [2201.10843], [2211.08961] |
| Nonlinear flow models | stress tensor, \(\mathbf V=\nabla \mathbf u\), \(\mathbf B=\nabla\phi\) | constitutive, momentum, incompressibility, and phase-field residuals [1809.01344], [1009.5711] |

For scalar elliptic problems, the natural trial space is typically \(H^1(\Omega)\times H(\operatorname{div};\Omega)\), or a constrained version encoding Dirichlet and Neumann data [2204.07227]. For time-harmonic Maxwell-type systems, the corresponding spaces involve \(H(\operatorname{curl})\) and \(H_0(\operatorname{curl})\) [2409.20264]. For simultaneous space-time formulations of parabolic problems, the natural space is a graph space combining \(L_2(I;H_0^1(\Omega))\), \(L_2(Q)^d\), and a divergence term in \(L_2(Q)\) [2208.10824]. For instationary Stokes, the solution space couples \(L_2(I;\mathbb H^1(\Omega))\), stress variables in \(L_2\), a space-time momentum residual, and an \(H^1\)-in-time divergence term [2201.10843].

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 \(H^2\)-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 [1911.02109].

## 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 \(H^1\)-conforming Lagrange spaces, while fluxes or stresses are approximated in \(H(\operatorname{div})\)-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 [2211.08961], [1809.01344], [1805.04147], [2412.04651].

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
\[
\tau = \min\left\{\frac12,\,(p+1)\sqrt{\frac{2}{b_0}\gamma^{-1}\ln N}\right\},
\]
so that thin boundary and corner layers are resolved while the formulation remains symmetric [1909.08598]. 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 [2208.10824], [2201.10843].

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
\[
\|F-\mathcal A v\|_{L^2(\Omega)}^2 = \sum_{K\in T}\|F-\mathcal A v\|_{L^2(K)}^2
\]
is used explicitly as a local error indicator [2409.20264]. For fully space-time parabolic PDEs, adaptive refinement on conforming simplicial meshes is driven by local indicators \(\eta(K;{\bf f},{\bf u}^\delta)=\|{\bf f}-G{\bf u}^\delta\|_{L(K)}\), and plain convergence of the adaptive loop is proved [2005.11000]. 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 [2211.08961].

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 [1711.06099].

## 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 \(J_\ast\) and \(J\) are uniformly equivalent to a parameter-dependent norm \(\|\cdot\|_t\), with constants independent of the perturbation parameter \(t\in(0,1]\). From this follow unique solvability, quasi-optimality,
\[
\|(u-u_h,M-M_h)\|_t \le C \inf_{(v_h,N_h)\in X_h}\|(u-v_h,M-N_h)\|_t,
\]
and a residual estimator that is efficient and reliable uniformly in the singular limit [2211.08961].

Analogous norm-equivalence statements underlie several space-time formulations. For instationary incompressible Stokes with slip boundary conditions, the operator \(\bar{\mathbf G}\) satisfies
\[
\|\bar{\mathbf G}(\mathbf u,\mathbf w,p)\|_{\mathscr F} \eqsim \|(\mathbf u,\mathbf w,p)\|_{\mathscr Z\times L_2(I\times\Omega)},
\]
which yields well-posedness, quasi-best approximation for any conforming finite-element triple, and an a posteriori estimator
\[
\eta_\delta=\|F-\bar{\mathbf G}(\mathbf u_\delta,\mathbf w_\delta,p_\delta)\|_{\mathscr F}
\]
that is reliable and efficient up to the continuity and inf-sup constants \(\mathfrak M,\mathfrak m\) [2201.10843]. For general second-order parabolic PDEs, the space-time first-order operator \(G\) is shown to be a linear isomorphism between the graph space \(U\) and the data space \(L\), which again makes the least-squares residual equivalent to the \(U\)-error [2208.09616].

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 \(L^2\) for the scalar variable [1805.04147]. For the space-time heat equation, Aubin–Nitsche-type estimates show that the scalar field error converges faster in \(L^2(Q)\) than in the full graph norm, and the divergence residual enjoys a higher-order conservation property [2412.04651]. For time-dependent heat and Stokes problems discretized by FOSLS, the discrete energy-law residual is shown theoretically and numerically to converge like \(\mathcal O(h^{2p})\), which is one order better in squared-error scaling than the usual \(\mathcal O(h^p)\) convergence of the field approximation [1709.00385].

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 [1407.4558].

## 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 \(\beta\), and numerical experiments report parameter robustness across perturbation levels \(10^{-6},10^{-8},10^{-10},10^{-12}\) [1909.08598]. For singularly perturbed Darcy/Brinkman flow, standard discrete spaces can lock as \(t\to 0\); the remedy is an enriched pseudostress space
\[
X_h^+ = \mathbb{P}^{k+1}_0(\mathcal T_h)^d \times \left(\mathbf{RT}^k(\mathcal T_h)+\mathbb{P}^{k+1}_*(\mathcal T_h)I\right),
\]
which removes the locking phenomenon [2211.08961].

FOSLS has also been adapted to indefinite wave propagation. For the Helmholtz equation at high wave number \(k\), the first-order least-squares formulation remains Hermitian positive definite, and the \(hp\)-version is shown to be quasi-optimal in \(L^2\) under the explicit resolution conditions that \(kh/p\) is sufficiently small and \(p\) is at least \(\mathcal O(\log k)\) [1409.3362]. 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 \(f\in H^{-1}(\Omega)\), standard lowest-order FOSLS is not directly defined because \(\operatorname{div}\tau+f\) need not lie in \(L^2\). The modified method replaces \(f\) by a bounded projector \(Q_h^\star f\in\mathbb P_0(\mathcal T)\), yielding a regularized FOSLS that remains quasi-optimal in weaker norms and restores an \(L^2\)-estimate for the scalar variable; on convex domains with \(f\in L^2(\Omega)\), the paper proves
\[
\|u-u_h\|\lesssim h^2\|f\|
\]
for the regularized scheme [2210.14063].

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 [1809.01344]. In two-phase flow based on the energetic variational approach, the Navier–Stokes/Allen–Cahn system is rewritten with gradient variables \(\mathbf V=\nabla \mathbf u\) and \(\mathbf B=\nabla\phi\); 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 [1009.5711].

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 [2208.10824], [2208.09616]. For instationary Stokes, the well-posedness theory relies on domain assumptions such as convexity or \(C^2\) 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 \(L_2(I;L_{2,0}(\Omega))\) or when no-slip conditions are used in the analogous formulation [2201.10843].

## 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 \(H^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 [1911.02109].

A later deep FOSLS method for elliptic PDEs makes the first-order structure explicit by introducing \(\sigma=A\nabla u\) and minimizing
\[
L(u,\sigma)=\|\sigma-A\nabla u\|_{L^2(\Omega)}^2+\|\operatorname{div}\sigma-Bu+f\|_{L^2(\Omega)}^2
\]
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,
\[
L_N(u,\sigma)=\frac{|\Omega|}{N}\sum_{k=1}^N \left(|\sigma(x_k)-A\nabla u(x_k)|^2+|\operatorname{div}\sigma(x_k)-Bu(x_k)+f(x_k)|^2\right).
\]
A central result is almost-sure \(\Gamma\)-convergence of the sampled loss to the continuous FOSLS functional, together with minimizer convergence under equi-coercivity and quasi-minimality hypotheses [2204.07227].

The neural extension is not limited to unconstrained elliptic equations. For the obstacle problem, a deep FOSLS variant simultaneously approximates the displacement \(u\), the flux \(\mathbf p\), and multiplier information through \(\lambda\) or \(\gamma=\operatorname{Div}\mathbf p+\lambda\). The admissible set is enforced structurally by writing
\[
u(x)=g(x)+d(x)\,a(v(x)), \qquad \gamma(x)=\operatorname{Div}\mathbf p(x)+a(\eta(x)),
\]
so that \(u\ge g\) and \(\gamma-\operatorname{Div}\mathbf p\ge 0\) hold by construction. The paper proves coercivity, local Lipschitz continuity, and almost-sure \(\Gamma\)-convergence for the Monte Carlo loss, and reports numerical experiments in dimensions \(10\) and \(20\), as well as on a non-Lipschitz slit domain [2508.19412].

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 [2409.20264].

The most recent robust deep variant in the supplied literature addresses transmission problems with discontinuous coefficients. There the continuous weighted loss
\[
\mathcal L(u,q)=\|\kappa^{-1/2}q+\kappa^{1/2}\nabla u\|_{(L^2)^d}^2 + 2(C_\kappa^P)^2\|\operatorname{Div}q-f\|_{L^2}^2
\]
is shown to be equivalent to a weighted energy norm with constants
\[
c_1=\frac18,\qquad c_2=2,\qquad \frac{c_2}{c_1}=16,
\]
independent of \(\Omega\) and the material coefficient \(\kappa\). 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 [2604.17549].

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.

Source: https://www.emergentmind.com/topics/first-order-system-least-squares-fosls