Hybrid Variational–MCMC Inference

• Hybrid Variational–MCMC Inference is a method that combines projection-based reduced-order modeling (e.g., POD) with on‐the-fly Krylov corrections to efficiently address large-scale, nonlinear problems.
• The approach adaptively enriches its low-rank subspace using Krylov or Newton–Krylov iterations, ensuring robust accuracy even when system dynamics evolve sharply.
• Applications span sequences of varying linear systems and nonlinear mechanics (e.g., damage or fracture) while achieving computational savings and controlled error thresholds.

Hybrid Variational–MCMC Inference integrates projection-based model reduction—most notably through Proper Orthogonal Decomposition (POD)—with Krylov subspace and Newton–Krylov iterative correction techniques to enable efficient, adaptive numerical solutions for large-scale, often highly nonlinear, problems. These approaches combine the variational construction of reduced-order models, which are efficient but may lose accuracy under strong system changes, with on-the-fly or periodically triggered Krylov subspace enrichment to restore fidelity. The result is a tunable, resource-efficient solver pipeline, applicable to sequences of linear systems with varying right-hand sides, as well as highly nonlinear mechanical systems exhibiting localized topological transformations such as crack initiation or severe damage. The central theme of these methods is their hybridization: the low-rank variational subspace provides initial efficiency, while correction and enrichment via Krylov (or Newton–Krylov) directions ensure robust accuracy even in regimes with changing or ill-represented dynamics (Carlberg et al., 2015, Kerfriden et al., 2011).

1. The Rationale for Hybrid Variational–Krylov and Variational–MCMC Inference

The key motivation for this hybridization arises from the limitations of pure reduced-order models (ROMs) when confronted with evolving right-hand sides or the emergence of localized nonlinearities and topological changes. In the classical framework, projection-based reduction (e.g., POD-Galerkin) is extremely efficient provided that all system responses are well captured by the precomputed basis. However, this premise fails in many practical applications:

• For linear systems $A_j x_j = b_j$, Krylov-subspace recycling can efficiently solve sequences when $A_j$ and $b_j$ change, but classical recycling strategies (e.g., deflation) are suboptimal for inexact solves.
• In nonlinear settings, especially those with structural damage, cracks, or local plasticity, the actual solution trajectories can deviate sharply from the span of the initial snapshots. This renders naive ROMs both inaccurate and inflexible (Carlberg et al., 2015, Kerfriden et al., 2011).

Hybrid inference schemes address this by coupling variational model reduction with adaptive Krylov or Newton–Krylov corrections, updating or enriching the basis using directions computed during the solution process itself.

2. Algorithmic Framework and Methodological Construction

Hybrid Variational–Krylov inference proceeds through two fundamental components: a variational stage that yields an efficient, low-rank subspace and an iterative Krylov-based enrichment or correction stage that restores accuracy adaptively.

POD-Augmented Krylov Recycling for Sequences of Linear Systems

For sequences of symmetric positive-definite (SPD) linear systems $A_j x_j = b_j$, the Adaptive POD–Krylov Enrichment algorithm introduces a three-stage hybrid solution scheme:

1. Direct Subspace Solution: Solve directly in the dominant POD modes $Y_1$ (dimension $r_1 \ll n$):

$(Y_1^T A_j Y_1) y_1 = Y_1^T b_j,\ \ x_1 = Y_1 y_1$

2. Reduced Augmented CG: Iteratively solve over the full POD subspace $Y_j$ (dimension $r$) to tolerance $\tau_2$:

$$(Y_j^T A_j Y_j) y_2 = Y_j^T (b_j - A_j x_1),\ \ x_2 = x_1 + Y_j y_2$$

3. Full-Space Augmented PCG: Augment with $Q = [Y_1, Y_j]$ and run augmented PCG to meet the final tolerance $\eta_j$:

$$A_j \Delta x = b_j - A_j x_2,\ \ \Delta x \in \operatorname{span}(Q) \oplus \text{Krylov}(A_j, Q^T r)$$

Snapshots (Krylov/search directions from previous solves) are accumulated and compressed into the POD basis using a goal-oriented truncation, balancing approximation error with storage cost. Noise from slow system variation is managed by selecting appropriate metrics ($\Theta = A_{j-1}$ for energy norm; $\Theta = C^T C$ for output norm) (Carlberg et al., 2015).

Adaptive POD–Krylov Enrichment for Nonlinear Problems

The algorithmic loop for nonlinear mechanics (damage, plasticity, fracture) integrates POD Galerkin models as a baseline and enriches the reduced basis adaptively:

• After each Newton iteration, both reduced ($\|V^T r\|$) and full ($\|r\|$) residuals are monitored.
• If the reduced model converges but the full residual remains high, a Newton–Krylov correction is performed on the orthogonal complement of the current reduced basis.
• This correction produces a new enrichment vector, which is orthonormalized and appended to the reduced basis. The Newton iteration continues without restart in the enlarged basis (Kerfriden et al., 2011).

A summary pseudo-code for the nonlinear variant appears in (Kerfriden et al., 2011). Hyperreduction (sampling operator $E$) can be employed to reduce the dimensionality of expensive operations within the enrichment process.

3. Mathematical Foundations

The mathematical structure of hybrid variational–MCMC inference is grounded in low-rank subspace optimization and orthogonalized iterative correction:

• Pod basis formation: Given snapshots $S \in \mathbb{R}^{n \times s}$ and weights $W = \operatorname{diag}(w_1, \ldots, w_s)$ with metric $\Theta$, perform an eigendecomposition of the weighted Gramian $G = W S^T \Theta S W$. The POD basis vectors $\Phi_r$ span the subspace minimizing

$$U_r = \arg\min_{\dim(\mathcal{A}) = r} \sum_{i=1}^s w_i^2 \| (I - P_\mathcal{A}) s_i \|_\Theta^2$$

with truncation based on an energy criterion, $$\sum_{i=1}^r \sigma_i^2 \geq (1 - \varepsilon) \sum_{i=1}^s \sigma_i^2$$ (Carlberg et al., 2015).

• Krylov enrichment: For nonlinear systems, corrections are computed in the $J$-orthogonal complement of the reduced subspace (where $J$ is the Jacobian). The correction is updated via preconditioned CG:

$P J P \Delta u_K = -P r$

with $P = I - V(V^T J V)^{-1} V^T J$ ensuring the search direction is orthogonal to $V$ in the $J$-inner product. The new direction is normalized and appended to $V$ (Kerfriden et al., 2011).

• Dimension control and truncation: To prevent unbounded basis growth, modes can be truncated adaptively and old modes discarded or periodically recompressed by SVD.

4. Computational and Numerical Performance

Hybrid inference offers substantial efficiency benefits for large-scale and sequential problems:

• For SPD linear systems, tests on problems such as contact mechanics (n ≈ $27\times 10^3$, $p=47$ systems) and torsion (n ≈ $39\times 10^3$, $p \approx 49$) revealed the following:
  • POD(100,0) (all modes in stage 1) matches untruncated methods in iteration count but with dramatically reduced orthogonalization cost.
  • Splitting between direct and iterative stages (e.g., POD(5,95)) can lower direct solve cost further while maintaining fast convergence in the final PCG stage.
  • Output-oriented POD ($\Theta = C^T C$) yields maximum convergence speed in the output norm for inexact solves.
  • The dominant cost is typically preconditioner application (when using expensive smoothers such as AMG), and recycling with hybrid schemes can drastically cut such applications (Carlberg et al., 2015).
• For nonlinear mechanics, the patch-assembly technique localizes correction costs. With parameters ($\varepsilon = 10^{-1}$, $\varepsilon_R = 10^{-6}$), the number of CG-based enrichments is modest (≤ 32), achieving errors below 1%. Cost–accuracy trade-offs can be continuously managed via the user-defined $\varepsilon$, interpolating between pure ROM and full-space Newton–Krylov (Kerfriden et al., 2011).

5. Application Domains and Extensions

The primary application domains for these hybrid inference strategies are:

• Sequences of linear systems with varying right-hand sides or system matrices (parametric or temporal evolution).
• Highly nonlinear mechanics: damage, fracture, plasticity, especially where local topological transitions occur.
• Contexts where the system Jacobian or operator varies slowly but the right-hand side experiences significant change.
• Problems where hyperreduction (sampling-restricted computations) is needed due to costly constitutive or stiffness integration.

Corrective formulations also extend efficiently to parametric studies (Monte Carlo variations in loading or material), delivering substantial computational savings compared to full-space solvers and maintaining fidelity even for sharply evolving solutions (Carlberg et al., 2015, Kerfriden et al., 2011).

6. Control Parameters, Best Practices, and Limitations

A central control parameter is the enrichment threshold $\varepsilon$, which modulates the balance between computational efficiency and solution accuracy:

• $\varepsilon \to 0$ recovers full Newton–Krylov (maximal accuracy at maximal cost).
• $\varepsilon \to \infty$ gives pure ROM (fastest but limited robustness).

Key implementation practices include:

• Setting the Krylov (CG) correction tolerance to $\varepsilon$ to avoid oversolving corrections.
• Limiting enrichment steps (e.g., 2–5 per timestep) and periodically recompressing the reduced basis by SVD to curb dimensional growth.
• Localizing stiffness reassembly (patch-assembly) in regions of large residuals for further cost reductions, particularly in damage localization scenarios.

Hyperreduced variants can exhibit large errors without correction (over 100% dissipation-error in some structural mechanics), while adaptive hybrid correction keeps errors below 5% at tolerable computational overheads (Kerfriden et al., 2011).

7. Context, Distinctions, and Future Directions

Hybrid variational–MCMC (understood here as hybrid variational–Krylov or Newton–Krylov) inference strategies distinguish themselves from classical recycling (deflation, harmonic Ritz enrichment, etc.) by:

• Focusing explicitly on inexact tolerance objectives and localized output error minimization.
• Selecting reduced basis vectors to optimize in energy- or output-norms, instead of targeting global spectra.
• Seamlessly integrating model-order reduction with on-the-fly iterative solvers, eliminating the need for restarts and supporting adaptation to evolving solution features.

A plausible implication is the potential of further extending such frameworks to multi-level or domain-decomposition schemes, in which corrections are restricted to process zones (e.g., crack tips), thereby further boosting scalability (Carlberg et al., 2015, Kerfriden et al., 2011).