Multiscale Solution Decomposition (MSD)

• Multiscale Solution Decomposition (MSD) is a framework representing a solution as a sum of distinct components linked to different scales and features.
• It isolates easily approximable coarse parts from oscillatory or singular fine details, enabling localized basis construction and stable time integration.
• MSD is widely applied in PDE solvers, inverse problems, neural network architectures, and hardware-efficient computation for scalable and accurate approximations.

Multiscale Solution Decomposition (MSD) denotes a family of constructions in which a target solution is represented as a sum of components associated with different scales, frequencies, singular structures, subspaces, or hierarchical levels. Across the literature cited here, the decomposition may act on the solution itself, on finite element trace or subspace variables, on temporal singularities, on oscillatory phases, or on computational state variables. The shared objective is to isolate the part that is easy to approximate or solve from the part that carries oscillation, stiffness, fine-scale heterogeneity, or localized detail, thereby enabling localized basis construction, improved regularity, more stable time integration, or hardware-efficient computation (Li et al., 2021, Liu et al., 21 Sep 2025, Henning et al., 2012, Liang et al., 29 Nov 2025).

1. Canonical forms of decomposition

A recurrent MSD pattern is an additive split in which each term is tailored to a distinct regime. In the SD$^2$NN formulation for elliptic multiscale PDEs, the solution is decomposed as

$u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$

where $u_c$ is the smooth or coarse component and $u_f$ is the oscillatory or fine component (Li et al., 2021). In nonlocal-in-time problems, MSD extracts explicitly known singular terms and leaves a smoother unknown remainder; for the fractional relaxation equation,

$$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$

with the remainder $v$ satisfying a modified equation with smoother right-hand side (Liu et al., 21 Sep 2025). In the nonlinear Klein–Gordon setting, the decomposition absorbs high-frequency temporal oscillations into explicit phases,

$$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$

where $z$ is a slowly varying envelope and $r$ is a remainder term (Liang et al., 29 Nov 2025). In ill-posed inverse problems, the Multiscale Hierarchical Decomposition Method (MHDM) builds

$x_n = \sum_{j=0}^{n} u_j,$

with each increment $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$0 fitted to the current residual under regularization (Kindermann et al., 2023).

Setting	Representative decomposition	Main role
Elliptic multiscale PDEs	$u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$1	Separate smooth and oscillatory parts
Nonlocal-in-time equations	$u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$2 explicit singular series	Remove initial-time singularity
Nonlinear Klein–Gordon equation	explicit phases $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$3 envelope $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$4 remainder	Absorb high-frequency time oscillation
Ill-posed problems	$u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$5	Recover details across scales

This breadth shows that MSD is not restricted to a single numerical paradigm. In some works it is a decomposition of function space; in others it is a decomposition of singular structure, of oscillatory phase, or of the computational representation itself. A plausible implication is that “scale” in MSD is problem-dependent: frequency bands in neural PDE solvers, coarse/fine subspaces in FEM, singular orders near $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$6 in fractional dynamics, or hierarchical abstractions in control and data analysis.

2. Subspace decomposition in elliptic and multiscale finite element methods

A major classical lineage of MSD is solution-space decomposition for elliptic problems. The localized orthogonal decomposition method splits the fine finite element space as

$u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$7

where $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$8 is a detail space induced by a quasi-interpolation operator and $u(\bm{x}) = u_c(\bm{x}) + u_f(\bm{x}),$9 is its $u_c$0-orthogonal complement. Multiscale basis functions are obtained by correcting coarse nodal basis functions through fine-scale problems, and localization to patches of diameter of order $u_c$1 yields linear convergence of the $u_c$2-error with respect to the coarse mesh size without assumptions on the oscillations of the coefficients (Henning et al., 2012).

Several later developments refine this logic for rough or high-contrast media. Madureira and Sarkis’ hybrid localized spectral decomposition method uses a primal hybrid formulation and a decomposition of the discrete trace space into coarse constants and zero-average fluctuations. The resulting non-local problem has exponentially decaying solutions, which justifies localization, while local generalized eigenvalue problems enrich the space by slowly decaying modes in high-contrast settings and lead to optimal a priori error estimates that mitigate contrast effects (Madureira et al., 2017). The spectral LOD method for rough, high-contrast diffusion problems likewise builds local auxiliary spaces from fine-scale eigenproblems on each coarse element, defines a projection kernel, and then constructs an $u_c$3-orthogonal complement as the ideal multiscale space. Its error bounds are explicit and computable from algorithmic outputs, and the method is designed so that performance is similar to standard finite element methods for the homogeneous Dirichlet Poisson problem on smooth or convex domains (Brenner et al., 8 Nov 2025).

Other variants emphasize basis construction and interface handling. The multiscale discontinuous Petrov–Galerkin method combines classical oversampling multiscale basis functions with a discontinuous Petrov–Galerkin formulation, uses a homogenization-based decomposition

$u_c$4

and eliminates the resonance error term proportional to $u_c$5 that affects standard multiscale approximations (Fei et al., 2017). In discontinuous Galerkin local orthogonal decomposition for convection–diffusion problems, the trial and test spaces are spanned by corrected basis functions computed on localized patches of size $u_c$6, and the convergence rates are independent of the variation in the coefficients (Elfverson, 2015).

Residual-based variational multiscale modeling in a DG framework makes the decomposition explicit at the PDE level by writing

$u_c$7

with discontinuous coarse and fine scales. The coarse weak form then contains two distinct fine-scale contributions: a residual-based volumetric term and independent interface terms. In one dimension, classical DG methods such as interior penalty can be rederived by selecting particular fine-scale interface models, while for advection–diffusion the analysis shows that upwind numerical fluxes can be interpreted as an ad hoc remedy for missing volumetric fine-scale terms (Stoter et al., 2017).

A further branch uses hierarchical interface decomposition for massive parallelism. In recursive multiscale mixed methods, the global interface problem is replaced by a family of small interface systems on nested subdomains, implemented with the Multiscale Robin Coupled Method of Guiraldello et al. Numerical studies report excellent strong and weak scalability on very large subsurface-flow problems, including cases with billions of discretization cells (Abreu et al., 2020).

3. Singularity extraction and temporal splitting

In nonlocal-in-time problems, the salient “scale” is often the singular behavior near the initial time rather than spatial oscillation. The MSD framework for fractional relaxation equations, Volterra integral equations, subdiffusion, integrodifferential equations, and diffusion-wave equations separates a series of known singular terms from the original solution so that the remaining unknown becomes smoother. The paper proves, for the fractional relaxation equation, that the remainder $u_c$8 gains regularity relative to the original solution $u_c$9; for suitable $u_f$0, derivatives of $u_f$1 that are unbounded for $u_f$2 become bounded at $u_f$3. This regularity restoration makes smoothness assumptions used by many classical schemes appropriate for the remainder problem and improves accuracy or stability relative to classical methods (Liu et al., 21 Sep 2025).

The computational consequences are explicit. For the fractional relaxation equation with $u_f$4, numerical experiments report that classical L1 on a uniform mesh is approximately first-order, $u_f$5, whereas MSD with $u_f$6 is approximately $u_f$7, nearly optimal for $u_f$8 (Liu et al., 21 Sep 2025). For the integrodifferential equation, separating just one singular term enables the trapezoidal convolution quadrature method to achieve $u_f$9 accuracy, compared to $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$0 without MSD (Liu et al., 21 Sep 2025).

A different temporal use of MSD appears in splitting methods for nonstationary problems. Instead of splitting the operator, the framework in “Splitting methods for solution decomposition in nonstationary problems” decomposes the solution itself as $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$1 and derives coupled equations for the components in which the time derivatives of the solution components are retained explicitly. The authors construct unconditionally stable splitting schemes for a first-order evolution equation in a finite-dimensional Hilbert space and present the method as a general framework combining temporal splitting with spatial decomposition (Efendiev et al., 2020).

The same viewpoint is specialized to multiscale spatial reduction in “Temporal Splitting algorithms for non-stationary multiscale problems”. There, a GMsFEM reduced model

$$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$2

is decomposed into subspaces, producing a block system for the coefficients and a three-level implicit–explicit temporal scheme. For block-diagonal splits with $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$3 blocks, unconditional stability is obtained under the simple parameter conditions $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$4 and $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$5 (Efendiev et al., 2020). The paper’s numerical results show errors close to the full backward Euler solution together with substantial reduction in computational time.

4. Neural MSD and frequency-separated learning

In neural PDE solvers, MSD is largely motivated by spectral bias. The SD$$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$6NN architecture for elliptic multiscale PDEs combines one low-frequency normal DNN submodule with one or several high-frequency MscaleDNN submodules, and reconstructs the output as

$$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$7

where $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$8 is a scaling hyper-parameter, typically small, e.g. $$u(t) = v(t) + I_t^\alpha \sum_{i=0}^{n-1} (-\lambda I_t^\alpha)^i f(t),$$9. Orthogonality between the low- and high-frequency subnetworks is encouraged by a penalty on their inner product, and the fine submodule uses a trigonometric activation

$v$0

in the first hidden layer (Li et al., 2021).

The reported numerical behavior is explicitly decomposition-oriented. In a 1D benchmark with $v$1, SD$v$2NN2 achieves a relative error of approximately $v$3, compared with $v$4 for MscaleDNN, $v$5 for WWP, and $v$6 for a standard DNN (Li et al., 2021). FFTs of the submodule outputs confirm that each submodule captures its assigned frequency band, and the architecture is reported to remain effective for nonlinear problems, high-dimensional domains, and perforated 3D cubes (Li et al., 2021).

Neural multiscale decomposition for the nonlinear Klein–Gordon equation addresses a different oscillatory regime: strong temporal oscillation with carrier frequency $v$7. A multiscale time integrator absorbs the oscillation into explicit phase factors, reducing the learning task to a nonlinear Schrödinger equation with wave operator plus a remainder equation with small initial data (Liang et al., 29 Nov 2025). The physical approximation error between the nonlinear Klein–Gordon solution and the NLSW approximation is $v$8 in $v$9, and as $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$0 the remainder becomes negligible (Liang et al., 29 Nov 2025).

The neural implementation is two-stage: first train a network for the envelope $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$1, then a second network for the remainder $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$2, with reconstruction by the explicit phase formula. To mitigate propagation failure, the method adds a gated gradient correlation correction that enforces temporal coherence across nearby time points (Liang et al., 29 Nov 2025). Comparative experiments at $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$3 report rMAE/rRMSE values of $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$4 for PINNs, $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$5 for FLS, $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$6 for RoPINNs, and $$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$7 for NeuralMD (Liang et al., 29 Nov 2025).

5. Hierarchical decomposition for ill-posed problems

In inverse problems, MSD takes the form of iterative regularized detail recovery. The Multiscale Hierarchical Decomposition Method starts from a regularized estimate

$$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$8

then recursively solves

$$u(\mathbf{x},t) = e^{i t / \varepsilon^2} z(\mathbf{x}, t) + e^{-i t / \varepsilon^2} \overline{z}(\mathbf{x}, t) + r(\mathbf{x}, t),$$9

with $z$0 and $z$1. The resulting expansion $z$2 is interpreted as a multiscale decomposition in which coarse features are recovered first and finer details are added as $z$3 increases (Kindermann et al., 2023).

The framework extends beyond total variation denoising to linear ill-posed problems in Banach spaces, with convex and certain non-convex penalties, and it admits adaptive regularization functionals. The paper establishes convergence results for the data fidelity term and derives a multiscale identity for the data norm. It also shows that Bregman iteration can be expressed as an adaptive MHDM (Kindermann et al., 2023).

An important distinction from classical regularization is that MHDM iterates do not generally coincide with one-step variational Tikhonov regularization. Necessary and sufficient conditions are given for agreement, and one-dimensional total variation denoising is identified as a case where the iterations do agree (Kindermann et al., 2023). This corrects a common simplification: hierarchical decomposition and single-parameter Tikhonov regularization may be closely related, but they are not generically identical.

6. Generalizations beyond PDE solvers

The term MSD also appears in settings where the “solution” is not a PDE field. In multiscale Markov decision problems, the state space is partitioned into clusters linked by bottlenecks, each cluster is compressed into a coarse MDP, and the procedure is applied recursively. The hierarchy decouples subtasks, improves convergence rates locally and globally, and creates transfer opportunities for policies and potential operators across related problems (Bouvrie et al., 2012). Here the multiscale decomposition is a hierarchy of decision processes rather than a spatial basis.

For higher-order tensors, the multiscale higher-order singular value decomposition first computes a global low-rank approximation, then partitions the residual tensor into subtensors via modewise clustering, and recursively applies local HoSVDs. The decomposition

$z$4

is designed for data that are not globally low-rank but are locally low-rank after partitioning (Ozdemir et al., 2017).

Operator-adapted wavelet decomposition on polygonal mesh hierarchies carries MSD back into FEM, but with strict scale decoupling. The finest space is represented as a direct sum of $z$5-orthogonal scaling and detail spaces, each scale is solved independently, and the final solution is obtained by superposition without recomputing coarser levels. The method uses adaptive polygonal coarsening and nearly linear sparse linear-algebra operations (Şık et al., 17 Dec 2025).

At the hardware level, “Multi-Scale Dequant” uses MSD for LLM inference by decomposing BF16 activations into multiple low-precision components that can be multiplied directly with quantized weights. For INT8 weights, two-pass INT8 decomposition achieves near 16 effective bits; for MXFP4 weights, two-pass MXFP4 decomposition yields near 6.6 effective bits with error bound $z$6 per block; and the method avoids dequantization on the GEMM critical path while reducing KV cache HBM traffic by up to 2.5 times in attention (Zheng et al., 13 May 2026).

These examples indicate that MSD is best understood as a structural principle rather than a single algorithm. The decomposition target may be a PDE solution, a residual singularity, a control hierarchy, a tensor residual tree, a wavelet-scale expansion, or even an activation path in quantized inference. What remains invariant is the constructive separation of scales so that each component can be modeled, solved, transferred, or executed under assumptions that would be false or inefficient for the original undivided problem.