Wiener Chaos Expansion Algorithm

• Wiener Chaos Expansion is a spectral method that represents square-integrable random variables or fields using an orthogonal basis of Gaussian coordinates, typically via Hermite or Wick polynomials.
• The algorithm truncates the infinite expansion to transform complex stochastic differential equations into deterministic coefficient systems through projection and numerical quadrature.
• Its versatility is demonstrated in finance, control, and stochastic simulations, where it accelerates computations and allows error control through chaos order, basis adaptation, and neural hybrid integrations.

The Wiener Chaos Expansion (WCE) algorithm denotes a class of spectral procedures in which a square-integrable random variable, process, or random field is represented in an orthogonal basis generated by Gaussian coordinates, typically Hermite or Wick polynomials, and then truncated to a finite chaos order for computation. Since the seminal work of Wiener, the chaos expansion has evolved to a powerful methodology for studying a broad range of stochastic differential equations, backward stochastic differential equations, stochastic partial differential equations, optimal stopping problems, and inverse problems. In the cited literature, WCE is used both as an analytic decomposition and as a numerical algorithm: after truncation, one computes deterministic coefficient functions, scalar products in $L^2$, or projected “propagator” systems, thereby replacing part of the original stochastic problem by deterministic algebra, quadrature, or embarrassingly parallel Monte Carlo averages (Gorji, 2018, Delgado-Vences et al., 27 Mar 2026).

1. Mathematical structure of the expansion

A standard formulation begins with an isonormal Gaussian process $W=\{W(h):h\in\mathfrak H\}$ on a separable Hilbert space $\mathfrak H=L^2([0,T])$, together with an orthonormal basis $\{e_i\}_{i\ge 1}$. For a multi-index $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$ with $|\alpha|=\sum_i \alpha_i<\infty$, one defines the Fourier–Hermite functional

$$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$

where $H_n$ is the probabilists’ Hermite polynomial of degree $n$. Then $\{\xi_\alpha:\alpha\in\mathcal J\}$ is an orthonormal basis of $W=\{W(h):h\in\mathfrak H\}$0, and any square-integrable process $W=\{W(h):h\in\mathfrak H\}$1 admits the full expansion

$W=\{W(h):h\in\mathfrak H\}$2

A parallel discretized Brownian formulation uses a time grid $W=\{W(h):h\in\mathfrak H\}$3 and Gaussian increments

$W=\{W(h):h\in\mathfrak H\}$4

with multivariate Hermite monomials

$W=\{W(h):h\in\mathfrak H\}$5

In that setting, any $W=\{W(h):h\in\mathfrak H\}$6 has the Wiener chaos expansion

$W=\{W(h):h\in\mathfrak H\}$7

The orthogonality relations appear in several equivalent normalizations. One recurring convention is

$W=\{W(h):h\in\mathfrak H\}$8

while normalized tensor bases satisfy

$W=\{W(h):h\in\mathfrak H\}$9

For random fields, the same structure is written pointwise in $\mathfrak H=L^2([0,T])$0: $\mathfrak H=L^2([0,T])$1 Wick powers and Wick products are embedded naturally in this hierarchy. For Gaussian $\mathfrak H=L^2([0,T])$2,

$\mathfrak H=L^2([0,T])$3

and these Wick powers lie in the corresponding chaos levels (Lelong, 2019, Shi et al., 9 Mar 2026, Delgado-Vences et al., 2024).

2. Truncation, coefficient recovery, and approximation error

Practical WCE algorithms truncate the infinite expansion. A standard finite-order set is

$\mathfrak H=L^2([0,T])$4

with cardinality $\mathfrak H=L^2([0,T])$5, which grows combinatorially in $\mathfrak H=L^2([0,T])$6 and $\mathfrak H=L^2([0,T])$7. The truncated approximation takes the form

$\mathfrak H=L^2([0,T])$8

or, in discrete Brownian settings,

$\mathfrak H=L^2([0,T])$9

Coefficient computation is a projection step. For random fields,

$\{e_i\}_{i\ge 1}$0

with the expectations approximated by Monte Carlo sampling in $\{e_i\}_{i\ge 1}$1 or by high-order quadrature in the Gaussian variables $\{e_i\}_{i\ge 1}$2. In pseudo-spectral settings one evaluates the model at quadrature nodes and forms the corresponding weighted inner products.

The $\{e_i\}_{i\ge 1}$3 tail controls truncation error. By orthonormality,

$\{e_i\}_{i\ge 1}$4

In the diffusion-bridge setting, after truncating both the maximal polynomial order $\{e_i\}_{i\ge 1}$5 and the number of basis functions $\{e_i\}_{i\ge 1}$6, one has the mean-square bound

$\{e_i\}_{i\ge 1}$7

In parameter-estimation formulations, the same tail identity is used operationally: one typically increases $\{e_i\}_{i\ge 1}$8 until the neglected variance falls below a tolerance. This suggests a generic numerical principle across the literature: WCE accuracy is controlled by a joint balance among chaos order, the number of retained Gaussian modes, and the numerical accuracy of the projection used to recover coefficients (Delgado-Vences et al., 27 Mar 2026, Delgado-Vences et al., 2024, Shi et al., 9 Mar 2026).

3. Recurrent algorithmic patterns

Although implementations differ by problem class, several algorithmic motifs recur.

A first motif is algebraic conditional expectation. In the truncated chaos space for Bermudan options,

$\{e_i\}_{i\ge 1}$9

so taking conditional expectation is equivalent to dropping all multi-indices that involve non-observable increments. This replaces least-squares regression by chaos projection in the policy iteration of a Longstaff–Schwartz-type scheme.

A second motif is forward Picard iteration with chaos projection for BSDEs. One sets

$$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$0

replaces $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$1 by $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$2, and computes

$$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$3

$$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$4

For BSDEs with jumps, the same structure is extended by adding the compensated Poisson random measure and the corresponding jump Malliavin derivative.

A third motif is the propagator system. After inserting a truncated chaos ansatz into an SDE or SPDE and projecting onto each chaos mode, one obtains a deterministic Volterra integral, ODE, or PDE hierarchy for the coefficient functions. In the parameter-estimation framework, this infinite hierarchy is called the propagator. In the stochastic generalized Kuramoto–Sivashinsky equation and stochastic Maxwell equations, the projection yields coupled deterministic PDEs for the chaos coefficients $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$5 or $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$6, which are then evolved by standard deterministic solvers.

A fourth motif is independent Monte Carlo scalar products. In the Bermudan option algorithm, each chaos coefficient $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$7 is an independent Monte Carlo average. The paper exploits this to obtain an embarrassingly parallel algorithm with two levels of parallelism, and reports that on up to 512 cores the parallel efficiency remains $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$8, with total clock time reduced from $$\alpha=(\alpha_1,\alpha_2,\dots)\in\mathcal J\subset\mathbb N_0^\infty$$9 to $|\alpha|=\sum_i \alpha_i<\infty$0 (Lelong, 2019, Briand et al., 2012, Ji et al., 2 Aug 2025).

4. Basis adaptation, dimension reduction, and transformed coordinates

A central practical difficulty of WCE is the rapid growth of the retained basis. One response is basis adaptation. In gradient-informed basis adaptation for Legendre chaos expansions, one considers a computational model $|\alpha|=\sum_i \alpha_i<\infty$1 with random input $|\alpha|=\sum_i \alpha_i<\infty$2, forms the gradient-covariance matrix

$|\alpha|=\sum_i \alpha_i<\infty$3

and computes its eigendecomposition

$|\alpha|=\sum_i \alpha_i<\infty$4

If $|\alpha|=\sum_i \alpha_i<\infty$5, then the first eigenvector $|\alpha|=\sum_i \alpha_i<\infty$6 defines a one-dimensional active subspace through

$|\alpha|=\sum_i \alpha_i<\infty$7

Because the induced marginal density $|\alpha|=\sum_i \alpha_i<\infty$8 may be nonstandard, the method introduces the auxiliary germ

$|\alpha|=\sum_i \alpha_i<\infty$9

and constructs a normalized univariate Legendre expansion

$$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$0

The reported savings are explicit. For a quadratic test with $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$1, a full 2nd-order PCE via sparse grid (level-2) uses 221 runs, whereas a 1D active-subspace Legendre expansion of order 20 uses only 33 runs and attains errors $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$2. In a multiphase-flow example with $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$3 uncertain parameters, the 1-D active-subspace Legendre surrogate yields bounded outputs, accurate PDFs, and $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$4 fewer runs overall compared to full 5-D, 4th-order expansion (801 runs) (Tsilifis, 2016).

A different reduction strategy is the logarithmic gradient transformation for Itô processes. Under smoothness and ellipticity assumptions, the Itô diffusion

$$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$5

is recast into

$$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$6

where $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$7 and $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$8 is the law density. The transformed process has the same time-marginals, but no Brownian increment: its only randomness lives in $$\xi_\alpha(\omega)=\prod_{i=1}^\infty H_{\alpha_i}\!\bigl(W(e_i)\bigr),$$9. A plausible implication is that the stochastic dimensionality is frozen at the initial condition, which is precisely the feature exploited by the subsequent Hermite-chaos discretization (Gorji, 2018).

5. Finance, control, and stochastic simulation

WCE algorithms have become particularly prominent in computational finance and stochastic simulation.

For path-dependent Bermudan options, the chaos-based policy iteration replaces the standard least-squares regression of Longstaff–Schwartz by a Wiener chaos expansion. The method can deal with a non Markovian setting, and in all tested examples $H_n$0 gives already good accuracy, within $H_n$1–$H_n$2. Numerical experiments on Bermudan put, basket put, Asian and moving-average options, and a put in the Heston model show that the chaos approach recovers Longstaff–Schwartz or benchmark prices within Monte Carlo error (Lelong, 2019).

For BSDEs, WCE is combined with Picard iterations, forward conditional-expectation formulas, and Malliavin derivatives for $H_n$3. The total mean-square error is decomposed into Picard iteration error, chaos-truncation error, basis-truncation error, and Monte Carlo error. The fully discrete estimate in one formulation is

$H_n$4

The jump-diffusion extension preserves the same architecture, replacing pure Wiener chaos by a mixed Brownian–Poisson chaos basis (Briand et al., 2012, Geiss et al., 2015).

For multi-period mean-variance portfolio allocation with transaction costs, the set of admissible portfolios is represented by a finite dimensional Wiener chaos expansion. Terminal discounted wealth is parameterized by chaos coefficients $H_n$5, the trading rule $H_n$6 is obtained by martingale projection, and the resulting finite-dimensional problem is optimized by stochastic gradient descent. The paper states that this numerical method is able to find an optimal strategy for the allocation problem subject to transaction costs (Cousin et al., 2023).

For diffusion bridges, WCE is applied to a proposal SDE with the same endpoints. The resulting method is fast and in every attempt constructs a diffusion bridge, which means there are no rejection rates. In the Ornstein–Uhlenbeck example with $H_n$7, $H_n$8, $H_n$9, $n$0, $n$1, $n$2, and $n$3, QQ-plots and Kolmogorov–Smirnov tests show no detectable difference to the exact OU bridge sampler of Bladt–Sørensen, and CPU times scale linearly in $n$4 (Delgado-Vences et al., 2024).

6. Scientific computing, neural hybrids, and generalizations

In stochastic PDEs and stochastic wave systems, WCE frequently converts the original random evolution into a deterministic coefficient system. For the stochastic generalized Kuramoto–Sivashinsky equation driven by Brownian motion forcing, the WCE ansatz leads to deterministic propagator PDEs for the coefficients $n$5. The reported absolute error is $n$6, where $n$7 is small constant and $n$8 is the time variabe, and the relative error is order $n$9 or less (Nijimbere, 2019). For stochastic Maxwell equations driven by Wiener process, the proposed algorithm can reduce the original stochastic system to the deterministic case and separate the randomness in the computation; in a one-dimensional test, CPU time is $\{\xi_\alpha:\alpha\in\mathcal J\}$0 for WCE versus $\{\xi_\alpha:\alpha\in\mathcal J\}$1 for Monte Carlo, while third and fourth moments coincide with MC up to $\{\xi_\alpha:\alpha\in\mathcal J\}$2 relative error (Ji et al., 2 Aug 2025).

Recent work integrates WCE with learned surrogates. In neural-operator formulations for singular SPDEs, the finite chaos features are injected into the backbone neural operator by feature-wise linear modulation,

$\{\xi_\alpha:\alpha\in\mathcal J\}$3

with the WCE-FiLM-NO architecture using Wick features up to order $\{\xi_\alpha:\alpha\in\mathcal J\}$4 for the dynamic $\{\xi_\alpha:\alpha\in\mathcal J\}$5 model. The cited study states that the resulting WCE-FiLM-NO shows excellent performance on $\{\xi_\alpha:\alpha\in\mathcal J\}$6, as measured by relative $\{\xi_\alpha:\alpha\in\mathcal J\}$7 loss, out-of-distribution $\{\xi_\alpha:\alpha\in\mathcal J\}$8 loss, and autocorrelation score, all without the help of renormalisation factor (Shi et al., 9 Mar 2026). A related approach approximates each deterministic propagator in a truncated chaos expansion by a neural network; deterministic and random networks are both used, with theoretical error bounds that separate the $\{\xi_\alpha:\alpha\in\mathcal J\}$9 network-approximation term from the Brownian-mode, time-basis, and chaos-order truncation tails (Neufeld et al., 2024).

WCE has also been generalized beyond the classical Hermite basis. Fractional Wiener chaos replaces integer-order Hermite polynomials by the power-normalised parabolic cylinder function

$W=\{W(h):h\in\mathfrak H\}$00

which extends the classical Hermite polynomials to non-integral order $W=\{W(h):h\in\mathfrak H\}$01, preserves martingale properties, and yields a fractional Wiener–Itô expansion in $W=\{W(h):h\in\mathfrak H\}$02 (Boguslavskaya et al., 2023). In rough-path and Gaussian-process analysis, chaos expansions of signature coordinates are written as

$W=\{W(h):h\in\mathfrak H\}$03

with kernels assembled from pairings and covariance derivatives; the combinatorial factor

$W=\{W(h):h\in\mathfrak H\}$04

makes the blow-up explicit (Ferrucci et al., 2022). For optimal prediction of Gaussian stationary-increment processes, one constructs a special orthonormal basis of the Gaussian space such that each basis element is either measurable or independent with respect to the given samples; conditioning then reduces to dropping all chaos terms whose support is not contained in the “past” index set (Alpay et al., 2014).

A persistent misconception is that WCE is a single fixed solver. The cited literature instead shows a family of algorithms: projection-based, propagator-based, pseudo-spectral, regression-free conditional expectation schemes, basis-adapted reductions, and neural hybrids. Another persistent misconception is that WCE removes all high-dimensional difficulty. The sources consistently report combinatorial growth in the retained basis, quadrature or kernel singularities in some settings, and the continued relevance of truncation order, mode selection, sparse structure, active subspaces, or learned coefficient models. This suggests that the defining feature of the Wiener Chaos Expansion algorithm is not a unique implementation, but a common spectral mechanism for converting Gaussian-driven randomness into a structured orthogonal expansion amenable to deterministic computation, controlled truncation, and application-specific acceleration (Shi et al., 9 Mar 2026, Neufeld et al., 2024, Boguslavskaya et al., 2023, Ferrucci et al., 2022, Alpay et al., 2014).