The Rauch–Tung–Striebel recursion gives fixed-interval Gaussian smoothing laws

At time $T$, the conditioning data for the forward filter are exactly $y_{1:T}$. By [F3], the terminal filter kernel $N(m_{T\mathrm{\mid }T},P_{T\mathrm{\mid }T})$ is an everywhere conditional-law version of $X_T$ given these data, with $P_{T\mathrm{\mid }T}⪰0$. It is therefore the base case for the fixed-interval smoothing kernels in [F4].

Suppose the Gaussian kernel ${\psi }_{k+1\mathrm{\mid }T}(y,\cdot )$ is an RCD of $X_{k+1}$ given $Y_{1:T}=y$. Define ${\psi }_{k\mathrm{\mid }T}(y,A)=\int Q_k(y_{1:k},x,A){\psi }_{k+1\mathrm{\mid }T}(y,dx)$ for each Borel state set $A$. This is a probability kernel by [F5]–[F6]. To verify its RCD identity, define $\nu (D)=\int \text{\hspace{0.17em}⁣}\int 1_D(y,x){\psi }_{k+1\mathrm{\mid }T}(y,dx){\mathbb{P}}_Y(dy)$ for Borel $D$ in the product of data and state spaces. Measurability follows from [F6]; for pairwise disjoint $D_n$, monotone convergence gives countable additivity, and $\nu$ has total mass one. Its rectangle values equal those of the actual joint law of $(Y_{1:T},X_{k+1})$ by the RCD identity for ${\psi }_{k+1\mathrm{\mid }T}$. By [F8], the measures agree on every Borel set. Thus for any Borel data set $B$ and state set $A$, ${\int }_B{\psi }_{k\mathrm{\mid }T}(y,A){\mathbb{P}}_Y(dy)={\int }_{B\times {\mathbb{R}}^{d_{k+1}}}{Q}_k(y_{1:k},x,A){\mathbb{P}}_{(Y_{1:T},X_{k+1})}(dy,dx)=\mathbb{P}(X_k\in A,Y_{1:T}\in B)$, where the last equality is [F5] applied to the conditioning event $B\times {\mathbb{R}}^{d_{k+1}}$. Hence ${\psi }_{k\mathrm{\mid }T}$ is the required RCD.

Write $P^{-}=P_{k+1\mathrm{\mid }k}$ and $V_k=P_{k\mathrm{\mid }k}-C_k{P}^{-}{C}_k^T$ as in [F5]. By [F13], row substitution gives $L_k{U}_k=F_{k+1}{P}_{k\mathrm{\mid }k}$ and $L_k^T{C}_k^T=U_k$. Since $P^{-}=L_k{L}_k^T$, these equations yield $P^{-}{C}_k^T=F_{k+1}{P}_{k\mathrm{\mid }k}$ and hence $C_k=P_{k\mathrm{\mid }k}{F}_{k+1}^T(P^{-}{)}^{-1}$, exactly the gain in [F5]. On the product probability space carrying independent $U\sim N(0,V_k)$ and $W\sim N(m_{k+1\mathrm{\mid }T},P_{k+1\mathrm{\mid }T})$, the integral in step 2.1 is the law of $m_{k\mathrm{\mid }k}+C_k(W-m_{k+1\mathrm{\mid }k})+U$ by [F9]. The covariance matrices of $U$ and $W$ are positive semidefinite by [F5] and the induction hypothesis, so both Gaussian laws exist by [A1, F11]. The product coordinates are independent by [F12]; [F10] makes their stack jointly Gaussian, and its affine image is Gaussian by the projection characterization in [F11]. Independence makes the cross-covariance terms vanish, so its mean is $m_{k\mathrm{\mid }k}+C_k(m_{k+1\mathrm{\mid }T}-m_{k+1\mathrm{\mid }k})$ and its covariance is $V_k+C_k{P}_{k+1\mathrm{\mid }T}{C}_k^T=P_{k\mathrm{\mid }k}+C_k(P_{k+1\mathrm{\mid }T}-P^{-})C_k^T$. This is exactly the smoother update in the statement and algorithm. The first covariance expression is positive semidefinite, so the updated covariance is positive semidefinite as well.

Starting from step 1.1, apply steps 2.1–3.1 for $k=T-1,T-2,\dots ,0$. The induction constructs an everywhere smoothing kernel at each time and identifies its mean and covariance with the algorithm output. When $T=0$ only the base case is needed. Positive definiteness is used only for the forecast inverses in [F5] and the forward innovation solves in [F2]; no singular forecast is inverted.