Score-Based Diffusion Models in Infinite Dimensions: A Malliavin Calculus Perspective

TL;DR

This paper extends score-based diffusion models to infinite-dimensional Hilbert spaces using Malliavin calculus, deriving a formula for the score function applicable to SPDEs with spatially correlated noise.

Contribution

It introduces an infinite-dimensional score formula for SPDE-driven diffusion models, accommodating complex noise structures and preserving geometric properties.

Findings

Derived a closed-form score function for infinite-dimensional SPDEs.

Validated the score formula numerically for linear SPDEs in 1D and 2D.

Extended the analysis beyond finite-dimensional models using operator-theoretic methods.

Abstract

We study score-based diffusion modelling in infinite-dimensional separable Hilbert spaces through Malliavin calculus, extending the analysis of generative models beyond the finite-dimensional setting. The forward diffusion process is formulated as a linear stochastic partial differential equation (SPDE) driven by space--time coloured noise with a trace-class covariance operator, ensuring well-posedness in arbitrary spatial dimensions. Building on Malliavin calculus and an infinite-dimensional extension of the Bismut--Elworthy--Li formula, we derive a closed-form expression for the logarithmic derivative of the transition measure along Cameron--Martin directions, which serves as the natural infinite-dimensional analogue of the score function. Our operator-theoretic approach preserves the intrinsic geometry of Hilbert spaces and accommodates general trace-class operators, thereby incorporating spatially correlated noise without assuming semigroup invertibility. We validate the derived score formula numerically for several classes of linear SPDEs in both one and two spatial dimensions using spectral methods.

Tables

Table 1. Table 1: One-dimensional SPDE classes used for numerical validation of the Malliavin score formula ( 41 ). All operators are of the form A = f ​ ( − Δ ) A=f(-\Delta) , ensuring diagonalisation in the sine eigenbasis. Stability requires a k < 0 a_{k}<0 for all k ≥ 1 k\geq 1 .

SPDE	Operator $A$	Eigenvalues $a_k$	Parameters
Heat Equation	$\Delta$	$-{\left(k\pi \right)}^2$	—
Ornstein–Uhlenbeck	$\Delta -\kappa I$	$-{\left(k\pi \right)}^2-\kappa$	$\kappa =2$
Scaled Diffusion	$\nu \Delta$	$-\nu {\left(k\pi \right)}^2$	$\nu =0.1$
Fractional Laplacian	$-{\left(-\Delta \right)}^{\alpha }$	$-{\left(k\pi \right)}^{2\alpha }$	$\alpha =0.75$

Table 2. Table 2: Two-dimensional SPDE classes used for the score validation. The first two operators are second-order; the last two are fourth-order.

SPDE	Operator $A$	Eigenvalues $a_{k_1,k_2}$	Parameters
Stochastic Heat Equation	$\nu \Delta$	$-\nu \left(k_1^2+k_2^2\right)$	$\nu =1$
Ornstein–Uhlenbeck	$\nu \Delta -\kappa I$	$-\nu \left(k_1^2+k_2^2\right)-\kappa$	$\nu =1,\kappa =2$
Stochastic Biharmonic	$-\mu {\Delta }^2$	$-\mu {\left(k_1^2+k_2^2\right)}^2$	$\mu =1$
Swift–Hohenberg	$\left(r-1\right)I-2\Delta -{\Delta }^2$	$r-{\left(1-k_1^2-k_2^2\right)}^2$	$r=-0.5$

Equations

$$du(t)=Au(t)dt+Q^{1/2}d{W}_t,u(0)=u_0\in H,$$

$$⟨\eta ,\zeta {⟩}_{H_W}={\int }_0^t⟨\eta (r),\zeta (r){⟩}_{U}dr.$$

$$\Vert T{\Vert }_{\mathrm{HS}}^2=\sum _{k=1}^{\infty }\Vert T{e}_k{\Vert }_Y^2,$$

$$⟨h_1,h_2{⟩}_{{\mathcal{H}}_t}:=⟨{\gamma }^{-1/2}{h}_1,{\gamma }^{-1/2}{h}_2{⟩}_H,$$

$$W_t=\sum _{k=1}^{\infty }{f}_k{B}_t^k,$$

$${\Vert Q^{1/2}\Vert }_{\mathrm{HS}}^2=\sum _{k=1}^{\infty }{\Vert Q^{1/2}{f}_k\Vert }_H^2<\infty ,$$

$${\partial }_{t}u\left(t,x\right)=Au\left(t,x\right)+\eta \left(t,x\right),u\left(0,x\right)=u_0\left(x\right),x\in V$$

$$u(t)=S(t)u_0+{\int }_0^{t}S(t-s)Q^{1/2}d{W}_s,$$

$${\int }_0^t{\Vert S(s)Q^{1/2}\Vert }_{\mathrm{HS}}^{2}ds<\infty ,$$

$${\Vert S(s)Q^{1/2}\Vert }_{\mathrm{HS}}^2=\sum _{j=1}^{\infty }{\Vert S(s)Q^{1/2}{f}_j\Vert }_H^2$$

$$u(t)=S(t)u_0+{\int }_0^{t}S(t-s)Q^{1/2}d{W}_s.$$

$$D_{u_0}u(t)=S(t).$$

$$d{Y}_t=A{Y}_{t}dt,Y_0=I_H\Rightarrow Y_t=S(t).$$

$${\Vert Y_{t}v\Vert }_H={\Vert S(t)v\Vert }_H\le {\Vert S(t)\Vert }_{\mathcal{L}(H)}{\Vert v\Vert }_H$$

$${\gamma }_t={\int }_0^{t}S(s)Q^{1/2}{\left(Q^{1/2}\right)}^{\ast }S(s{)}^{\ast }ds,$$

$${\gamma }_t=Y_t{C}_t{Y}_t^{\ast },C_t={\int }_0^t{Y}_r^{-1}{Q}^{1/2}{\left(Q^{1/2}\right)}^{\ast }{\left(Y_r^{-1}\right)}^{\ast }dr,$$

$$u(t)=S(t)u_0+{\int }_0^{t}S(t-s)Q^{1/2}d{W}_s,$$

$${\int }_0^t{\Vert S(s)Q^{1/2}\Vert }_{\mathrm{HS}}^{2}ds<\infty ,$$

$${\mathcal{D}}_r\left[S(t)u_0\right]=0.$$

$${\mathcal{D}}_{r}u(t)=S(t-r)Q^{1/2}{1}_{[0,t]}(r),$$

$${\gamma }_t={\int }_0^t{\mathcal{D}}_{r}u(t){\left({\mathcal{D}}_{r}u(t)\right)}^{\ast }dr,$$

$${\gamma }_t={\int }_0^{t}S(t-r)Q^{1/2}{\left(S(t-r)Q^{1/2}\right)}^{\ast }dr.$$

$${\left(S(t-r)Q^{1/2}\right)}^{\ast }={\left(Q^{1/2}\right)}^{\ast }S(t-r{)}^{\ast }.$$

$${\mathcal{D}}_{r}u(t)({\mathcal{D}}_{r}u(t){)}^{\ast }=S(t-r)Q^{1/2}{\left(Q^{1/2}\right)}^{\ast }S(t-r{)}^{\ast },$$

$${\gamma }_t={\int }_0^{t}S(t-r)Q^{1/2}{\left(Q^{1/2}\right)}^{\ast }S(t-r{)}^{\ast }dr.$$

$${\gamma }_t={\int }_0^{t}S(s)Q^{1/2}{\left(Q^{1/2}\right)}^{\ast }S(s{)}^{\ast }ds.$$

$${\gamma }_t={\int }_0^{t}S(t-r)Q^{1/2}{\left(Q^{1/2}\right)}^{\ast }S(t-r{)}^{\ast }dr=Y_t{C}_t{Y}_t^{\ast },$$

$$C_t={\int }_0^t{Y}_r^{-1}{Q}^{1/2}{\left(Q^{1/2}\right)}^{\ast }{\left(Y_r^{-1}\right)}^{\ast }dr.$$

$$F=\left(⟨u(t),h_1{⟩}_H,\dots ,⟨u(t),h_m{⟩}_H\right)\in {\mathbb{R}}^m.$$

$${\gamma }_F={\left(⟨h_i,{\gamma }_t{h}_j{⟩}_H\right)}_{i,j=1}^m,$$

Full text

††volume: Volume x, Issue x, 2026

1 Introduction

Recent advances in diffusion generative models (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2021) have substantially advanced data synthesis. These models rely on the so-called score function (Hyvärinen, 2005, 2007; Song et al., 2019; Ni, 2025), the gradient of the log-density, to generate an iterative denoising process that reverses a forward diffusion process, achieving state-of-the-art performance in many practical applications such as image and audio generation, and design of new molecular structures. Classical score-based diffusion models are defined in finite, yet often very high-dimensional, spaces. Extending these models to infinite dimensions, such as those arising in functional data analysis, presents significant mathematical challenges. For instance, infinite-dimensional spaces lack a Lebesgue measure, making traditional density gradients ill-defined. Instead, one must work with logarithmic derivatives (or Fomin derivatives) of measures along appropriate directions, which are well-defined only in the Cameron–Martin space associated with the underlying Gaussian measure. Moreover, for diffusion processes governed by stochastic partial differential equations (SPDEs), the noise must often be regularised to guarantee well-posedness of the initial-boundary value problem. In the case of abstract SPDEs of the form

[TABLE]

where $H$ is a separable Hilbert space, $A$ is a densely defined unbounded linear operator (e.g., the Laplacian with Dirichlet or Neumann boundary conditions), and $W_{t}$ is a cylindrical Wiener process, this regularisation is achieved by taking $Q^{1/2}:U\to H$ to be Hilbert-Schmidt111More precisely, $W_{t}$ is a cylindrical Wiener process on a separable Hilbert space $U$ (hence not $U$-valued), and $Q^{1/2}:U\to H$ is a Hilbert-Schmidt operator. The induced covariance operator $Q:=Q^{1/2}(Q^{1/2})^{*}\in\mathcal{L}_{1}(H)$ is then trace-class on $H$. (Da Prato and Zabczyk, 2014), so that $Q^{1/2}dW_{t}$ defines spatially correlated (coloured) noise.

The main objective of this paper is to study score-based diffusion modelling in infinite dimensions using Malliavin calculus (Malliavin, 1978). Our analysis extends the recent work (Mirafzali et al., 2025a); (Mirafzali et al., 2025b) from finite-dimensional diffusion generative models – where it has proved effective in deriving exact closed-form expressions for the score function – to the infinite-dimensional setting of linear SPDEs with additive Gaussian noise. Our work is closely related to that of (Greco, 2025), who studied a similar problem in infinite-dimensional score-based diffusion generative models using Malliavin calculus. In their formulation, the score function is identified (via time reversal) as a Malliavin derivative and corresponds to a conditional expectation, and the approach uses Gaussian Ornstein–Uhlenbeck noising defined via Dirichlet forms and is built around $\Gamma$-calculus. Likewise, (Pidstrigach et al., 2025) employ Malliavin calculus for conditioned diffusions to derive a Tweedie-like formula for the score function, although their focus remains on finite-dimensional conditioning mechanisms rather than on score computation in infinite dimensions. On the other hand, our approach relies on an operator-theoretic derivation of the logarithmic derivative (directional score) for linear SPDEs of the form (1) with a general Hilbert-Schmidt operator $Q^{1/2}:U\to H$, without resorting to any spatial discretisation. By leveraging techniques for differentiating Hilbert-valued processes (Nualart and Nualart, 2018) together with infinite-dimensional extensions of the Bismut–Elworthy–Li formula (Bakhtin and Mattingly, 2007; Elworthy and Li, 1994), we obtain closed-form expressions for the logarithmic derivative in the infinite-dimensional setting. Other recent contributions formulate diffusion models directly on infinite-dimensional function spaces—either by perturbing functions with a Gaussian process specified by a covariance kernel and deriving a resolution–independent discretised algorithm (Lim et al., 2025), by casting the dynamics as stochastic evolution equations in Hilbert spaces and implementing them via spatial discretisation (Lim et al., 2023), or by developing an infinite-dimensional formulation with dimension–independent guarantees (Pidstrigach et al., 2024). None of these approaches employ Malliavin calculus.

This paper is organised as follows. In Section 2 we establish notation and define the relevant spaces and topologies. In Section 3 we present our methodology, including the operator-theoretic derivation of the logarithmic derivative for linear SPDEs with space–time coloured noise, the Malliavin covariance operator, and the infinite-dimensional Bismut–Elworthy–Li formula. In Section 4 we verify the score formula numerically in both one and two spatial dimensions. In Section 5, we discuss our main findings and outline possible avenues for future work. For the reader’s convenience, we provide A where we review the relevant mathematical background, including infinite-dimensional diffusion processes generated by linear SPDEs and Malliavin calculus.

2 Notation and Preliminaries

We begin by establishing notation and defining the relevant Hilbert spaces, operator spaces, and topologies used throughout this paper. This section is essential for ensuring mathematical consistency.

2.1 Hilbert Spaces and Inner Products

Throughout this paper, we work with the following separable Hilbert spaces:

• •

$H$: the state space, typically $L^{2}(V)$ for some spatial domain $V\subset\mathbb{R}^{d}$, equipped with inner product $\langle\cdot,\cdot\rangle_{H}$ and norm $\|u\|_{H}=\sqrt{\langle u,u\rangle_{H}}$;

• •

$U$: an auxiliary separable Hilbert space on which the cylindrical Wiener process is defined, equipped with inner product $\langle\cdot,\cdot\rangle_{U}$ and norm $\|\cdot\|_{U}$;

• •

$H_{W}:=L^{2}([0,t];U)$: the space of square-integrable $U$-valued functions on $[0,t]$, equipped with the inner product

[TABLE]

All closures of subsets of these Hilbert spaces are taken with respect to the corresponding Hilbert space norm unless otherwise specified.

2.2 Operator Spaces

We denote by:

• •

$\mathcal{L}(X,Y)$: the space of bounded linear operators from Hilbert space $X$ to Hilbert space $Y$, with operator norm $\|\cdot\|_{\mathcal{L}(X,Y)}$;

• •

$\mathcal{L}_{\mathrm{HS}}(X,Y)$: the space of Hilbert-Schmidt operators from $X$ to $Y$, with norm

[TABLE]

where $\{e_{k}\}_{k=1}^{\infty}$ is any orthonormal basis of $X$ (the sum is independent of the choice of basis);

• •

$\mathcal{L}_{1}(H)$: the space of trace-class operators on $H$, consisting of operators $T\in\mathcal{L}(H)$ for which $\operatorname{Tr}(|T|)<\infty$.

For a bounded operator $T\in\mathcal{L}(X,Y)$, we denote its adjoint by $T^{*}\in\mathcal{L}(Y,X)$.

2.3 Derivatives

To avoid confusion, we use distinct notation for different types of derivatives:

• •

$\nabla\phi$: the Fréchet derivative of a functional $\phi:H\to\mathbb{R}$, identified via the Riesz representation theorem with an element of $H$;

• •

$\mathcal{D}_{r}F$: the Malliavin derivative of a random variable $F$ at time $r\in[0,t]$, which takes values in $U$ for scalar-valued $F$, or in $\mathcal{L}_{\mathrm{HS}}(U,H)$ for $H$-valued $F$;

• •

$D_{u_{0}}u(t)$: the Fréchet derivative of the SPDE solution with respect to the initial condition, which for the linear SPDE (1) equals the semigroup $S(t)$ (see Section 3, Eq. (9)).

2.4 The Cameron–Martin Space

A central concept in infinite-dimensional analysis is the Cameron–Martin space, which provides the natural domain for logarithmic derivatives. For a Gaussian measure $\mu$ on $H$ with covariance operator $\gamma$, the Cameron–Martin space is defined as $\mathcal{H}_{t}:=\operatorname{Ran}(\gamma^{1/2})$, equipped with the inner product

[TABLE]

where $\gamma^{-1/2}$ denotes the inverse of $\gamma^{1/2}$ on its range. When $\gamma$ is trace-class (as is the case for our Malliavin covariance operator), the Cameron–Martin space $\mathcal{H}_{t}$ is strictly smaller than $H$ (indeed, it has $\mu$-measure zero), but it is precisely the space of directions along which the Gaussian measure admits a logarithmic derivative; see (Bogachev, 2015; Da Prato and Zabczyk, 2014).

While $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$ may appear “extremely small” compared to $H$, it is in fact the maximal space on which the score is well-defined. For typical choices of $A$ and $Q$ (e.g., $A=\Delta$ with Dirichlet boundary conditions and $Q$ diagonal in the eigenbasis of $-A$), $\mathcal{H}_{t}$ consists of functions with higher Sobolev regularity than generic elements of $H=L^{2}$. In practical implementations using finite-dimensional discretisations or projections onto the first $n$ eigenmodes, the projected Cameron–Martin space becomes all of $\mathbb{R}^{n}$, and our operator formulae correspond exactly to the continuum limit.

3 Methodology

In this section, we derive a closed-form expression for the score function associated with the solution of linear SPDEs of the form (1), driven by space–time (trace class) coloured noise. Our analysis is based on the following key well-posedness assumptions:

$H$ is a separable Hilbert space, typically $L^{2}(V)$ for some spatial domain $V\subset\mathbb{R}^{d}$ (e.g., $V=(0,1)$ in one dimension), equipped with the inner product $\langle u,v\rangle_{H}$ and norm $\left\|u\right\|_{H}=\sqrt{\left\langle u,u\right\rangle_{H}}$; 2. 2.

$A:D(A)\subset H\to H$ is a densely defined, unbounded linear operator (e.g., the Laplacian $\Delta$ with Dirichlet or Neumann boundary conditions) generating a strongly continuous semigroup $S(t)=e^{tA}$ on $H$, satisfying $\left\|S(t)u\right\|_{H}\leq Me^{\omega t}\left\|u\right\|_{H}$ for some $M\geq 1$, $\omega\in\mathbb{R}$; 3. 3.

$W_{t}$ is a cylindrical Wiener process on an auxiliary separable Hilbert space $U$; we do not assume $U=H$. The process is defined on the filtered probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ and formally expressed via an orthonormal basis $\{f_{k}\}_{k=1}^{\infty}$ of $U$ as

[TABLE]

where $\{B_{t}^{k}\}_{k=1}^{\infty}$ are independent standard Brownian motions222As is well known, the series (5) may not converge in $U$ but defines a cylindrical process in a larger space.. 4. 4.

$Q^{1/2}:U\to H$ is a Hilbert-Schmidt operator with norm

[TABLE]

where $\{f_{k}\}_{k=1}^{\infty}$ is an orthonormal basis of $U$. 5. 5.

$u_{0}\in H$ is a deterministic initial condition.

The strong form, for intuition, is

[TABLE]

with $\eta(t,x)=Q^{1/2}\xi(t,x)$ and appropriate boundary conditions (e.g., $u(t,x)=0$ on $\partial V$). The mild solution is

[TABLE]

where the stochastic integral is the Hilbert space Itô integral in $L^{2}(\Omega;H)$. For well-definedness, $Q^{1/2}:U\to H$ is Hilbert-Schmidt, and we require

[TABLE]

where

[TABLE]

for an orthonormal basis $\{f_{j}\}_{j=1}^{\infty}$ of $U$.

Our objective now is to derive the Malliavin covariance operator of the SPDE solution $u(t)$. This operator measures the stochastic sensitivity of $u(t)$ with respect to perturbations in the space–time coloured noise across the entirety of $H$, providing insights into the regularity and density properties of the law of $u(t)$ (e.g., via the absolute continuity with respect to a Gaussian measure in $H$). With the noise term $Q^{1/2}dW_{t}$, this sensitivity is with respect to the coloured noise, reflecting the spatial covariance structure induced by $Q^{1/2}$, unlike the white noise case where spatial correlations are absent. This, in turn, will allow us to compute a closed-form expression for the logarithmic derivative of the transition measure of infinite-dimensional diffusion processes satisfying (1).

3.1 Malliavin Covariance Operator

To compute the Malliavin covariance operator associated with $u(t)$, we express it in terms of the first variation process, which for state-independent diffusion (i.e., $Q^{1/2}$), coincides with the semigroup $S(t)$. This mirrors the finite-dimensional stochastic differential equation (SDE) case $dx_{t}=b(x_{t})dt+\sigma dW_{t}$, where the first variation process for state-independent diffusion coincides with the solution of a deterministic linear differential equation driven by the Jacobian of the drift, i.e., $Y_{t}={\partial x_{t}}/{\partial x_{0}}$ satisfies $dY_{t}=b^{\prime}(x_{t})Y_{t}dt$, $Y_{0}=I$. For the linear SPDE (1) we have the formal solution

[TABLE]

The Fréchet derivative of $u(t)$ with respect to the initial condition $u_{0}$ is (Venturi and Dektor, 2021)

[TABLE]

Hence, the first variation process for the SPDE (1) satisfies

[TABLE]

Note that for each $v\in H$

[TABLE]

which is finite for all $t\geq 0$.

Hereafter we leverage this linearity to define the Malliavin covariance operator of $u(t)$ directly on $H$. The following theorem provides the precise formulation of this result.

Theorem 3.1**.**

Consider the linear SPDE (1) with the assumptions at the beginning of Section 3. Define the first variation process $Y_{t}=S(t)$, i.e., the semigroup generated by $A$. The Malliavin covariance operator $\gamma_{t}:H\to H$ of the solution $u(t)$ at time $t>0$ is given by

[TABLE]

where $S(s)^{*}$ is the adjoint of $S(s)$ in $\mathcal{L}(H)$. The operator $\gamma_{t}$ is positive, self-adjoint, and trace-class on $H$.

If $A$ generates a strongly continuous group $\{S(t)\}_{t\in\mathbb{R}}$, then $S(t)$ is invertible for all $t\in\mathbb{R}$ and the Malliavin covariance operator can be alternatively expressed as

[TABLE]

where $Y_{r}^{-1}=S(-r)$ and the integral $C_{t}$ is a trace-class operator on $H$.

Proof 3.2**.**

We begin by establishing existence of the mild solution to our SPDE (1). Since the equation is driven by space–time coloured noise, the noise operator $Q^{1/2}:U\to H$ is a Hilbert-Schmidt operator. The mild solution is given by

[TABLE]

where the stochastic integral is taken in $H$ with respect to the cylindrical Wiener process $W_{t}$. The well-posedness of this integral requires the integrability condition

[TABLE]

which is satisfied under our assumptions. This condition ensures that $u(t)\in L^{2}(\Omega;H)$ and that the stochastic integral is well-defined.

To derive the Malliavin covariance operator, we first compute the Malliavin derivative $\mathcal{D}_{r}$ (at time $r$) of $u(t)$. This derivative measures the sensitivity of the solution $u(t)$ with respect to perturbations in the noise $W_{r}$. Since $u_{0}$ is deterministic, we have

[TABLE]

For the stochastic integral term in (12), the Malliavin derivative at time $r\in[0,t]$ is given by the integrand evaluated at the perturbation time (see A.1)

[TABLE]

where $\mathbf{1}_{[0,t]}(r)$ is the indicator function ensuring causality. This result follows from the standard theory of Malliavin calculus for stochastic integrals: if $\Phi(s)$ is a deterministic integrand, then $\mathcal{D}_{r}\left(\int_{0}^{t}\Phi(s)dW_{s}\right)=\Phi(r)\mathbf{1}_{[0,t]}(r)$. In our case we have $\Phi(s)=S(t-s)Q^{1/2}$, yielding the expression (14). Note that for each $r\in[0,t]$, $\mathcal{D}_{r}u(t)\in\mathcal{L}_{\mathrm{HS}}(U,H)$. Next, define the Malliavin covariance operator

[TABLE]

where $(\mathcal{D}_{r}u(t))^{*}:H\to U$ denotes the adjoint and $\mathcal{D}_{r}u(t)(\mathcal{D}_{r}u(t))^{*}\in\mathcal{L}(H)$. Substituting (14) into (15) we obtain

[TABLE]

Since $Q^{1/2}:U\to H$ is a Hilbert-Schmidt operator with adjoint $(Q^{1/2})^{*}:H\to U$, we have

[TABLE]

This gives us

[TABLE]

which in turn allows us to write (15) as

[TABLE]

To obtain the final form of the Malliavin covariance operator, we perform the change of variables $s=t-r$. This gives

[TABLE]

The covariance operator $\gamma_{t}$ is positive, self-adjoint, and trace-class on $H$. The trace-class property follows from the fact that $Q^{1/2}$ is Hilbert-Schmidt, $S(s)$ and $S(s)^{*}$ are bounded operators, and the integral is over a finite interval $[0,t]$. The integrability condition (6) ensures these properties are preserved.

For the alternative representation of $\gamma_{t}$ when $A$ generates a $C_{0}$-group, we observe that if $S(t)$ forms a strongly continuous group, then $S(t)$ is invertible for all $t\in\mathbb{R}$ with $S(t)^{-1}=S(-t)$. Setting $Y_{t}=S(t)$ and $Y_{r}^{-1}=S(-r)$, we can factor the Malliavin covariance operator as

[TABLE]

where

[TABLE]

The integral $C_{t}$ in (17) represents the accumulated noise covariance transformed by the inverse semigroup, while the outer factors $Y_{t}$ and $Y_{t}^{*}$ propagate this covariance to time $t$.

To illustrate the connection between the infinite-dimensional setting discussed in Theorem 3.1 and its finite-dimensional counterpart, let $\{h_{i}\}_{i=1}^{m}\subset H$ be a set of linearly independent vectors, and define

[TABLE]

The Malliavin covariance matrix of $F$ can be computed as

[TABLE]

which provides a finite-dimensional representation of the infinite-dimensional operator $\gamma_{t}$. Note that the Malliavin covariance operator $\gamma_{t}$ encapsulates both the spatial correlation induced by $Q^{1/2}$ and the temporal evolution governed by the semigroup $S(t)$, providing a characterisation of the stochastic variability of the solution in the infinite-dimensional Hilbert space $H$.

3.2 The Cameron–Martin Space and Logarithmic Derivatives

Before proceeding to the Bismut formula, we must address a fundamental issue: in infinite dimensions, there is no Lebesgue measure on $H$, and hence the notion of a “density” $p_{u(t)}(u)$ with respect to which one computes a gradient requires careful interpretation. The appropriate framework is that of logarithmic derivatives (also called Fomin derivatives) of measures.

Let $\mu_{t}$ denote the law of $u(t)$ on $H$. Since $u(t)$ is Gaussian with mean $m_{t}:=S(t)u_{0}$ and covariance operator $\gamma_{t}$, the measure $\mu_{t}$ is a Gaussian measure $\mathcal{N}(m_{t},\gamma_{t})$ on $H$. The Cameron–Martin space associated with $\mu_{t}$ is

[TABLE]

equipped with the Cameron–Martin inner product (4) (with $\gamma=\gamma_{t}$). Note that $\gamma_{t}^{-1/2}:\mathcal{H}_{t}\to H$ is unbounded when viewed as an operator on $(H,\|\cdot\|_{H})$, but when $\mathcal{H}_{t}$ is equipped with the Cameron–Martin norm $\|h\|_{\mathcal{H}_{t}}:=\|\gamma_{t}^{-1/2}h\|_{H}$, the map $\gamma_{t}^{-1/2}:(\mathcal{H}_{t},\|\cdot\|_{\mathcal{H}_{t}})\to(H,\|\cdot\|_{H})$ is an isometry.

Definition 3.3** (Logarithmic derivative).**

Let $\mu_{t}$ be the law of $u(t)$ on $H$. For $h\in\mathcal{H}_{t}$, the logarithmic derivative of $\mu_{t}$ along $h$ is the function $\beta_{h}:H\to\mathbb{R}$ defined $\mu_{t}$-almost everywhere by the integration-by-parts formula

[TABLE]

for all $\phi\in C_{b}^{1}(H)$ (bounded continuously Fréchet differentiable functions).

The logarithmic derivative $\beta_{h}(u)$ serves as the infinite-dimensional analogue of the directional derivative $\langle\nabla\log p(u),h\rangle$ in finite dimensions. The key insight is that $\beta_{h}$ is well-defined only for directions $h\in\mathcal{H}_{t}$; for $h\notin\mathcal{H}_{t}$, the measure $\mu_{t}$ is not quasi-invariant under translation by $h$, and no logarithmic derivative exists. Importantly, this restriction is not a limitation but rather reflects the intrinsic geometry of Gaussian measures in infinite dimensions. Translations along directions outside $\mathcal{H}_{t}$ move the measure to a mutually singular measure, so there is no meaningful notion of a “score” in those directions. This is the precise infinite-dimensional analogue of the fact that in finite dimensions, the score $\nabla\log p(x)$ can only be computed at points where $p(x)>0$. For practitioners, this means when discretising to $n$ modes, the projected Cameron–Martin space $\Pi_{n}\mathcal{H}_{t}$ fills out $\mathbb{R}^{n}$ for any finite $n$, and the restriction becomes invisible. Our infinite-dimensional formulae give the continuum limit that these discretised scores converge to.

3.3 Bismut Formula

In this section, we extend the Bismut formula originally developed for finite-dimensional stochastic differential equations (Bismut, 1984; Elworthy and Li, 1994; Elworthy, 1982; Nualart and Nualart, 2018) to the infinite-dimensional setting of SPDEs of the form (1). This generalisation is then used to express the logarithmic derivative of the transition measure $\mu_{t}$, i.e., the infinite-dimensional analogue of the score function.

For directions $h\in\operatorname{Ran}(\gamma_{t})$ where we can construct an explicit covering vector field, the Bismut formula provides a stochastic representation for the logarithmic derivative. Define $G(u):=\mathbb{E}[\delta(v_{h})\mid\sigma(u(t))]$, which is a $\sigma(u(t))$-measurable random variable. The Bismut formula states that

[TABLE]

where $\delta(v_{h})$ is the Skorokhod integral of the covering vector field $v_{h}$, defined as the adjoint of the Malliavin derivative operator $\mathcal{D}_{r}$ (see A.2). See (Nualart and Nualart, 2018) for an exhaustive treatment. Unlike the Itô integral, the Skorokhod integral extends to non-adapted processes, capturing the effects of stochastic perturbations in the infinite-dimensional noise structure. The direction $h\in\operatorname{Ran}(\gamma_{t})$ corresponds to an admissible perturbation along which the explicit covering construction applies. By density of $\operatorname{Ran}(\gamma_{t})$ in $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$ (in the Cameron–Martin norm), the logarithmic derivative extends uniquely to all $h\in\mathcal{H}_{t}$ as an element of $L^{2}(\mu_{t})$; see Theorem 3.11.

Our next step is to define the covering vector field $v_{h}(r)$ appearing in the Skorokhod integral $\delta(v_{h})$ within the SPDE framework. This definition accounts for the spatially correlated noise introduced by $Q^{1/2}$, ensuring consistency with the SPDE and its mild solution in $H$.

Definition 3.4** (Covering Vector Field).**

For each $h\in\operatorname{Ran}(\gamma_{t})$, define $\tilde{h}:=\gamma_{t}^{-1}h\in H$ (noting that $\gamma_{t}^{-1}:\operatorname{Ran}(\gamma_{t})\to H$ is densely defined but unbounded). The covering vector field $v_{h}:[0,t]\to U$ is defined by

[TABLE]

Since $\operatorname{Ran}(\gamma_{t})$ is dense in $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$ with respect to the Cameron–Martin norm, the covering construction and the resulting score formula extend by continuity to all $h\in\mathcal{H}_{t}$. For practical computations (e.g., finite-rank truncations where $\gamma_{t}^{(n)}$ has finite rank), one has $\operatorname{Ran}(\gamma_{t}^{(n)})=\operatorname{Ran}((\gamma_{t}^{(n)})^{1/2})$, so this distinction disappears in approximations.

This definition ensures that $v_{h}$ is a $U$-valued process reflecting the noise’s covariance structure, with $Q^{1/2}$ modulating spatial dependence and $\gamma_{t}^{-1}h$ adjusting $h$ per the solution’s stochastic dependence. In the following theorem we prove that (22) satisfies the covering property

[TABLE]

where $H_{W}=L^{2}([0,t];U)$ is the space of square-integrable $U$-valued functions with inner product (2). Here, $\langle\mathcal{D}u(t),v_{h}\rangle_{H_{W}}$ denotes the $H$-valued Bochner integral

[TABLE]

where $\mathcal{D}_{r}u(t)\in\mathcal{L}_{\mathrm{HS}}(U,H)$ acts on $v_{h}(r)\in U$ to produce an element of $H$.

Theorem 3.5** (Covering property).**

Assume 1–5 and (6). Let $u(t)$ be the mild solution of (1), let $\gamma_{t}$ be given by (11), and let $H_{W}=L^{2}([0,t];U)$. For each $h\in\operatorname{Ran}(\gamma_{t})$, define $v_{h}$ by (22) with $\tilde{h}=\gamma_{t}^{-1}h$. Then the covering property holds:

[TABLE]

Proof 3.6**.**

Consider the space $H_{W}=L^{2}([0,t];U)$, where $U$ is the Hilbert space representing the input noise. In the context of the Malliavin derivative $\mathcal{D}_{r}u(t)$, which in this linear additive setting is a deterministic function of $r$ with values in $\mathcal{L}_{\mathrm{HS}}(U,H)$, and the vector field $v_{h}\in H_{W}$, the pairing $\left\langle\mathcal{D}u(t),v_{h}\right\rangle_{H_{W}}$ is defined as in (24).

By (14), $\mathcal{D}_{r}u(t)=S(t-r)Q^{1/2}\mathbf{1}_{[0,t]}(r)$. Let $h\in\operatorname{Ran}(\gamma_{t})$ and write $\tilde{h}=\gamma_{t}^{-1}h$. Substituting the covering vector field (22) into the pairing (24) we obtain

[TABLE]

Here, $S(t-r)Q^{1/2}$ maps $U$ to $H$, and $\left(Q^{1/2}\right)^{*}S(t-r)^{*}\tilde{h}$ is a vector in $U$ ensuring the composition is well‑defined. To evaluate (25), we perform the change of variables $s=t-r$. This yields

[TABLE]

From Theorem 3.1, we know that the Malliavin covariance operator can be expressed as in (11). Thus, the pairing simplifies to

[TABLE]

Theorem 3.7** (Solution structure and minimal norm).**

Under the assumptions of Theorem 3.5, define the subspace $\mathcal{V}:=\{(\mathcal{D}u(t))^{*}z:z\in H\}\subset H_{W}$ and fix $h\in\operatorname{Ran}(\gamma_{t})$. Then:

There is a unique $v\in\mathcal{V}$ with $\langle\mathcal{D}u(t),v\rangle_{H_{W}}=h$, namely $v_{h}$ as defined in (22). 2. 2.

Among all solutions $v\in H_{W}$ of the covering equation $\langle\mathcal{D}u(t),v\rangle_{H_{W}}=h$, the covering vector field $v_{h}$ is the unique solution of minimal $H_{W}$-norm. 3. 3.

Every solution in $H_{W}$ is of the form $v_{h}+w$ with $w\in\mathcal{V}^{\perp}$, and

[TABLE]

Proof 3.8**.**

(1) Let $v\in\mathcal{V}$ solve the covering equation. Then $v=(\mathcal{D}u(t))^{*}z$ for some $z\in H$, and

[TABLE]

For $h\in\operatorname{Ran}(\gamma_{t})$, there exists unique $\tilde{h}=\gamma_{t}^{-1}h\in H$ with $\gamma_{t}\tilde{h}=h$. Then $\gamma_{t}z=h$ implies $z-\tilde{h}\in\ker(\gamma_{t})$. Since $\gamma_{t}=\mathcal{D}u(t)(\mathcal{D}u(t))^{*}$ is positive self-adjoint, $\ker\gamma_{t}=\ker(\mathcal{D}u(t))^{*}$. Thus

[TABLE]

as $(\mathcal{D}u(t))^{*}$ annihilates $\ker\gamma_{t}$. Therefore $v=v_{h}$, proving uniqueness in $\mathcal{V}$.

(2) Let $T:H_{W}\to H$ be $Tv:=\langle\mathcal{D}u(t),v\rangle_{H_{W}}$ (Bochner integral). Then $\ker T=\mathcal{V}^{\perp}$ because for any $z\in H$,

[TABLE]

Thus every solution $v$ can be written uniquely as $v=v_{\parallel}+v_{\perp}$ with $v_{\parallel}\in\mathcal{V}$, $v_{\perp}\in\mathcal{V}^{\perp}$ and $Tv_{\perp}=0$, $Tv_{\parallel}=h$. Among all such decompositions,

[TABLE]

is minimised iff $v_{\perp}=0$. By part (1), $v_{\parallel}=v_{h}$ is unique, so the minimal-norm solution is uniquely $v_{h}$.

(3) If $v$ is any solution, then $T(v-v_{h})=0$, so $v-v_{h}\in\ker T=\mathcal{V}^{\perp}$. Thus $v=v_{h}+w$ for some $w\in\mathcal{V}^{\perp}$. The norm identity follows from orthogonality $\mathcal{V}\perp\mathcal{V}^{\perp}$.

3.4 Skorokhod Integral

In this section we show that, for the linear SPDE (1) with state-independent colouring $Q^{1/2}$ (Hilbert–Schmidt), the Skorokhod integral simplifies to the more tractable Itô integral. To this end, let us first recall the Bismut formula (21), which for $h\in\operatorname{Ran}(\gamma_{t})$ states that

[TABLE]

expressing the logarithmic derivative in terms of the conditional expectation of the Skorokhod integral $\delta(v_{h})$. The Skorokhod integral $\delta:\mathrm{Dom}\,(\delta)\subset L^{2}(\Omega;H_{W})\to L^{2}(\Omega)$ is defined as the adjoint of the Malliavin derivative $\mathcal{D}$ via the duality relation

[TABLE]

for all $F\in\mathbb{D}^{1,2}$ and $v\in\mathrm{Dom}\,(\delta)$; see A.2 for details. For deterministic or adapted integrands $v_{h}\in H_{W}$ satisfying $\mathbb{E}[\int_{0}^{t}\|v_{h}(r)\|_{U}^{2}\,dr]<\infty$, the Skorokhod integral coincides with the Itô integral:

[TABLE]

Our goal now is to show that, due to the linearity of the drift term $Au(t)$ in the SPDE (1) and the state-independence of the noise $Q^{1/2}dW_{t}$, the covering vector field $v_{h}(r)$ is adapted to the filtration $\{\mathcal{F}_{r}\}_{r\geq 0}$. In fact, $v_{h}(r)$ is deterministic, which implies adaptedness trivially. This simplifies the Skorokhod integral to an Itô integral with respect to the cylindrical Wiener process $W_{t}$

[TABLE]

A substitution of (22) into (30) yields

[TABLE]

where $\tilde{h}=\gamma_{t}^{-1}h\in H$ for $h\in\operatorname{Ran}(\gamma_{t})$.

In the next theorem we show that the covering vector field $v_{h}(r)$ associated with the linear SPDE (1) is indeed deterministic and adapted to the filtration $\{\mathcal{F}_{r}\}_{r\geq 0}$. Consequently, the Skorokhod integral (29) coincides with the Itô integral (31).

Theorem 3.9**.**

The covering vector field (22) for the linear SPDE (1) is a deterministic, $U$-valued process that is adapted to the filtration $\{\mathcal{F}_{r}\}_{r\geq 0}$, thereby satisfying the necessary conditions for Itô integration in the infinite-dimensional setting. Consequently, the Skorokhod integral (29) coincides with the Itô integral

[TABLE]

Moreover, $v_{h}\in\mathrm{Dom}\,(\delta)$ with $\mathbb{E}[\|\delta(v_{h})\|^{2}]=\|v_{h}\|_{H_{W}}^{2}<\infty$.

Proof 3.10**.**

Consider the covering vector field $v_{h}(r)$ defined in (22), where $\tilde{h}=\gamma_{t}^{-1}h\in H$ for $h\in\operatorname{Ran}(\gamma_{t})$. To evaluate the properties of $v_{h}(r)$ in the space $U$, we first notice that $\tilde{h}\in H$ is deterministic, and $S(t-r)^{*}:H\to H$ can be bounded as

[TABLE]

Since $\left(Q^{1/2}\right)^{*}:H\to U$ is bounded (adjoint of a Hilbert-Schmidt operator), with

[TABLE]

we have

[TABLE]

This norm is a deterministic function of $r$, finite for all $r\in[0,t]$, and continuous due to the strong continuity of $S(s)^{*}$. Since $\left(Q^{1/2}\right)^{*}$, $S(t-r)^{*}$, and $\tilde{h}\in H$ are all deterministic, we have that $v_{h}(r)$ is deterministic. A deterministic process is trivially adapted to any filtration, including $\{\mathcal{F}_{r}\}_{r\geq 0}$, because $v_{h}(r)$ is measurable with respect to the trivial $\sigma$-algebra $\{\emptyset,\Omega\}$, which is contained in $\mathcal{F}_{r}$ for all $r\geq 0$. Thus, $v_{h}(r)$ is both deterministic and $\{\mathcal{F}_{r}\}$-adapted, satisfying the prerequisites for Itô integration.

Finally, we show that the Skorokhod integral (29) reduces to the Itô integral (30). A key result in (Nualart, 2006) states that if $v_{h}\in\text{Dom}(\delta)$ and satisfies both $v_{h}(r)$ being $\{\mathcal{F}_{r}\}$-adapted and $\mathbb{E}\left[\int_{0}^{t}\|v_{h}(r)\|_{U}^{2}\,dr\right]<\infty$, then $\delta(v_{h})$ coincides with the Itô integral. We have established already that $v_{h}(r)$ is deterministic and thus adapted to $\{\mathcal{F}_{r}\}$. For square-integrability, compute the expectation

[TABLE]

since $v_{h}(r)$ is deterministic, reducing the expectation to the deterministic norm squared. Estimating,

[TABLE]

Integrating over $[0,t]$

[TABLE]

and using

[TABLE]

we obtain

[TABLE]

which is finite, as $\left\|Q^{1/2}\right\|_{\mathcal{L}_{\mathrm{HS}}(U,H)}$, $M$, $\left\|\tilde{h}\right\|_{H}$, and $t$ are all finite. Thus, $v_{h}\in L^{2}([0,t]\times\Omega;U)$, satisfying the integrability condition. Since $v_{h}(r)$ is adapted and square-integrable, the Skorokhod integral reduces to the Itô integral (30) (or (31)). Next, recall that the cylindrical Wiener process $W_{t}$ admits the representation (5), where $\{B^{k}_{t}\}_{k=1}^{\infty}$ are independent Brownian motions. Substituting (5) into (30) yields

[TABLE]

This integral is well-defined if

[TABLE]

Recalling the definition of $v_{h}(r)$ (22)

[TABLE]

which we have shown to be integrable. Thus,

[TABLE]

where the integral is well-defined as an Itô integral, with each component deterministic, confirming the reduction. This completes the proof.

3.5 Logarithmic Derivative (Score Function)

In this section, we derive a closed-form expression for the logarithmic derivative associated with the solution of linear SPDEs of the form (1), driven by space–time coloured noise. In particular, we leverage the Bismut formula discussed in Section 3.3 and Theorem 3.9 to express the logarithmic derivative in terms of the Malliavin covariance operator and the first variation process $Y_{t}=S(t)$.

Theorem 3.11** (Logarithmic derivative / Score function).**

Consider the linear SPDE (1) with the assumptions at the beginning of Section 3, and let the Malliavin covariance operator $\gamma_{t}$ be injective333For instance, $\gamma_{t}$ is injective if the closed linear span of $\bigcup_{0\leq s\leq t}\operatorname{Ran}(S(s)Q^{1/2})$ equals $H$ (equivalently, the controllability Gramian has trivial kernel).. Let $\mu_{t}=\mathcal{N}(S(t)u_{0},\gamma_{t})$ denote the Gaussian law of $u(t)$ on $H$. Then for each $h\in\operatorname{Ran}(\gamma_{t})$, the logarithmic derivative of $\mu_{t}$ along $h$ is given by

[TABLE]

where $\gamma_{t}^{-1}:\operatorname{Ran}(\gamma_{t})\to H$ is the (unbounded, densely defined) inverse of $\gamma_{t}$ on its range.

The map $h\mapsto\beta_{h}$ is continuous from $(\mathcal{H}_{t},\|\cdot\|_{\mathcal{H}_{t}})$ into $L^{2}(\mu_{t})$. By density of $\operatorname{Ran}(\gamma_{t})$ in $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$, the logarithmic derivative extends uniquely to all $h\in\mathcal{H}_{t}$ as an element of $L^{2}(\mu_{t})$, given by the eigenseries formula (40).

Proof 3.12**.**

We derive the logarithmic derivative using the Bismut–Elworthy–Li formula within the infinite-dimensional Malliavin calculus framework, adapting the methodology to the coloured noise case driven by $Q^{1/2}dW_{t}$.

In the proof of Theorem 3.1 we established that the Malliavin derivative of the solution $u(t)$ with respect to the noise at time $r\in[0,t]$ is

[TABLE]

where $\mathbf{1}_{[0,t]}(r)$ is the indicator function ensuring causality. Here $\mathcal{D}_{r}u(t)\in\mathcal{L}_{\mathrm{HS}}(U,H)$. The Malliavin covariance operator $\gamma_{t}$ (see Eq. (15)) is positive, self-adjoint, and trace-class under the stated assumptions on $Q^{1/2}$ and $S(t)$, as it arises from the integral of positive semi-definite operators in the context of SPDEs with additive noise. The injectivity assumption on $\gamma_{t}$ (equivalently, $\ker\gamma_{t}=\{0\}$) ensures that all finite-dimensional projections $\Pi u(t)$ have non-degenerate covariance matrices, and hence smooth densities with respect to Lebesgue measure on the range of $\Pi$. Since $\gamma_{t}$ is trace-class (hence compact), the inverse $\gamma_{t}^{-1}$ is unbounded as an operator on $H$, even when $\gamma_{t}$ is injective. Injectivity together with the trace-class property implies that $\gamma_{t}$ has dense range (i.e., $\overline{\operatorname{Ran}(\gamma_{t})}=H$), but the eigenvalues of $\gamma_{t}$ accumulate at zero, precluding any uniform positive lower bound. This is why $\gamma_{t}^{-1}:\operatorname{Ran}(\gamma_{t})\to H$ is densely defined but unbounded. For linear SPDEs with additive noise, Malliavin differentiability of $u(t)$ follows directly from the explicit representation (7), and no Hörmander-type condition is needed. The injectivity of $\gamma_{t}$ is equivalent to the approximate controllability condition: the closed linear span of $\bigcup_{0\leq s\leq t}\operatorname{Ran}(S(s)Q^{1/2})$ equals $H$. This ensures that noise propagates to all directions in $H$. For a direction $h\in\operatorname{Ran}(\gamma_{t})$, let $\tilde{h}=\gamma_{t}^{-1}h\in H$. The covering vector field (Definition 3.4) is

[TABLE]

which is a $U$-valued, deterministic process. By Theorem 3.5, this satisfies the covering property $\langle\mathcal{D}u(t),v_{h}\rangle_{H_{W}}=h$.

Since $v_{h}(r)$ is deterministic and adapted, as shown in Theorem 3.9, the Skorokhod integral coincides with the Itô integral. The mild solution representation (7) gives

[TABLE]

where $z$ is a centred Gaussian random variable in $H$ with covariance $\gamma_{t}$.

For any $\xi\in H$, we have

[TABLE]

Taking $\xi=\tilde{h}=\gamma_{t}^{-1}h$:

[TABLE]

By the Bismut formula (21), defining $G:=\mathbb{E}[\delta(v_{h})\mid\sigma(u(t))]$, we have $\beta_{h}(u(t))=-G$ as $\sigma(u(t))$-measurable random variables. Since $\delta(v_{h})=\langle z,\tilde{h}\rangle_{H}$ (by (39)) and $z=u(t)-S(t)u_{0}$, the conditional expectation simplifies:

[TABLE]

since $\langle u(t)-S(t)u_{0},\tilde{h}\rangle_{H}$ is already $\sigma(u(t))$-measurable. Thus

[TABLE]

which is formula (36).

For the extension to all $h\in\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$: since $\operatorname{Ran}(\gamma_{t})$ is dense in $\mathcal{H}_{t}$ with respect to the Cameron–Martin norm, we extend by continuity in $L^{2}(\mu_{t})$. Let $\{e_{k}\}_{k=1}^{\infty}$ be an orthonormal basis of $H$ consisting of eigenvectors of $\gamma_{t}$, with $\gamma_{t}e_{k}=\lambda_{k}e_{k}$ and $\lambda_{k}>0$ (by injectivity). For $h\in\mathcal{H}_{t}$, write $h=\sum_{k=1}^{\infty}h_{k}e_{k}$ where $\sum_{k=1}^{\infty}h_{k}^{2}/\lambda_{k}<\infty$ (the Cameron–Martin condition). The logarithmic derivative is

[TABLE]

where the series converges in $L^{2}(\mu_{t})$.

In finite dimensions where $H=\mathbb{R}^{n}$, the covariance matrix $\gamma_{t}$ is positive definite (under our injectivity assumption), so $\gamma_{t}^{-1}\in\mathcal{L}(\mathbb{R}^{n})$ is a bounded operator and $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t})=\mathbb{R}^{n}$. In this case, formula (36) can be written equivalently as

[TABLE]

which is the standard score function for a Gaussian distribution with mean $S(t)u_{0}$ and covariance $\gamma_{t}$. In infinite dimensions, such a formula is not well-posed because $\gamma_{t}^{-1}$ is unbounded and $(u-S(t)u_{0})\notin\operatorname{Ran}(\gamma_{t})$ generically. The correct formulation is (36), where the unbounded inverse $\gamma_{t}^{-1}$ acts on the direction $h\in\operatorname{Ran}(\gamma_{t})$, not on the random state.

One might ask: if $u(t)$ is Gaussian, why employ Malliavin calculus? The answer is twofold. First, the Bismut formula expresses the score as a conditional expectation $\beta_{h}=-\mathbb{E}[\delta(v_{h})\mid\sigma(u(t))]$, which provides an intrinsic representation that does not rely on the explicit Gaussian form; this representation extends naturally to settings where the law of $u(t)$ is non-Gaussian. Second, working with Cameron–Martin spaces and the eigenseries representation (40) ensures an intrinsic infinite-dimensional formulation that is mathematically rigorous and independent of any discretisation scheme. The Malliavin/Bismut machinery presented here could serve as a starting point for extensions to nonlinear SPDEs with multiplicative noise, where the law of $u(t)$ is no longer Gaussian and the score becomes a genuine conditional expectation; such extensions, however, remain the subject of future work.

4 Numerical Results

In this section, we validate the closed-form Malliavin score formula

[TABLE]

we obtained in Theorem 3.11 through numerical simulations of SPDEs in one- and two-dimensional spatial domains. Specifically, in Section 4.1, we consider one-dimensional SPDEs on a bounded domain with Dirichlet boundary conditions, discretised via spectral Galerkin truncation in the Laplacian eigenbasis. In Section 4.2, we extend the validation to two-dimensional SPDEs on a periodic domain, discretised using a Fourier spectral method. These experiments confirm that the score formula applies to both second- and fourth-order operators, independently of the spatial dimension, boundary conditions, and choice of spectral basis.

4.1 One Dimension

All numerical experiments in this section are conducted on one-dimensional linear SPDEs of the form (1) defined on the spatial domain $(0,1)$ with homogeneous Dirichlet boundary conditions. We use second-order finite-differences to approximate the directional derivative of the log-density, against which we validate the closed-form formula (41).

4.1.1 Galerkin discretisation

We consider the state space $H=L^{2}(0,1)$ with homogeneous Dirichlet boundary conditions, for which the Laplacian $\Delta$ admits the orthonormal eigenbasis $\{\varphi_{k}\}_{k=1}^{\infty}$ given by

[TABLE]

For operators $A$ that are functions of the Laplacian, i.e., $A=f(-\Delta)$ for some function $f:\mathbb{R}_{+}\to\mathbb{R}$, the eigenbasis (42) diagonalises $A$ with eigenvalues $a_{k}=f(\lambda_{k})$.

The spectral Galerkin approximation truncates the expansion to $N$ modes. In eigenspace coordinates, the SPDE (1) decouples into $N$ independent scalar Ornstein–Uhlenbeck processes:

[TABLE]

where $u_{k}(t)=\langle u(t),\varphi_{k}\rangle_{H}$ are the Fourier coefficients, $\{B_{k}(t)\}_{k=1}^{N}$ are independent standard Brownian motions, and $q_{k}$ are the eigenvalues of the noise covariance operator $Q$. We take $q_{k}=k^{-2}$, ensuring trace-class regularity $\sum_{k=1}^{\infty}q_{k}<\infty$. The rapid algebraic decay of $q_{k}$ ensures that the Galerkin truncation to $N=64$ modes fully resolves the solution. The total solution variance $\mathrm{Tr}(\gamma_{t})$ is identical at $N=64$ and $N=128$ to six significant figures, with the residual noise variance in modes $k>64$ accounting for less than $0.5\%$ of the total. Increasing the truncation level produces no visible change in the solution fields or the score. The mild solution admits the representation

[TABLE]

which is Gaussian with mean $m_{k}(t)=e^{a_{k}t}u_{0,k}$, where $u_{0,k}=\langle u_{0},\varphi_{k}\rangle_{H}$, and variance

[TABLE]

with the convention $\gamma_{k}(t)=q_{k}t$ when $a_{k}=0$. The Malliavin covariance operator $\gamma_{t}$ is diagonal in the eigenbasis with eigenvalues (45).

4.1.2 1D SPDE classes

We validate the score formula across four classes of linear SPDEs. Table 1 summarises the operator structure and eigenvalues for each class.

The first three operators correspond to classical diffusion processes: Brownian diffusion (heat equation); the Ornstein–Uhlenbeck process (heat equation with linear damping leading to faster relaxation towards equilibrium); and scaled Brownian diffusion (heat equation with reduced diffusivity). The fractional Laplacian with $\alpha=0.75$ (order $2\alpha=1.5$) models anomalous subdiffusion.

4.1.3 Validation methodology

For each SPDE class, we validate the Malliavin score formula (41) against a central finite-difference approximation of the directional derivative. In the $N$-dimensional Galerkin truncation, the projected solution $u^{(N)}(t)\in\mathbb{R}^{N}$ admits a Lebesgue density $p_{t}^{(N)}$. For a direction $h\in H$ with $\|h\|_{H}=1$, the finite-difference approximation is

[TABLE]

with step size $\varepsilon=10^{-5}$. Since $u^{(N)}(t)\sim\mathcal{N}(m_{t},\gamma_{t})$ is Gaussian, the log-density (up to an additive constant) is

[TABLE]

and the finite-difference approximation (46) can be computed exactly in eigenspace. The Malliavin score formula in eigenspace coordinates reads

[TABLE]

where $h_{k}=\langle h,\varphi_{k}\rangle_{H}$ are the Fourier coefficients of the direction $h$. We take $h=\varphi_{1}$ (the first eigenmode) throughout, corresponding to the lowest-frequency perturbation.

For each SPDE class listed in Table 1, we simulate $P=4$ independent sample paths of the solution process $\{u(t)\}_{t\in[0,T]}$ with $T=1$ and initial condition $u_{0}=\varphi_{1}+\tfrac{1}{2}\varphi_{2}+\tfrac{1}{4}\varphi_{3}$. The truncation level is $N=64$ modes, and we evaluate the score at $50$ uniformly spaced time points $t_{j}\in[0.02,1]$. At each time point and along each sample path, we compute:

The Malliavin score $\beta_{h}(u(t_{j}))$ via formula (48); 2. 2.

The finite-difference score $\beta_{h}^{\mathrm{FD}}(u(t_{j}))$ via formula (46).

The absolute error $|\beta_{h}(u)-\beta_{h}^{\mathrm{FD}}(u)|$ quantifies the discrepancy between the two methods.

4.1.4 Results

Figure 1 presents the numerical validation results for the four SPDE classes discussed in Section 4.1.2.

It is seen that for the heat equation, Ornstein–Uhlenbeck, scaled diffusion, and fractional Laplacian SPDEs, the Malliavin and finite-difference scores exhibit excellent agreement, with errors in the range $10^{-11}$ to $10^{-7}$. This is consistent with the expected truncation error of the central finite-difference scheme: for a function $f$ with bounded fourth derivative, $|f^{\prime}(x)-(f(x+\varepsilon)-f(x-\varepsilon))/(2\varepsilon)|=\mathcal{O}(\varepsilon^{2})$, giving $\mathcal{O}(10^{-10})$ for $\varepsilon=10^{-5}$. The slightly larger observed errors arise from amplification by the curvature of the log-density, which scales with $\gamma_{k}^{-1}$. These results demonstrate that Theorem 3.11 applies to any linear SPDE of the form (1) where $A$ is diagonalisable in a common eigenbasis with the noise covariance $Q$. The specific physics encoded in the operator $A$, whether diffusion, damping, or higher-order smoothing, enters only through the eigenvalues $a_{k}$, which determine the Malliavin covariance (45). The score formula itself is structurally identical across all cases.

4.2 Two Dimensions

The one-dimensional experiments of Section 4.1 employed spectral Galerkin discretisation in the sine eigenbasis, which diagonalises the operator $A$ by construction. To confirm that the Malliavin score formula applies independently of the spatial dimension, boundary conditions, and choice of spectral basis, we now validate it on the two-dimensional periodic domain $[0,2\pi]^{2}$ using a Fourier spectral method with both second- and fourth-order operators.

4.2.1 Fourier spectral discretisation

We consider the uniform grid on $[0,2\pi]$ given by $x_{j}=2\pi j/(N+1)$ for $j=0,\ldots,N$ with $N=48$, and construct the associated two-dimensional tensor-product grid on $[0,2\pi]\times[0,2\pi]$, yielding $49\times 49=2{,}401$ degrees of freedom. Since the operators in Table 2 depend only on $|k|^{2}=k_{1}^{2}+k_{2}^{2}$, they are diagonal in the Fourier basis, with analytically known eigenvalues $a_{k_{1},k_{2}}$. All computations—including eigenvalue evaluation, covariance construction, sampling, and score evaluation—are performed mode-by-mode via the two-dimensional FFT, resulting in a Fourier spectral (Galerkin) discretisation on the periodic domain. As an independent cross-validation, we construct the physical-space differentiation matrix $D\in\mathbb{R}^{(N+1)\times(N+1)}$ for the odd-point trigonometric interpolant

[TABLE]

and assemble the two-dimensional Laplacian and biharmonic operators via Kronecker products

[TABLE]

where $D^{(2)}=D^{2}$ and $D^{(4)}=(D^{2})^{2}$. The numerically computed spectrum of these matrices agrees with the analytical eigenvalues to within $\mathcal{O}(10^{-4})$, confirming the correctness of the spectral implementation.

4.2.2 2D SPDE classes

We validate the score formula across four linear SPDEs on $[0,2\pi]^{2}$ with periodic boundary conditions, encompassing both second- and fourth-order operators.

The stochastic heat equation and Ornstein–Uhlenbeck process are second-order operators representing classical diffusion. The stochastic biharmonic equation is fourth-order and models surface diffusion and thin-film dynamics. The linearised Swift–Hohenberg equation combines second- and fourth-order terms and arises in pattern formation. Noise covariance eigenvalues are $q_{k_{1},k_{2}}=(1+k_{1}^{2}+k_{2}^{2})^{-2}$, which is trace-class in two dimensions.

4.2.3 Spectral properties of the forcing term and resolution study

Before validating the score formula, we examine the spectral properties of the coloured noise to confirm that the Fourier truncation to $N+1=49$ modes per direction adequately resolves the forcing. In Figure 2 we plot four independent realisations of the noise field $Q^{1/2}\xi(x,y)$ at a fixed time. The samples exhibit the characteristic structure of coloured noise: random spatial patterns with $\mathcal{O}(1)$ amplitude and a smooth appearance, reflecting the rapid spectral decay $q_{k_{1},k_{2}}=(1+|k|^{2})^{-2}\sim|k|^{-4}$ at high wavenumbers. The PDE solution is smoother still. For the stochastic heat equation, the semigroup $S(t)=e^{t\Delta}$ further damps mode $k$ by a factor $e^{-\nu|k|^{2}t}$, so the per-mode variance of the stochastic convolution decays as $\gamma_{k}(t)\sim|k|^{-2s-2}$ for large $|k|$, where $s$ is the exponent in the noise covariance $q_{k}=(1+|k|^{2})^{-s}$ (here $s=2$). For the fourth-order operators the decay is $|k|^{-2s-4}$. Thus, the coloured noise that ensures well-posedness of the SPDE also acts as a spectral filter, and the semigroup provides further smoothing. To verify that the Fourier truncation is sufficient, we compare the total solution variance $\mathrm{Tr}(\gamma_{t})$ at two resolutions, $49\times 49$ modes (used in our experiments) and $99\times 99$ modes. For the stochastic heat equation at $t=1$, the values are $\mathrm{Tr}(\gamma_{t})=1.597227$ and $1.597228$, respectively—identical to six significant figures. The variance carried by modes beyond $|k|>24$ (i.e., those added by the higher resolution) accounts for less than $0.0001\%$ of the total. Likewise, the noise variance beyond $|k|>24$ is less than $0.1\%$ of the total. These results confirm that with the spectral decay $s=2$, the truncation to $49$ modes per direction fully resolves both the forcing and the solution, and increasing the resolution produces no change in the computed fields or score errors.

4.2.4 Validation methodology

Since the SPDE (1) is linear with additive Gaussian noise, the solution $u(t)$ is Gaussian with known mean and covariance in Fourier space. Each Fourier mode $\hat{u}_{k_{1},k_{2}}(t)$ is an independent complex Gaussian with mean

[TABLE]

and variance

[TABLE]

where $\hat{u}_{0,k_{1},k_{2}}$ denotes the Fourier coefficient of the initial condition. We take $u_{0}(x,y)=\cos x+\tfrac{1}{2}\cos y+\tfrac{1}{4}\cos(x+y)$ and evaluate at $t=1$. The modes are sampled analytically via the FFT, avoiding time-stepping errors entirely. For each SPDE class, we compute two quantities at each grid point $(x_{i},y_{j})$:

The Malliavin score from Theorem 3.11. Note that in Fourier space we have $\hat{s}_{k_{1},k_{2}}=-(\hat{u}_{k_{1},k_{2}}-\hat{m}_{k_{1},k_{2}})/\gamma_{k_{1},k_{2}}$. This can be mapped to physical space via the inverse FFT. 2. 2.

A per-mode finite-difference approximation of the score in Fourier space. For each mode $(k_{1},k_{2})$, the log-density contribution is $-|\hat{c}_{k_{1},k_{2}}|^{2}/(2\gamma_{k_{1},k_{2}})$, where $\hat{c}_{k_{1},k_{2}}=\hat{u}_{k_{1},k_{2}}-\hat{m}_{k_{1},k_{2}}$. Writing $\hat{c}=\alpha+i\beta$ (suppressing mode indices), the central finite difference is applied to the real and imaginary parts separately

[TABLE]

Using the algebraic identity $(x+\varepsilon)^{2}-(x-\varepsilon)^{2}=4\varepsilon x$, which is exact for quadratics, this simplifies to $\hat{s}^{\mathrm{FD}}_{k_{1},k_{2}}=-\hat{c}_{k_{1},k_{2}}/\gamma_{k_{1},k_{2}}$ for any $\varepsilon>0$, yielding zero truncation error and avoiding catastrophic cancellation in high-frequency modes. The physical-space FD score is then obtained via the inverse FFT of $\hat{s}^{\mathrm{FD}}$.

The pointwise error field $|s_{\mathrm{Mall}}(x,y)-s_{\mathrm{FD}}(x,y)|$ is computed from a single inverse FFT of the Fourier-space difference, ensuring that the physical-space error reflects only the per-mode discrepancy.

4.2.5 Results

Figure 3 displays the stochastic component $(u-S(t)u_{0})(x,y)$ for each of the four SPDEs at $t=1$. The heat equation and Ornstein–Uhlenbeck fields show mid-frequency fluctuations, while the fourth-order operators (biharmonic, Swift–Hohenberg) produce smoother fields due to stronger damping of high-frequency modes.

The visual smoothness of the stochastic component is a direct consequence of the trace-class noise assumption that underpins the well-posedness of (1). For the stochastic heat equation, the per-mode variance of the stochastic convolution decays as $\gamma_{k_{1},k_{2}}(t)\sim q_{k_{1},k_{2}}/(2\nu|k|^{2})\sim|k|^{-2s-2}$ for large $|k|^{2}=k_{1}^{2}+k_{2}^{2}$, where $s$ is the exponent in the noise covariance $q_{k_{1},k_{2}}=(1+|k|^{2})^{-s}$. With $s=2$, this gives $\gamma_{k}\sim|k|^{-6}$, which by the Sobolev embedding theorem in two dimensions places the sample paths almost surely in $H^{\alpha}$ for $\alpha<s=2$, and hence in $C^{0}([0,2\pi]^{2})$. The same mechanism applies to all four operators, with the fourth-order cases exhibiting even faster spectral decay. This regularity is characteristic of SPDEs driven by spatially correlated noise with trace-class covariance, where the coloured noise that ensures well-posedness in $H$ simultaneously regularises the solution (Da Prato and Zabczyk, 2014; Ferrante and Sanz-Solé, 2006). In contrast, space–time white noise ($Q=I$, $s=0$) does not produce function-valued solutions in $d\geq 2$ (Hairer, 2009). Thus, the observed smoothness is not a numerical artefact but a genuine feature of the coloured-noise regime studied in this paper.

Figure 4 shows the corresponding pointwise score error field $|s_{\mathrm{Mall}}(x,y)-s_{\mathrm{FD}}(x,y)|$. The error is spatially unstructured and lies at machine precision throughout: $\mathcal{O}(10^{-10})$ for the second-order operators and $\mathcal{O}(10^{-9})$ for the fourth-order operators, where the larger eigenvalues (${\sim}10^{6}$) amplify floating-point rounding. In Fourier space, the maximum error across all modes is $\mathcal{O}(10^{-11})$ for second-order and $\mathcal{O}(10^{-10})$ for fourth-order operators. These results confirm that Theorem 3.11 applies to arbitrary linear SPDEs of the form (1) on periodic domains in two spatial dimensions, independently of the order of the differential operator and the discretisation method.

5 Summary

We studied score-based diffusion models in infinite-dimensional separable Hilbert spaces using Malliavin calculus. By formulating the forward diffusion process as a linear SPDE driven by space–time coloured noise with a trace-class covariance operator, we ensured mathematical well-posedness across arbitrary spatial dimensions. Our derivation of the logarithmic derivative of the transition measure, the natural infinite-dimensional analogue of the score function, uses Malliavin calculus and an infinite-dimensional generalisation of the Bismut–Elworthy–Li formula, yielding a closed-form expression along Cameron–Martin directions without relying on finite-dimensional projections or approximations. This operator-theoretic approach preserves the intrinsic structure of Hilbert spaces and accommodates general trace-class operators, incorporating spatially correlated noise without assuming semigroup invertibility.

A key insight of our analysis is that the score (logarithmic derivative) is naturally defined only along directions in the Cameron–Martin space $\mathcal{H}_{t}=\operatorname{Ran}(\gamma_{t}^{1/2})$, which is strictly smaller than $H$. While this may appear restrictive, it is in fact the maximal domain on which the score is meaningful; translations outside $\mathcal{H}_{t}$ move the Gaussian measure to a mutually singular measure. In practical discretisations, this restriction becomes invisible as the projected Cameron–Martin space fills the finite-dimensional approximation space.

We validated the score formula numerically for four classes of linear SPDEs in one spatial dimension (spectral Galerkin discretisation with Dirichlet boundary conditions) and four classes in two spatial dimensions (Fourier spectral discretisation with periodic boundary conditions), the latter including both second- and fourth-order operators. In all cases, the Malliavin score agrees with finite-difference approximations to machine precision.

Acknowledgements.

Daniele Venturi was supported by the U.S. Department of Energy (DOE) under grant DE–SC0024563.

Appendix A Malliavin Calculus

In this appendix, we provide a brief overview of Malliavin calculus for linear SPDEs of the form (1). To this end, let $W_{t}$ be a cylindrical Wiener process on a separable Hilbert space $U$. For any $\Phi\in L^{2}([0,t];\mathcal{L}_{\mathrm{HS}}(U,H))$ we define the stochastic integral $\int_{0}^{t}\Phi(s)\,dW_{s}$ as

[TABLE]

where $\{f_{j}\}$ is an orthonormal basis of $U$ and $B_{j}(s)=\left\langle W_{s},f_{j}\right\rangle_{U}$ are independent Brownian motions. The integral (50) is well-defined in $L^{2}(\Omega;H)$ and it satisfies Itô’s isometry

[TABLE]

For the SPDE (1), we have $\Phi(s)=S(t-s)Q^{1/2}$, and the mild solution

[TABLE]

Lemma A.1**.**

Under condition (6), the stochastic convolution $W_{A}(t)=\int_{0}^{t}S(t-s)Q^{1/2}\,dW_{s}$ is well-defined in $L^{2}(\Omega;H)$ with

[TABLE]

where $\gamma_{t}$ is the Malliavin covariance operator (16).

Proof A.2**.**

By the Itô isometry and the cyclic property of trace

[TABLE]

This completes the proof.

Let us now define cylindrical Wiener processes, the Cameron–Martin space, and cylindrical functionals of the Wiener process. For each $u\in U$ the cylindrical Wiener process on $U$ is the real-valued Brownian motion

[TABLE]

As is well known, the variance and covariance of $W_{t}(u)$ are, respectively

[TABLE]

Definition A.3** (Cameron–Martin space).**

The Cameron–Martin space is defined as

[TABLE]

with inner product $\left\langle h,g\right\rangle_{H_{W}}=\int_{0}^{t}\left\langle h(s),g(s)\right\rangle_{U}\,ds.$

For $h\in H_{W}$, the Wiener integral

[TABLE]

is a Gaussian random variable with mean $\mathbb{E}[W(h)]=0$ and variance $\mathbb{E}[W(h)^{2}]=\left\|h\right\|_{H_{W}}^{2}$.

A.1 Cylindrical functionals and Malliavin Derivatives

Let $\mathcal{S}$ denote smooth cylindrical functionals of the form

[TABLE]

where $h_{i}\in H_{W}$ and $f\in C_{\mathrm{pol}}^{\infty}(\mathbb{R}^{n})$ (smooth functions with polynomial growth derivatives).

Definition A.4** (Malliavin derivative).**

For $F\in\mathcal{S}$, the Malliavin derivative is the $H_{W}$-valued random variable

[TABLE]

Beyond the score computation pursued in this paper, Malliavin derivatives have been used to approximate polynomial nonlinearities in nonlinear SPDEs via Wick–Malliavin expansions (Venturi et al., 2013).

Lemma A.5** (Malliavin derivative of the solution to the SPDE (1)).**

Let $\mathbb{D}^{1,2}$ be the Sobolev space defined as the completion of $\mathcal{S}$ under the norm

[TABLE]

The mild solution of the linear SPDE (1), i.e., (7), belongs to $\mathbb{D}^{1,2}(H)$ and its Malliavin derivative is given by

[TABLE]

Proof A.6**.**

Consider the perturbation $W^{\varepsilon}(s)=W(s)+\varepsilon\int_{0}^{s}\eta(\tau)\,d\tau$ for $\eta\in H_{W}$. The perturbed solution is

[TABLE]

Taking the Fréchet derivative

[TABLE]

By the Riesz representation theorem, $\mathcal{D}_{r}u(t)=S(t-r)Q^{1/2}\mathbf{1}_{[0,t]}(r)$. Finally, we verify $u(t)\in\mathbb{D}^{1,2}(H)$. To this end, we notice that

[TABLE]

by condition (6).

Hereafter, we characterise the Malliavin derivative of a function of the SPDE solution $u(t)$.

Theorem A.7** (Chain rule).**

Let $\phi\in C_{b}^{1}(H)$ with bounded Fréchet derivative. Then for $F=\phi(u(t))$

[TABLE]

Proof A.8**.**

By the chain rule for Fréchet derivatives

[TABLE]

Since $\mathcal{D}_{r}u(t)=S(t-r)Q^{1/2}\mathbf{1}_{[0,t]}(r)$ takes values in $\mathcal{L}_{\mathrm{HS}}(U,H)$, we need the adjoint action. For $v\in U$ we have

[TABLE]

Therefore $\mathcal{D}_{r}F=\left(Q^{1/2}\right)^{*}S(t-r)^{*}\nabla\phi(u(t))\mathbf{1}_{[0,t]}(r)$.

A.2 Skorokhod Integral

For $\eta\in H_{W}$, let us define the operator $\mathcal{D}u(t):H_{W}\to H$ as

[TABLE]

with adjoint $\mathcal{D}u(t)^{*}:H\to H_{W}$

[TABLE]

Definition A.9** (Skorokhod integral).**

The Skorokhod integral $\delta:\mathrm{Dom}(\delta)\subset L^{2}(\Omega;H_{W})\to L^{2}(\Omega)$ is the adjoint of $\mathcal{D}$

[TABLE]

for all $F\in\mathbb{D}^{1,2}$ and $v\in\mathrm{Dom}(\delta)$.

It can be shown that for deterministic $v\in H_{W}$, the Skorokhod integral reduces to the Itô integral

[TABLE]

Proposition A.10** (Integration by parts).**

Let $\gamma_{t}$ be the Malliavin covariance operator (11). For $F=\left\langle u(t),h\right\rangle_{H}$ with $h\in\mathrm{Ran}(\gamma_{t})$, and $v(r)=(Q^{1/2})^{*}S(t-r)^{*}\gamma_{t}^{-1}h\,\mathbf{1}_{[0,t]}(r)$ we have

[TABLE]

We have seen that the solution to the SPDE (1) is Gaussian with mean $S(t)u_{0}$ and covariance operator (11). This covariance operator satisfies the following recursion.

Lemma A.11** (Covariance recursion).**

For all $s,r\geq 0$

[TABLE]

Consequently, $S(t)\mathrm{Ran}(\gamma_{s})\subset\mathrm{Ran}(\gamma_{t+s})$ for all $t,s\geq 0$.

Proof A.12**.**

By direct computation,

[TABLE]

The range inclusion follows immediately.