High-Order Linearization Methods

• High-Order Linearization is a family of techniques that replaces complex nonlinear interactions with linear dynamics in enriched coordinate systems while preserving higher-order details.
• These methods are applied across diverse areas such as ODEs, PDEs, optimization algorithms, and kernel-based representation learning to enhance analytical and numerical performance.
• They balance improved fidelity and computational complexity by leveraging approaches like Carleman liftings, variational equations, and Fiedler pencils for structured approximations.

Searching arXiv for the most relevant papers on high-order linearization and closely related uses of the term. High-order linearization denotes a family of constructions rather than a single universally accepted method. Across the literature, it refers to procedures that replace nonlinear, higher-derivative, or higher-order interaction structure by linear or linearly organized dynamics while retaining information beyond first-order tangent approximation. In the sources considered here, the phrase covers algorithmic point-linearization of arbitrary-order ordinary differential equations, high-order local-linearization integrators, multilinear polarization schemes for higher-order equations, higher-order variational equations, higher-order corrections to Newton-type optimization, kernel linearization followed by higher-order pooling in representation learning, direct linearization of higher-derivative gravity, finite-order truncations of exact Volterra feedback linearizers for nonlinear PDEs, Carleman liftings to monomial coordinates, and commutation results between higher-order linearization and homogenization (Gerdt et al., 2015, Cruz et al., 2012, Brooks, 2023, Armstrong et al., 2019).

1. Terminological scope and recurring structures

The term “high-order” is used in several distinct senses. In some works it refers to the order of the differential equation itself, as in point-linearization of $n$-th order ODEs or Kahan-type discretizations for $n$-th order equations (Gerdt et al., 2015, Hone et al., 2019). In others it refers to the order of approximation, as in second-, third-, and fourth-order corrections to Newton, Gauss–Newton, and Levenberg–Marquardt methods (Brooks, 2023). In yet others it refers to higher-order statistics or tensor structure, as in order-$r$ pooled descriptors built after kernel linearization (Cherian et al., 2017). A further use appears in higher-derivative field theories, where “higher-order frame” means the original $f(R)$ formulation rather than an Einstein-frame reduction (Zhong et al., 2016).

A recurring distinction is between exact and approximate linearization. Exact point transformation of ODEs, exact nonlinear backstepping for hyperbolic PDEs, and exact infinite-dimensional Carleman embeddings belong to the first category, although they may be impractical or infinite-dimensional in use (Gerdt et al., 2015, Krstic, 5 Jul 2026, Itani et al., 2021). Truncated Volterra feedback laws, finite-order Carleman models, and higher-order Taylor-like optimizer updates belong to the second category: they no longer remove the nonlinearities exactly, but the uncanceled part begins at a controlled higher order (Krstic, 5 Jul 2026, Akiba et al., 7 May 2026, Brooks, 2023). This suggests that many “high-order linearization” methods are best understood as structured defect-control procedures rather than as one-shot exact reductions.

A second recurring distinction concerns what is being linearized. In some papers the object is the differential equation itself, in others the flow jet, a kernel, a nonlinear least-squares path, or a system matrix. The literature therefore does not support a single canonical definition. Instead, the common thread is the replacement of nonlinear composition laws by linear dynamics in a richer coordinate system, with the retained order determining how much nonlinear structure survives in explicit form.

2. Ordinary differential equations and dynamical systems

For scalar ODEs, one major line of work studies point-linearization of arbitrary-order equations. A nonlinear equation

$$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$

is sought so that the transformed equation takes the Laguerre normal form

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$

The key result is that the linearizability question can be reduced to the consistency of an overdetermined nonlinear PDE system in $f$, $g$, and $A_i$, together with the constraint $(A_i)_x g_y-(A_i)_y g_x=0$; differential Thomas decomposition is then used as the computational engine, and a symmetry-dimension pretest provides a fast necessary obstruction in high order (Gerdt et al., 2015). In this framework, high-order linearization is exact if the determining PDE system is consistent.

A different use of the term appears in Local Linearization–Runge–Kutta methods. There the nonlinear IVP

$n$0

is split at each step into an exactly solvable local linearized problem plus an auxiliary nonlinear remainder equation. The linear part is treated by local linearization, while the remainder is approximated by an explicit RK scheme, yielding updates of the form

$n$1

The point emphasized in the paper is that the method’s “high-order” character does not come from a higher-order Taylor linear model of the vector field; it comes from solving the remainder equation with an order-$n$2 RK method. Under $n$3, the global error is $n$4, and the resulting schemes preserve equilibria, preserve linearization at equilibria, and can be implemented as A-stable explicit integrators (Cruz et al., 2012).

In the Kahan-type literature, high-order linearization is instead algebraic multilinearization. For systems

$n$5

with $n$6, the $n$7-th derivative is replaced by the $n$8-th forward difference and each degree-$n$9 monomial is replaced by its complete symmetric multilinearization across $r$0 adjacent time levels. The construction reduces to classical Kahan discretization when $r$1, is birational, and is “order-preserving” in the specific sense that an $r$2-th order ODE becomes an $r$3-th order difference equation with the same $r$4-dimensional initial-value space (Hone et al., 2019). The paper explicitly warns that “order-preserving” here does not refer to numerical order of accuracy.

Equivalent linearization gives yet another meaning. For weakly nonlinear oscillators, the equation is rewritten as an effective linear oscillator with frequency $r$5, and resonant fundamental harmonics are eliminated iteratively. The paper’s central point is that equivalent linearization is not confined to first-order frequency corrections: when first-order averages vanish, one should average over the unperturbed period, compute the first waveform correction, substitute it back, and eliminate the fundamental again. This produces nonzero $r$6 corrections for systems such as

$r$7

and for the van der Pol oscillator, for which the paper obtains

$r$8

(Chattopadhyay et al., 2016). The misconception corrected there is that equivalent linearization is only a first-order device.

Higher-order variational equations linearize not the vector field directly but the jet dynamics of the flow along a particular solution $r$9. Writing $f(R)$0 and $f(R)$1 for the $f(R)$2-th derivatives of the flow, the nonlinear hierarchy

$f(R)$3

is embedded into a linear system on the graded symmetric algebra,

$f(R)$4

The $f(R)$5-th truncation $f(R)$6 acts on $f(R)$7, with explicit block matrices, fundamental matrices, and monodromy formulas (Simon, 2013). In this setting, high-order linearization means passing from nonlinear jet recursion to a linear system on symmetric powers.

3. Optimization and kernelized representations

In nonlinear least-squares optimization, high-order linearization is used to improve Newton, Gauss–Newton, and Levenberg–Marquardt steps in narrow curved valleys. The paper defines a “natural optimization pathway” by

$f(R)$8

so that $f(R)$9 reproduces the Newton step at $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$0, and then differentiates this path to compute $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$1, $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$2, and $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$3. The second-order term is the geodesic acceleration

$$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$4

and higher corrections generate second-, third-, and fourth-order accelerated methods. On the synthetic benchmark $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$5 with $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$6, the paper reports $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$7 iterations for first order, $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$8 for second order, $$u=f(x,y),\qquad t=g(x,y),\qquad J=f_xg_y-f_yg_x\neq 0$$9 for third order, and $u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$0 for fourth order (Brooks, 2023). The relevant caution is that this is not a classical high-order Taylor minimization of the scalar objective; it is a higher-order expansion of a curved residual-reduction path.

A different construction arises in action recognition. There, a nonlinear sequence kernel on frame-level CNN classifier scores and temporal positions is linearized into explicit feature maps by finite pivot approximations of Gaussian kernels,

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$1

and then higher-order pooling is applied to the concatenated score-time maps. The resulting HOK descriptor is

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$2

with $u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$3 used in the experiments, followed by HOSVD-based power normalization and a linear SVM (Cherian et al., 2017). The paper is explicit that the linearized object is the kernel, whereas the high-order object is the pooling statistic. That distinction is central to the method’s interpretation.

4. Partial differential equations, higher-derivative field theories, and feedback linearization

In five-dimensional warped $u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$4 gravity, linearization in the original higher-derivative formulation is complicated by the term $u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$5, which produces higher-order scalar perturbation equations. The key device is the curvature gauge

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$6

under which the dangerous higher-derivative contribution disappears from the quadratic action. After scalar-vector-tensor decomposition and elimination of constraints, the scalar normal mode takes the second-order Schrödinger-type form

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$7

while the tensor and vector sectors are already second order (Zhong et al., 2016). The companion equation-of-motion treatment derives the full linearized field equations before gauge fixing and recovers the same scalar and tensor master equations, while also identifying an additional vector constraint absent from the quadratic-action derivation (Zhong et al., 2017). The recurrent misconception addressed by these papers is that working directly in the higher-order frame necessarily leaves one with irreducibly fourth-order physical perturbation equations.

For nonlinear first-order hyperbolic PDEs, exact feedback linearization is available through an infinite Volterra backstepping operator

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$8

which transforms

$u^{(n)}(t)+\sum_{i=0}^{n-3}A_i(t)u^{(i)}(t)=0.$9

to the transport target $f$0 when the exact boundary law $f$1 is used (Krstic, 5 Jul 2026). The paper’s contribution is to truncate the controller at finite order $f$2,

$f$3

and prove that the tail

$f$4

is high-order in the state amplitude, so that every fixed truncation order still yields local stabilization. The theorem gives forward invariance, practical stability, practical finite-time attractivity after one transport crossing time, and asymptotic stability via a class-$f$5 estimate; the certified region of attraction grows with $f$6, and the practical residual shrinks with $f$7 (Krstic, 5 Jul 2026). Here high-order linearization is explicitly approximate: the first $f$8 nonlinear orders are canceled, and the residual begins at order $f$9.

5. Homogenization and infinite-dimensional analytic linearization

In nonlinear stochastic homogenization, higher-order linearization becomes a hierarchy of linearized equations around nonlinear correctors, and the central statement is that this hierarchy commutes with homogenization. For the heterogeneous equation

$g$0

the paper defines multilinear forcing fields $g$1 from derivatives of $g$2, stationary higher-order correctors $g$3, and linearized coefficient fields

$g$4

The $g$5-th linearized corrector solves

$g$6

and the effective derivatives of the homogenized Lagrangian satisfy

$g$7

The paper proves simultaneously that $g$8, that the $g$9-th heterogeneous linearized problem homogenizes to the $A_i$0-th linearized problem of $A_i$1, and that higher linearization errors satisfy optimal large-scale estimates of order $A_i$2, yielding a Liouville theorem and an explicit heterogeneous Taylor expansion of arbitrary solutions (Armstrong et al., 2019). The phrase “higher-order linearization commutes with homogenization” is therefore literal.

In infinite-dimensional Hamiltonian dynamics near an elliptic fixed point, high-order linearization is formulated as symplectic conjugacy to the quadratic Hamiltonian

$A_i$3

for systems of the form

$A_i$4

on Gevrey sequence spaces $A_i$5 (Procesi et al., 2021). The paper first defines formal symplectic linearizability by the existence of a formal generator $A_i$6 such that

$A_i$7

equivalently by the vanishing of the unique formal Birkhoff normal form $A_i$8. Under the Bourgain-type Diophantine condition

$A_i$9

and analytic majorant regularity of the perturbation, the main theorem states that formal symplectic linearization implies analytic symplectic linearization (Procesi et al., 2021). This is an infinite-dimensional analogue of the classical principle that formal nonresonant linearization converges under suitable arithmetic control.

6. Lifting methods, structure-preserving system linearizations, and computational scaling

Carleman linearization converts polynomial nonlinear dynamics into linear dynamics on monomials. For lattice Boltzmann with BGK equilibrium, the collision operator becomes quadratic in the discrete populations $(A_i)_x g_y-(A_i)_y g_x=0$0, so the lifted state $(A_i)_x g_y-(A_i)_y g_x=0$1 containing monomials up to Carleman order $(A_i)_x g_y-(A_i)_y g_x=0$2 satisfies

$(A_i)_x g_y-(A_i)_y g_x=0$3

after truncation (Itani et al., 2021). The paper gives an upper bound on the number of Carleman variables,

$(A_i)_x g_y-(A_i)_y g_x=0$4

and shows analytically that when streaming is included the truncation error scales as

$(A_i)_x g_y-(A_i)_y g_x=0$5

which underlies its statement that the error improves exponentially with Carleman order (Itani et al., 2021). At the same time, linearizing the collision step sacrifices the exactness of streaming: nonlocal products such as $(A_i)_x g_y-(A_i)_y g_x=0$6 appear, so the local lifted basis is no longer closed under transport. The paper’s D1Q3 collision-only example shows that second-order lifting is already exact for that simplified setting.

A 2026 extension addresses the practical assembly bottleneck that emerges precisely in high-order Carleman truncations. For polynomial ODEs $(A_i)_x g_y-(A_i)_y g_x=0$7, the method uses a shift-and-lift architecture with monomial bases

$(A_i)_x g_y-(A_i)_y g_x=0$8

constructs a shifted affine polynomial model around a local center $(A_i)_x g_y-(A_i)_y g_x=0$9, and then forms the truncated lifted affine dynamics

$n$00

Distinct tuples $n$01 can produce the same target exponent $n$02, so the paper introduces symmetry-reduced monomial bases, packed exponent-key indexing

$n$03

and sparse triplet coalescing during assembly (Akiba et al., 7 May 2026). Its exactness proposition states that, for fixed $n$04 and fixed shift $n$05, key-resolved triplet assembly yields the same truncated affine pair $n$06 as term-by-term expansion in exact arithmetic, up to floating-point summation order. The same paper couples shift-and-lift with moving-center expansion, rebuilding $n$07 around the current state, and compares the resulting local high-order linearization with Jacobian linearization through fixed-step error, admissible step size, and cost-at-target-accuracy criteria (Akiba et al., 7 May 2026).

For higher-order linear time-invariant state-space systems

$n$08

with regular matrix polynomial

$n$09

another form of high-order linearization is obtained by Fiedler pencils of the Rosenbrock system matrix

$n$10

The paper constructs system Fiedler pencils $n$11 that are strictly system equivalent to $n$12, hence preserve transfer function, controllability, observability, and zero structure, including input, output, and input-output decoupling zeros (Behera, 2021). Companion forms arise as special cases, and proper generalized Fiedler pencils provide block-tridiagonal and Hermitian structured linearizations. The paper also proves that zero directions of the transfer function can be recovered from eigenvectors of the Fiedler pencils without performing arithmetic operations, in the specific sense that the inverse maps reduce to block extraction in the relevant linearizations (Behera, 2021).

Taken together, these lifting and pencil-based approaches make clear that high-order linearization is often less about a single approximation formula than about choosing an enlarged coordinate system in which the nonlinear problem becomes linear, sparse, recursive, or structurally equivalent. The principal tradeoffs are equally consistent across the literature: higher retained order improves fidelity, but increases state dimension, preprocessing complexity, closure sensitivity, or the burden of solving auxiliary determining systems.