“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Orthogonal polynomials from a stochastic Sturm-Liouville?

Questions

The question. Can you recover the corresponding family of orthogonal polynomials from a stochastic Sturm-Liouville equation?

Setup: what the words mean here

The classical (deterministic) Sturm-Liouville eigenvalue problem on an interval is

- (p (x) y^{'} (x))^{'} + q (x) y (x) = λ w (x) y (x),

with self-adjoint boundary conditions on the interval. Eigenfunctions are orthogonal with respect to the weight $w (x)$ , and the classical orthogonal polynomial families — Hermite, Laguerre, Jacobi, Legendre, Chebyshev — each correspond to a specific $(p, q, w)$ triple on a specific interval. So "Sturm-Liouville → orthogonal polynomials" is a one-way map already understood in the deterministic case: given the operator, the eigenfunctions are the polynomials.

The "stochastic" framing is what I'm trying to pin down. One natural interpretation: every 1D diffusion $d X_{t} = b (X_{t}) d t + σ (X_{t}) d W_{t}$ has an infinitesimal generator

L = b (x) \frac{d}{d x} + \frac{1}{2} σ^{2} (x) \frac{d ^{2}}{d x ^{2}}

that can be rewritten in Sturm-Liouville form by multiplying through by an appropriate weight. The orthogonal-polynomial family associated with the process arises as the eigenfunctions of this generator. Three canonical examples:

Ornstein-Uhlenbeck ( $b = - x$ , $σ = 1$ ) → Hermite polynomials.
Wright-Fisher / Jacobi diffusion on $[0, 1]$ → Jacobi polynomials.
CIR / squared-Bessel process on $(0, \infty)$ → Laguerre polynomials.

Under this reading the question becomes: given the SDE, can the orthogonal-polynomial family that diagonalizes its generator be derived in closed form, or at least systematically? A more aggressive reading of "stochastic Sturm-Liouville" would be one where $p, q, w$ are themselves random fields. That's a different and harder problem; for now I'm working on the diffusion-generator reading.

What I'm working on

The insight that pulled this together for me: the stationary distribution of the diffusion is exactly the weight function for the orthogonal polynomial family. So you don't need to write down the SL operator at all — simulate the SDE, look at where the samples accumulate, and that distribution IS the orthogonality weight.

Why this works in one line: the diffusion's generator $L$ , multiplied through by the stationary density $π (x)$ , is self-adjoint with weight $π$ . So $L$ is self-adjoint on $L^{2} (π)$ , and its eigenfunctions are orthogonal with respect to $π$ . Same statement, two angles: the weight in the SL form is the stationary density of the corresponding diffusion.

Concrete recipe:

Simulate the SDE for a long time. Discard a burn-in. Keep the samples.
The samples are distributed according to the stationary density $π (x)$ .
Any inner product against $π$ can be computed as a sample average — no histogram, no density estimate, just $⟨ f, g ⟩ \approx \frac{1}{N} \sum_{i} f (X_{i}) g (X_{i})$ .
Gram-Schmidt the monomials ${1, x, x^{2}, x^{3}, \dots}$ under this inner product.
The resulting polynomials are the orthogonal polynomial family for the diffusion. For Ornstein-Uhlenbeck: Hermite. For Jacobi diffusion: Jacobi. For CIR: Laguerre.

Numerical verification: OU → Hermite

Sanity-checking this end to end. Simulate $d X = - X d t + d W$ for $T = 10, 000$ with $d t = 0.01$ , keep 900k samples after burn-in, and run Gram-Schmidt:

"""
Recover Hermite polynomials from simulations of the Ornstein-Uhlenbeck process.

OU:   dX = -X dt + dW
Stationary distribution:  pi(x) = (1/sqrt(pi)) exp(-x^2),   i.e. N(0, 1/2)
This is exactly the Hermite weight (physicists' convention).

Procedure:
  1. Simulate OU for a long time.
  2. After burn-in, samples X_t are distributed according to pi.
  3. Inner products against pi are sample averages: <f,g> ~= mean(f(X_i) g(X_i))
  4. Gram-Schmidt {1, x, x^2, x^3, x^4} under this inner product.
  5. The result IS the orthogonal polynomial family — Hermite.
"""
import numpy as np

# 1. Simulate the OU process via Euler-Maruyama
rng     = np.random.default_rng(0)
T       = 10_000.0
dt      = 0.01
n_steps = int(T / dt)

dW = rng.standard_normal(n_steps) * np.sqrt(dt)
X  = np.empty(n_steps + 1)
X[0] = 0.0
for i in range(n_steps):
    X[i+1] = X[i] - X[i] * dt + dW[i]

samples = X[n_steps // 10:]   # drop burn-in

# 2. Inner products via sample averages (the simulation IS the measure)
def ip(f, g):
    return float(np.mean(f * g))

# 3. Gram-Schmidt on {1, x, x^2, x^3, x^4}
DEGREE = 4
P_vals, P_coeffs = [np.ones_like(samples)], [[1.0]]
for k in range(1, DEGREE + 1):
    cur, coeffs = samples ** k, [0.0] * k + [1.0]
    for j in range(k):
        c = ip(cur, P_vals[j]) / ip(P_vals[j], P_vals[j])
        cur = cur - c * P_vals[j]
        new = list(coeffs)
        for i in range(len(P_coeffs[j])):
            new[i] -= c * P_coeffs[j][i]
        coeffs = new
    P_vals.append(cur)
    P_coeffs.append(coeffs)

# 4. Compare to known monic physicists' Hermite
KNOWN = [
    [1.0],
    [0.0, 1.0],
    [-0.5, 0.0, 1.0],
    [0.0, -1.5, 0.0, 1.0],
    [0.75, 0.0, -3.0, 0.0, 1.0],
]
for k in range(DEGREE + 1):
    print(f"P_{k} recovered : {[round(c, 4) for c in P_coeffs[k]]}")
    print(f"H_{k} known     : {KNOWN[k]}")

OU simulation: T = 10000.0, dt = 0.01,  samples (after burn-in): 900,001
  Sample mean: +0.01197   (true: 0)
  Sample var : 0.49596    (true: 1/2)

Empirical pi(0): 0.56315   true pi(0) = 1/sqrt(pi) = 0.56419
Max |empirical - true| over histogram bins: 0.0079

  k  recovered                                       known monic Hermite
----------------------------------------------------------------------------------------------------
  0  +1.0000                                         +1.0000
  1  -0.0120  +1.0000                                +0.0000  +1.0000
  2  -0.4956  -0.0417  +1.0000                       -0.5000  +0.0000  +1.0000
  3  +0.0131  -1.4712  -0.0445  +1.0000              +0.0000  -1.5000  +0.0000  +1.0000
  4  +0.7202  -0.0258  -2.9256  +0.0231  +1.0000     +0.7500  +0.0000  -3.0000  +0.0000  +1.0000

Max |recovered coefficient - true coefficient|: 0.0744

Two things to read off. (a) The empirical stationary density matches $(1/ π) e^{- x^{2}}$ to within 0.008 over the histogram bins. (b) The polynomials produced by Gram-Schmidt match the monic physicists' Hermite polynomials, with coefficient errors growing from 0 at degree 0 to ~0.07 at degree 4 — accumulated Monte Carlo noise from roughly $1 0^{4}$ effectively independent samples (OU's autocorrelation time at $θ = 1$ is about 1 in time units, so 10⁶ correlated samples is ~10⁴ effective). The Hermite polynomials are recovered without ever writing down the SL operator. The simulation alone gives us the weight, and the weight gives us the polynomial family.

What this tells me about the original question: the answer is yes, you can recover the orthogonal polynomials. But the recovery is via the stationary distribution as an intermediate object, not directly from the SL form of the operator. The pipeline is SDE → simulate → empirical $π$ → Gram-Schmidt → polynomials. The fact that this works at all is the substance of the Karlin-McGregor / Wong classification — the stationary density of a "nice" 1D diffusion is one of the classical weight functions (Gaussian, Beta, Gamma) iff the SDE has the right form.

Open threads I want to think about: (1) for the random-coefficient reading of "stochastic SL" — where $p, q, w$ are themselves random — the polynomial family becomes random too. What's the right notion of "the orthogonal polynomials" in that case — pathwise families, ensemble average, something else? (2) Beyond the classical five — can the simulate-and-Gram-Schmidt recipe produce non-classical orthogonal polynomials from non-classical diffusions, and if so, what do they look like?

References I want to chase down

Karlin & McGregor (1957), "The differential equations of birth-and-death processes, and the Stieltjes moment problem" — the foundational connection between diffusions and orthogonal polynomials.
Wong (1964), "The construction of a class of stationary Markoff processes" — the classification of SDEs whose generators have classical-orthogonal-polynomial eigenfunctions.
Schoutens (2000), Stochastic Processes and Orthogonal Polynomials, Springer Lecture Notes in Statistics — book-length treatment.
Bakry & Émery (1985), "Diffusions hypercontractives" — curvature-dimension calculus that connects diffusion generators to functional inequalities on orthogonal polynomial expansions.
Ismail (2005), Classical and Quantum Orthogonal Polynomials in One Variable, Cambridge — encyclopedic reference for the classical theory.

Why this question is interesting to me

(placeholder — fill in)