InvarianceEquation

InvarianceEquation

Assemble the cohomological linear systems that arise in the parametrisation method, a reduced-order modelling technique for high-dimensional dynamical systems.

Nomenclature

Non-resonant master modes have trivial (zero) reduced dynamics and are excluded from the cohomological equations; their columns in C(s) are omitted.

Per-multiindex cohomological equation

For each multi-index k with superharmonic s = Σᵢ kᵢ λᵢ (λᵢ eigenvalues of the master modes), the cohomological equation has the block structure

[ L(s)  C(s) ] * [ W_k; f_res ] = RHS_k

where:

• L(s) (FOM × FOM) is the parametrisation operator (characteristic-matrix polynomial of the full-order model),
• C(s) (FOM × |R|) acts on the unknown reduced-dynamics coefficients f_res of the resonant master modes,
• RHS_k contains all known lower-order contributions and external forcing.

External forcing modes are not unknowns; their contributions are handled separately and appear on the right-hand side.

Construction of L(s) and C(s) via Horner

The operator L(s) is defined as

L(s) = Σ_{k=1}^{ORD+1} B[k] s^{k-1}

where B[k] are the coefficient matrices of the linear part of the full-order model (size FOM × FOM). L(s) is evaluated efficiently using Horner's method.

The operator acting on the reduced dynamics is

C(s) = Σ_{j=1}^{ORD} D[j] s^{j-1}

with pre-computed coefficient matrices D[j] (size FOM × NVAR) given by

D[j] = Σ_{k=j+1}^{ORD+1} B[k] · generalised_right_eigenmodes
							 · reduced_dynamics_linear^{k-(j+1)}

Here:

• generalised_right_eigenmodes: NVAR × FOM matrix collecting the generalised eigenvectors of the master modes and the external forcing modes.
• reduced_dynamics_linear: Jordan matrix of the linear part of the reduced dynamics.

The D[j] matrices are pre-computed once per order using a downward recurrence (similar to the Horner scheme in MasterModeOrthogonality).

Precomputation and assembly

The coefficients D[j] are pre-computed for all NVAR = ROM + N_EXT variables. However, when assembling the linear system for a given multi-index (with superharmonic s), only a subset of the columns of C(s) is used:

• For the left-hand-side matrix [L(s) C(s)], only the columns corresponding to the resonant master modes (size |R|) are extracted from C(s). Non-resonant master modes are omitted because their reduced dynamics is identically zero.
• The external forcing modes (N_EXT columns) are handled separately and do not appear as unknowns; their contributions are moved to the right-hand side via the operator -E(s) (see below).

Right-hand-side assembly

RHS_k is the sum of two independent contributions: lower-order terms from the cohomological equation, and external forcing terms. Both are evaluated using fused Horner passes that reuse intermediate matrices to minimise computational cost.

Lower-order RHS (cohomological coupling)

During the Horner evaluation of L(s), the intermediate matrices

L[j](s) = Σ_{k=j+1}^{ORD+1} B[k] · s^{k-(j+1)},   j = 1,…,ORD

are naturally available. Multiplying each L[j](s) by a pre-computed coupling vector ξ[j] (obtained from lower-order solution coefficients) gives the contribution of lower-order terms to the RHS:

RHS_lower = -Σ_{j=1}^{ORD} L[j](s) · ξ[j]

The negative sign arises because these terms originate from the left-hand side of the cohomological equation and are moved to the right-hand side. This accumulation is performed in the same Horner loop that computes L(s), avoiding recomputation of the L[j](s) intermediates.

External forcing RHS

For external forcing modes e = 1,…,N_EXT, the polynomial coefficients E_e[L] (FOM × 1 column vectors) are pre-computed such that

E_e(s) = Σ_{L=1}^{ORD} E_e[L] · s^{L-1}

is the contribution of forcing mode e to the cohomological equation when multiplied by its known amplitude external_dynamics[e]. The total external contribution is

RHS_ext = Σ_{e=1}^{N_EXT} E_e(s) · external_dynamics[e]

To evaluate this efficiently, the coefficients of all active (non-zero) external modes are first combined into a single vector polynomial:

g_L = Σ_{e active} external_dynamics[e] · E_e[L],   L = 1,…,ORD

Then g(s) = Σ_{L=1}^{ORD} g_L · s^{L-1} is evaluated in a single Horner pass. The result is added to the RHS accumulator. This fused approach avoids evaluating each E_e(s) independently and scales only with the number of active external modes.

The complete right-hand side is therefore

RHS_k = RHS_lower + RHS_ext

where both parts are computed using dedicated fused Horner passes that share the polynomial evaluation structure of the main operator L(s).

Module contents

assemble_cohomological_matrix_and_rhs!(M, rhs, s, linear_terms, C_coeffs, E_coeffs,
										resonance, lower_order_couplings,
										external_dynamics, g_buffer) → nothing

In-place variant: writes the invariance-equation system matrix and RHS directly into the caller-supplied M and rhs buffers. No heap allocation occurs.

M must have size FOM × n_sys where n_sys = FOM + nR. rhs must have length FOM. Both are overwritten on entry. g_buffer is the pre-allocated FOM-length scratch buffer for the external RHS.

assemble_cohomological_matrix_and_rhs(
		s, linear_terms, C_coeffs, E_coeffs,
		resonance, lower_order_couplings, external_dynamics)
-> (M, rhs)

Assemble the full cohomological system matrix and the complete right-hand side for one multi-index, in a minimal number of Horner passes.

System matrix

M = [ L(s) | C_{j₁}(s)  C_{j₂}(s)  … ]

where j₁ < j₂ < … are the resonant master modes (resonance[j] == true). M has size FOM × (FOM + nR) with nR = count(resonance).

• Columns 1 : FOM → L(s) block, the parametrisation operator.
• Columns FOM+1 : FOM+nR → one column C[j](s) per resonant master mode,

						   in increasing-`j` order.

Right-hand side

rhs = lower_order_rhs + rhs_ext

where

lower_order_rhs = -Σ_{j=1}^{ORD} L[j](s) · ξ[j]
rhs_ext         = -Σ_{e=1}^{N_EXT} external_dynamics[e] · E_e(s)

The L(s) block and lower_order_rhs are computed in a single fused Horner pass by evaluate_system_matrix_and_lower_order_rhs!. rhs_ext is then accumulated into the same vector by evaluate_external_rhs!.

Arguments

• s :: T – evaluation frequency s_k = Σᵢ kᵢ λᵢ.
• linear_terms :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear_terms[k] = B[k].
• C_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed reduced-dynamics coefficients from precompute_column_polynomials; C_coeffs[j] is FOM × ORD.
• E_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed external-forcing coefficients from precompute_column_polynomials; E_coeffs[e] is FOM × ORD.
• resonance :: SVector{ROM, Bool} – resonance[j] is true iff master mode j is resonant at the current multi-index.
• lower_order_couplings :: SVector{ORD, <:AbstractVector{T}} – coupling vectors ξ[j] for j = 1,…,ORD; obtained from MORFE.LowerOrderCouplings.compute_lower_order_couplings.
• external_dynamics :: AbstractVector{T} – known amplitudes of the N_EXT external forcing modes; typically sparse.

Returns

(M, rhs) where M is FOM × (FOM + nR) and rhs is a length-FOM vector.

build_sparse_L_and_rhs!(rhs, L_template, mappings, linear_terms, s, lower_order_couplings)
-> L_template

In-place Horner evaluation of the parametrisation operator using a pre-allocated union-pattern template and index mappings from precompute_sparse_L_template.

Avoids ALL per-monomial sparse arithmetic allocations regardless of whether the linear_terms share a common sparsity pattern. The template's nzval is overwritten on each call; its colptr/rowval are never modified.

evaluate_column!(c, s, r, C_coeffs) -> c

Evaluate the r-th reduced-dynamics operator column

C_r(s) = Σ_{L=1}^{ORD} C_coeffs[r][:, L] · s^{L-1}

in-place via Horner's method, overwriting the pre-allocated FOM-vector c.

c may be a plain Vector{T} or a column view view(M, :, col).

Horner recurrence

L runs from ORD-1 down to 1:

c  ←  C_coeffs[r][:, ORD]              (highest-degree coefficient)
for L = ORD-1, …, 1:
	c ← c · s + C_coeffs[r][:, L]

Column access C_coeffs[r][:, L] reads contiguous memory (Julia is column-major), so the loop touches sequential cache lines.

Arguments

• c :: AbstractVector{T} – output buffer (length FOM), overwritten.
• s :: T – evaluation frequency.
• r :: Int – 1-based master-mode index (1 ≤ r ≤ ROM).
• C_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed coefficients from precompute_column_polynomials; C_coeffs[r] is FOM × ORD.

Complexity

O(ORD · FOM)

evaluate_external_rhs!(rhs, s, external_dynamics, E_coeffs) -> rhs

Accumulate the external-forcing contribution to the cohomological right-hand side:

rhs += Σ_{e=1}^{N_EXT} external_dynamics[e] · E_e(s)

where E_e(s) = Σ_{L=1}^{ORD} E_coeffs[e][:, L] · s^{L-1} is the e-th external column polynomial (pre-computed by precompute_column_polynomials).

The sign flip was already absorbed into the pre‑computed coefficients E_coeffs, since they originate from the left‑hand side of the invariance equation and are moved to the right‑hand side, so the contribution is added to the RHS rather than subtracted.

Sparse exploitation

Only the non-zero entries of external_dynamics are processed. For periodic forcing of a few harmonics this is typically a small subset of N_EXT.

Combined Horner pass

Rather than evaluating each E_e(s) independently, the non-zero contributions are combined into a single degree-(ORD-1) vector polynomial

g(s) = Σ_{e active} external_dynamics[e] · E_e(s)

and evaluated in one Horner pass (ORD-1 scalar-vector multiplies plus ORD-1 saxpy operations over FOM), instead of N_EXT_active separate Horner passes.

Arguments

• rhs :: AbstractVector{T} – accumulator (length FOM), updated in-place.
• s :: T – evaluation frequency.
• external_dynamics :: AbstractVector{T} – known amplitudes of the N_EXT external forcing modes; typically sparse.
• E_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed external coefficients from precompute_column_polynomials; E_coeffs[e] is FOM × ORD.

Complexity

• O(N_EXT_active · FOM · ORD) for combining coefficients.
• O(FOM · ORD) for the single Horner evaluation.

evaluate_system_matrix_and_lower_order_rhs!(parametrisation_operator,
											lower_order_rhs,
											s, lower_order_couplings,
											linear_terms)
-> parametrisation_operator

Evaluate the parametrisation operator L(s) and accumulate the lower-order right-hand-side contributions in a single Horner pass, reusing the transient intermediate matrices that are available only during the polynomial evaluation.

Mathematical context

At step j of the Horner recurrence (before the scalar multiply by s), the intermediate matrix

L[j](s) = Σ_{k=j+1}^{ORD+1} B[k] · s^{k-(j+1)}

is available. Multiplying by the pre-computed coupling vector ξ[j] = lower_order_couplings[j] gives the contribution of lower-order solution terms at derivative order j to the right-hand side:

contribution[j] = -L[j](s) · ξ[j]

The negative sign arises because these terms originate from the left-hand side of the cohomological equation and are transposed to the right-hand side.

Summed over j = 1, …, ORD, the full lower-order RHS is

lower_order_rhs = -Σ_{j=1}^{ORD} L[j](s) · ξ[j]
				= -Σ_{j=1}^{ORD} ( Σ_{k=j+1}^{ORD+1} B[k] · s^{k-(j+1)} ) · ξ[j]

This computation must share the Horner loop with L(s): the L[j] intermediates are transient, and recomputing them would double the O(ORD · FOM²) work.

The coupling vectors are obtained from MORFE.LowerOrderCouplings.compute_lower_order_couplings applied to the lower-order multi-indices associated with each Horner step.

Arguments

• parametrisation_operator :: AbstractMatrix{T} – output buffer (FOM × FOM), overwritten with L(s) = Σ_{k=1}^{ORD+1} B[k] · s^{k-1}.
• lower_order_rhs :: AbstractVector{T} – accumulator (length FOM), updated in-place. Must be initialised to zero (or the desired starting value) by the caller.
• s :: T – evaluation superharmonic.
• lower_order_couplings :: SVector{ORD, <:AbstractVector{T}} – coupling vectors ξ[j] for j = 1,…,ORD; each element is an AbstractVector{T} of length FOM.
• linear_terms :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear_terms[k] = B[k].

Complexity

O(ORD · FOM²), shared with the L(s) evaluation.

precompute_column_polynomials(fom_matrices, generalised_right_eigenmodes,
							  reduced_dynamics_linear, ROM)
-> (C_coeffs, E_coeffs)

Pre-compute the polynomial coefficient arrays for the cohomological operator columns (master modes) and the external-forcing right-hand-side columns.

These arrays are computed once per polynomial order and reused for every multi-index at that order.

Return values

• C_coeffs :: Vector{Matrix{T}} of length ROM, where C_coeffs[r] is FOM × ORD. Column j of C_coeffs[r] is the degree-(j-1) coefficient of the r-th reduced-dynamics operator column:

  C_r(s) = Σ_{j=1}^{ORD} C_coeffs[r][:, j] · s^{j-1}

• E_coeffs :: Vector{Matrix{T}} of length N_EXT = NVAR - ROM, where E_coeffs[e] is FOM × ORD. Column j of E_coeffs[e] is the degree-(j-1) coefficient of the e-th external-forcing operator column:

  E_e(s) = Σ_{j=1}^{ORD} E_coeffs[e][:, j] · s^{j-1}

Arguments

• fom_matrices :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear matrices of the full-order model; fom_matrices[k+1] corresponds to B[k] (0-indexed in the ODE).
• generalised_right_eigenmodes :: AbstractMatrix{T} of size FOM × NVAR – generalised eigenvectors; columns 1:ROM are the master modes, columns ROM+1:NVAR are the external forcing modes.
• reduced_dynamics_linear :: AbstractMatrix{T} of size NVAR × NVAR – Jordan-form matrix of the linear part of the reduced dynamics on the SSM.
• ROM :: Int – number of master modes (dimension of the reduced-order model).

Recurrence

The output matrices are filled by a single downward Horner recurrence (j runs from ORD down to 1) using one FOM × NVAR working buffer D:

D ← B[ORD+1] * generalised_right_eigenmodes                               (j = ORD)
C_coeffs[r][:, j] ← D[:, r]       for r = 1…ROM
E_coeffs[e][:, j] ← D[:, ROM+e]   for e = 1…N_EXT

D ← D * reduced_dynamics_linear + B[j+1] * generalised_right_eigenmodes  (j = ORD-1, …, 1)
C_coeffs[r][:, j] ← D[:, r]       for r = 1…ROM
E_coeffs[e][:, j] ← D[:, ROM+e]   for e = 1…N_EXT

After step j, column j of every per-target matrix holds the exact degree-(j-1) coefficient

D[:, ·] = Σ_{k=j+1}^{ORD+1} B[k] * generalised_right_eigenmodes * reduced_dynamics_linear^{k-(j+1)}

Complexity

• Time: O(ORD · FOM · NVAR) (dominated by ORD matrix–matrix products)
• Storage: O(ORD · FOM · NVAR)

precompute_external_column_polynomials(fom_matrices, external_directions,
										reduced_dynamics_linear, D_master_steps)
-> E_coeffs

Φext-dependent half of [`precomputecolumnpolynomials](@ref). Computes the external-mode column polynomialsEe(s)reusing the pre-saved master Horner intermediatesDmastersteps(from [precomputemastercolumn_polynomials`](@ref)) instead of recomputing the master-column work.

Pass external_directions = zeros(FOM, N_EXT) to obtain the partial (Φext = 0) Ecoeffs needed for the initial external-forcing solve.

Arguments

• D_master_steps — from precompute_master_column_polynomials; length ORD, each FOM × ROM. D_master_steps[j] is the master Horner buffer at step j.

precompute_master_column_polynomials(fom_matrices, master_modes, Λ_master)
-> (C_coeffs, D_master_steps)

Φext-independent half of [`precomputecolumnpolynomials](@ref). Computes only the master-mode column polynomialsCr(s)and saves the intermediate FOM×ROM Horner buffer at every step for later reuse by [precomputeexternalcolumn_polynomials`](@ref).

Because Λ is upper-triangular, the master columns of the Horner buffer D form a closed subsystem under the recurrence D ← D·Λ + B[j+1]·Y: for column r ≤ ROM, (D·Λ)[:,r] = Σ_{k≤r} D[:,k]·Λ[k,r] depends only on master columns. Therefore C_coeffs is independent of the external directions Φ_ext.

Returns

• C_coeffs :: Vector{Matrix{T}} — same layout as from precompute_column_polynomials.
• D_master_steps :: Vector{Matrix{T}} — length ORD; D_master_steps[j] is the FOM×ROM master block of the Horner buffer at the step that wrote C_coeffs[:,j]. Used by precompute_external_column_polynomials to avoid recomputing master work.

precompute_sparse_L_template(linear_terms) -> (L_template, mappings)

Pre-allocate a SparseMatrixCSC{ComplexF64} with the union sparsity pattern of all linear_terms, and compute index mappings so that each entry of linear_terms[k] can be accumulated into the correct position of the template's nzval array.

Returns (L_template, mappings) where mappings[k][i] is the index into L_template.nzval that corresponds to position i of linear_terms[k].nzval.

Used once at context construction; the template is reused across all monomials by the in-place build_sparse_L_and_rhs! overload, eliminating per-monomial sparse arithmetic allocations even when the input matrices do not share a common pattern (e.g. when the damping matrix is C = α*M + β*K and Julia's sparse addition drops entries that cancel to exactly zero).