Block D — Optional: coherent sets from data

Step 1 (the objective is a lagged second moment). By the definition of the inner product, $\ip{g,\K_\tau f}=\E[g(X_t)\,(\K_\tau f)(X_t)]$. By the definition of $\K_\tau$, $(\K_\tau f)(X_t)=\E[f(X_{t+\tau})\mid X_t]$. By the "taking out what is known" rule for conditional expectation ($g(X_t)$ is a function of the conditioning variable), $g(X_t)\,\E[f(X_{t+\tau})\mid X_t]=\E[g(X_t)f(X_{t+\tau})\mid X_t]$. By the tower property (Block 0, interlude), $\E\big[\E[g(X_t)f(X_{t+\tau})\mid X_t]\big]=\E[g(X_t)f(X_{t+\tau})]$. Chaining the four equalities gives the claim. Likewise $\|f\|^2=\E[f(X_t)^2]=\E[f(X_{t+\tau})^2]$ by stationarity, and $f\perp1$ means $\E[f(X_t)]=0$ by the definition of the inner product with the constant function.

Step 1 by direct integration. Write $\K_\tau$ as an integral operator, $\K_\tau f(x)=\int p_\tau(y\mid x)f(y)\,dy$ (the conditional expectation of a process with transition density $p_\tau$), and the inner product as an integral against $\rho_0$: $$\ip{g,\K_\tau f}=\int g(x)\Big[\int p_\tau(y\mid x)f(y)\,dy\Big]\rho_0(x)\,dx=\iint g(x)f(y)\,\big[\rho_0(x)\,p_\tau(y\mid x)\big]\,dx\,dy=\E[g(X_t)f(X_{t+\tau})].$$ Second equality: $g(x)\rho_0(x)$ is constant in $y$, so it moves inside the inner integral; the order of integration is immaterial because the double integral converges absolutely (Fubini; $f,g\in L^2(\rho_0)$ suffices by Cauchy–Schwarz). Third: $\rho_0(x)p_\tau(y\mid x)$ is the joint density of $(X_t,X_{t+\tau})$ (chain rule for densities: joint = marginal of the first times conditional of the second), and integrating a function of two variables against their joint density is its expectation. The tower property of the first version is the statement that the inner integral may be done first.

Step 2 (singular system). $\K_\tau$ is a contraction on $L^2(\rho_0)$ (Problem 7) and, under the paper's compactness assumption, has a singular system: orthonormal families $\{e_i\}$ and $\{h_i\}$ and numbers $1=\sigma_0\ge\sigma_1\ge\sigma_2\ge\cdots\ge0$ with $\K_\tau e_i=\sigma_ih_i$ and $\K_\tau^\dagger h_i=\sigma_ie_i$, so that $\K_\tau f=\sum_i\sigma_i\ip{e_i,f}h_i$ for every $f$. The top pair is $e_0=h_0=1$ with $\sigma_0=1$: $\K_\tau1=1$ because the conditional expectation of a constant is that constant, and $\K_\tau^\dagger1=1$ because stationarity makes $\rho_0$ invariant (Problem 7). No singular value exceeds the operator norm, which is $\le1$ for a contraction, so $\sigma_0=1$ is the largest.

Step 3 (expand and bound). Write $f=\sum_i\alpha_ie_i+f_\perp$ with $f_\perp$ orthogonal to all $e_i$ (so $\K_\tau f_\perp=0$ by the expansion above), and $g=\sum_i\beta_ih_i+g_\perp$ similarly. Then, by the expansion of $\K_\tau$ and orthonormality of $\{h_i\}$, $$\ip{g,\K_\tau f}=\sum_i\sigma_i\alpha_i\beta_i .$$ The constraints $f\perp1$ and $g\perp1$ say $\alpha_0=\ip{e_0,f}=0$ and $\beta_0=0$, so the sum starts at $i=1$. Now bound: $\sum_{i\ge1}\sigma_i\alpha_i\beta_i\le\sigma_1\sum_{i\ge1}|\alpha_i\beta_i|$ (ordering of the singular values), $\le\sigma_1\big(\sum\alpha_i^2\big)^{1/2}\big(\sum\beta_i^2\big)^{1/2}$ (Cauchy–Schwarz in $\ell^2$), $\le\sigma_1\|f\|\,\|g\|=\sigma_1$ (Bessel's inequality and the unit-norm constraints). Equality holds for $f=e_1$, $g=h_1$, which satisfy all constraints. So the maximum is $\sigma_1$, attained at the subdominant singular pair.

Which is which. $g$ is evaluated at $X_t$, the present, so $\operatorname{sign}g(X_t)$ is a statement about the current state that is correlated with $f(X_{t+\tau})$, a statement about the future: $g$ is the predictive function and its positive set is the earlier member $A$ of the coherent pair. $f$ is the retrospective function and its positive set is $B$. In the paper's Eq. (8), $\K_\tau u_i=\lambda_iv_i$, so the function that is fed to $\K_\tau$ lives on the later slice: $f=u_1$, $g=v_1$.

Setup. With $[\Sigma_0]_{ij}=\E[\phi_i(X_t)\phi_j(X_t)]$ and $[\Sigma_\tau]_{ij}=\E[\phi_i(X_t)\phi_j(X_{t+\tau})]$ (Problem 13(b)), the objective is, by linearity of expectation, $\E[g(X_t)f(X_{t+\tau})]=\sum_{ij}b_ia_j\E[\phi_i(X_t)\phi_j(X_{t+\tau})]=\vec b^{\,\top}\Sigma_\tau\vec a$, and the constraints are $\vec a^{\,\top}\Sigma_0\vec a=1$, $\vec b^{\,\top}\Sigma_0\vec b=1$ (same rule); the mean-zero constraints hold automatically because the features have mean zero. The Lagrangian is $$L(\vec a,\vec b,\alpha,\beta)=\vec b^{\,\top}\Sigma_\tau\vec a-\tfrac{\alpha}{2}\big(\vec a^{\,\top}\Sigma_0\vec a-1\big)-\tfrac{\beta}{2}\big(\vec b^{\,\top}\Sigma_0\vec b-1\big).$$

Step 1 (gradients). Rules: $\nabla_{\vec b}(\vec b^{\,\top}M\vec a)=M\vec a$ and $\nabla_{\vec a}(\vec b^{\,\top}M\vec a)=M^\top\vec b$ (gradient of a linear form); $\nabla_{\vec b}\tfrac12\vec b^{\,\top}S\vec b=S\vec b$ for symmetric $S$ (gradient of a quadratic form). Hence $$\nabla_{\vec b}L=\Sigma_\tau\vec a-\beta\,\Sigma_0\vec b,\qquad \nabla_{\vec a}L=\Sigma_\tau^\top\vec b-\alpha\,\Sigma_0\vec a .$$

Step 2 (gradients as expectations). By the definition of $\Sigma_\tau$ and linearity of expectation, $\Sigma_\tau\vec a=\E\big[\vec\phi(X_t)\,\vec a\cdot\vec\phi(X_{t+\tau})\big]=\E[\vec\phi(X_t)\,f(X_{t+\tau})]$; by the definition of $\Sigma_0$, $\Sigma_0\vec b=\E[\vec\phi(X_t)\,g(X_t)]$. By the definition of the transpose, $\Sigma_\tau^\top\vec b=\E[\vec\phi(X_{t+\tau})\,g(X_t)]$, and by stationarity $\Sigma_0\vec a=\E[\vec\phi(X_{t+\tau})f(X_{t+\tau})]$. Gradient ascent with step $\eta$ is therefore $$\Delta\vec b=\eta\,\E\big[\vec\phi(X_t)\big(f(X_{t+\tau})-\beta\,g(X_t)\big)\big],\qquad \Delta\vec a=\eta\,\E\big[\vec\phi(X_{t+\tau})\big(g(X_t)-\alpha\,f(X_{t+\tau})\big)\big].$$

Step 3 (one-sample rule). Each update is an expectation of a quantity computable from one pair of samples, so replacing the expectation by that quantity gives an unbiased stochastic gradient (its expectation is the gradient by construction): $$\Delta\vec b=\eta\,\vec\phi(x_t)\big[f(x_{t+\tau})-\beta\,g(x_t)\big],\qquad \Delta\vec a=\eta\,\vec\phi(x_{t+\tau})\big[g(x_t)-\alpha\,f(x_{t+\tau})\big].$$ Interpretation. The predictive neuron's weights $\vec b$ change by its own input $\vec\phi(x_t)$ times the retrospective partner's output $f(x_{t+\tau})$ one horizon later: a Hebbian term across time. The subtracted piece is its input times its own output, an anti-Hebbian (Oja-type) decay that holds $\E[g^2]$ at one. The retrospective neuron's rule is the mirror image, with its partner's output one horizon earlier. Each neuron needs its partner's output at a time offset by $\tau$, which in a circuit means a delay line or an eligibility trace, but it needs nothing else: no kernel, no model of the dynamics. The paper's Section 4 gives the fixed point of this rule, not the rule.

Step 4 (the multipliers). At a stationary point both gradients vanish. Multiply $\Sigma_\tau\vec a=\beta\Sigma_0\vec b$ on the left by $\vec b^{\,\top}$: $\vec b^{\,\top}\Sigma_\tau\vec a=\beta\,\vec b^{\,\top}\Sigma_0\vec b=\beta$ (constraint). Multiply $\Sigma_\tau^\top\vec b=\alpha\Sigma_0\vec a$ on the left by $\vec a^{\,\top}$: $\vec a^{\,\top}\Sigma_\tau^\top\vec b=\alpha$; and $\vec a^{\,\top}\Sigma_\tau^\top\vec b=(\vec b^{\,\top}\Sigma_\tau\vec a)^\top=\vec b^{\,\top}\Sigma_\tau\vec a$ (a scalar equals its transpose). So $\alpha=\beta=\vec b^{\,\top}\Sigma_\tau\vec a$, the achieved objective, as in Problem 14(a).

Step 5 (the fixed point is CCA). Writing the common multiplier as $\rho$, the stationarity conditions read $\Sigma_\tau\vec a=\rho\,\Sigma_0\vec b$ and $\Sigma_\tau^\top\vec b=\rho\,\Sigma_0\vec a$: exactly the CCA conditions of Problem 14(a) with $X=\vec\phi(x_t)$, $Y=\vec\phi(x_{t+\tau})$, whose solutions are the singular pairs of $\Sigma_0^{-1/2}\Sigma_\tau\Sigma_0^{-1/2}$ (14(b)(i)), equivalently the eigenvectors of the matrices of Eq. (24) (13(e)). The top non-trivial stationary point is the global maximum, so the rule's target is the paper's filter. Medium on the dynamics of the rule itself: plain gradient ascent with the multipliers held at their stationary values does not by itself keep the iterates on the constraint set; the standard remedy (Oja 1982, and its CCA generalisations) replaces $\alpha,\beta$ by the current output products, which changes the transient but not the fixed points.

Complete basis. Let $\{\psi_k\}$ be a complete orthonormal set in $L^2(\rho_0)$. Every $f\in L^2(\rho_0)$ is $f=\sum_ka_k\psi_k$ with $a_k=\ip{\psi_k,f}$ (completeness), and $\|f\|^2=\sum_ka_k^2$ (Parseval). By linearity and continuity of $\K_\tau$ and of the inner product, $\ip{g,\K_\tau f}=\sum_{jk}b_ja_k\ip{\psi_j,\K_\tau\psi_k}=\vec b^{\,\top}\Sigma_\tau\vec a$ with $\Sigma_0=I$. So the map $f\mapsto\vec a$ is an isometric bijection between $L^2(\rho_0)$ and $\ell^2$ carrying the objective and constraints of Q1(i) to those of Q1(ii): the two problems have the same value and corresponding maximisers. Learning in a complete basis finds the exact singular pair.

Finite basis. For $f,g\in F$: $\ip{g,\K_\tau f}=\ip{P_Fg,\K_\tau f}$ ($g\in F$ so $P_Fg=g$) $=\ip{g,P_F\K_\tau f}$ ($P_F$ is self-adjoint, being an orthogonal projection) $=\ip{g,Cf}$ (definition of $C$). So the in-subspace problem is the singular-value problem of $C$ on $F\ominus\{1\}$, its value is $\sigma_1(C)$, and $\sigma_1(C)\le\sigma_1(\K_\tau)$ by Problem 15(a).

The example. For this OU process the stationary law is $N(0,1)$ and the Hermite polynomials are eigenfunctions, $\K_\tau\mathrm{He}_n=e^{-n\tau}\mathrm{He}_n$, orthogonal under $N(0,1)$ with $\|\mathrm{He}_n\|^2=n!$ (Problem 11). The process is reversible, so singular functions coincide with eigenfunctions and $\sigma_n=e^{-n\tau}$: the true subdominant pair is $f=g=\mathrm{He}_1=x$, value $e^{-\tau}$. With $F=\operatorname{span}\{\mathrm{He}_2\}$ the only admissible unit functions are $f,g=\pm\mathrm{He}_2/\sqrt2$ (mean zero holds because $\E[\mathrm{He}_2]=\E[X^2]-1=0$), and the in-subspace objective is $\ip{\mathrm{He}_2,\K_\tau\mathrm{He}_2}/2=e^{-2\tau}\|\mathrm{He}_2\|^2/2=e^{-2\tau}$ (eigen-relation, then the norm). Meanwhile the projection of the true pair onto $F$ is $\ip{\mathrm{He}_2,x}\,\mathrm{He}_2/2=0$, because $\ip{\mathrm{He}_2,x}=\E[X^3]-\E[X]=0$ for a standard normal. So the learner in this basis converges to a nonzero function whose "coherence" $e^{-2\tau}$ is genuine but which is the third singular function, not any approximation of the second; the projection of the second is zero. Which singular function a finite basis delivers depends on the basis, and no quantity computable inside $F$ reports the gap $e^{-\tau}-e^{-2\tau}$ to the true optimum (Problem 15(e)).

Empirical problem. Replacing each expectation in Q1(i) by a sample average (the plug-in estimator) gives $$\hat J(f,g)=\tfrac1N\sum_ng(x_n)f(y_n),\qquad \tfrac1N\sum_ng(x_n)^2=\tfrac1N\sum_nf(y_n)^2=1,\qquad \tfrac1N\sum_ng(x_n)=\tfrac1N\sum_nf(y_n)=0 .$$ Each of these depends on $g$ only through the vector of values $(g(x_n))_n$ and on $f$ only through $(f(y_n))_n$.

Degeneracy. By Cauchy–Schwarz in $\R^N$, $\hat J\le\big(\tfrac1N\sum g(x_n)^2\big)^{1/2}\big(\tfrac1N\sum f(y_n)^2\big)^{1/2}=1$. If the class can interpolate, pick any signs $s_n\in\{\pm1\}$ with $\sum_ns_n=0$ ($N$ even) and choose $g$ with $g(x_n)=s_n$ and $f$ with $f(y_n)=s_n$ (possible because the $x_n$ are distinct and the $y_n$ are distinct). Then the mean constraints hold, the second-moment constraints hold ($s_n^2=1$), and $\hat J=\tfrac1N\sum_ns_n^2=1$: the bound is attained. The signs were arbitrary, so the "perfectly coherent pair" found is one of $\binom{N}{N/2}$ equally good ones, none of which has anything to do with $\K_\tau$, whose true $\sigma_1$ is strictly less than $1$ for a noisy system (Block B, Q2(vi)). Interpolation is available to any class of dimension $\ge N$, so the estimator is meaningful only when the dictionary size $d$ is small compared with $N$, or an explicit penalty is added. The finite basis Section 4 starts from is therefore a regulariser, not only a limitation, and this is a legitimate reason to introduce it early.

Hierarchy. Let the upstream layer have learned filters $\vec b_1,\dots,\vec b_m$ and emit the ON/OFF pair $z_i^\pm=[\pm\vec b_i\cdot\vec\phi]_+$ (the paper's proposal that neurons rectify the positive and negative parts). A downstream neuron with these $2m$ features has function class $F_2=\operatorname{span}\{z_i^\pm\}$. It contains every upstream function: $\vec b_i\cdot\vec\phi=z_i^+-z_i^-$ (the identity $z=[z]_+-[-z]_+$ for any real $z$), so $\operatorname{span}\{\vec b_i\cdot\vec\phi\}\subseteq F_2$. It is strictly larger: $z_i^++z_i^-=|\vec b_i\cdot\vec\phi|$ is in $F_2$ and is not a linear function of $\vec\phi$. So each layer enlarges the dictionary with nonlinear functions of the input, and the downstream neuron runs the same rule as Q1(ii) on its own inputs: the approximation problem is attacked by depth, the learning problem stays local. Medium on how far this goes: the upstream filters were chosen to maximise their objective, not to be useful features for the layer above, so nothing guarantees that $F_2$ approaches the true nonlinear singular functions; the paper claims composability, not convergence.

Take $f=\vec c\cdot\vec\phi\in F$. By linearity of $\K_\tau$ (Problem 1), $\K_\tau f=\sum_jc_j\,\K_\tau\phi_j$, and each $\K_\tau\phi_j(x)=\E[\phi_j(X_{t+\tau})\mid X_t=x]$ is some new function of $x$ with no reason to be a combination of the $\phi_i$. Example: on a nonlinear system $\K_\tau x$ is a nonlinear function of $x$, so with the single feature $\phi(x)=x$ the image leaves $F$.

Step 1 (closest point = orthogonal projection). The closest element $g$ of a subspace $F$ to a given function $h$ is characterised by the residual $h-g$ being orthogonal to every element of $F$, equivalently to every basis function: $\ip{\phi_i,h-g}=0$ for all $i$. (Proof: for any $\delta\in F$, $\|h-g-\delta\|^2=\|h-g\|^2-2\ip{h-g,\delta}+\|\delta\|^2$, which is minimised at $\delta=0$ for all $\delta$ iff $\ip{h-g,\delta}=0$ for all $\delta\in F$.) This orthogonality condition is the Galerkin condition.

Step 2 (apply it to $h=\K_\tau f$). Seek $g=\vec c\,'\cdot\vec\phi$ with $\ip{\phi_i,\K_\tau f-g}=0$ for every $i$. Expanding both terms by linearity of the inner product, $$\sum_jc_j\ip{\phi_i,\K_\tau\phi_j}=\sum_jc'_j\ip{\phi_i,\phi_j}\qquad\text{for all }i,$$ which in matrix form is $\Sigma_\tau\vec c=\Sigma_0\,\vec c\,'$.

Step 3 (solve). Linear independence of the $\phi_i$ makes $\Sigma_0$ (a Gram matrix) invertible, so $\vec c\,'=\Sigma_0^{-1}\Sigma_\tau\,\vec c$. The map $\vec c\mapsto\vec c\,'$ is the matrix $K_\tau=\Sigma_0^{-1}\Sigma_\tau$, Eq. (22) (Supplement Eq. 51). Its $(i,j)$ entry answers: how much of $\phi_i$ appears in the projection of $\K_\tau\phi_j$.

The name comes from numerical PDEs, where Galerkin's method means: look for the solution in a finite-dimensional subspace and impose the equation only against test functions from that subspace. In the Koopman literature the same construction is called extended dynamic mode decomposition (EDMD; Williams, Kevrekidis & Rowley 2015, Klus et al. 2018), which are the paper's references for Eq. (22).

$[\Sigma_0]_{ij}=\ip{\phi_i,\phi_j}=\E[\phi_i(X_t)\phi_j(X_t)]$ is the definition. For $\Sigma_\tau$, insert the definition of $\K_\tau$ and use the two rules of Problem 7(d): $$[\Sigma_\tau]_{ij}=\E\big[\phi_i(X_t)\,\K_\tau\phi_j(X_t)\big]=\E\Big[\phi_i(X_t)\,\E[\phi_j(X_{t+\tau})\mid X_t]\Big]=\E\Big[\E\big[\phi_i(X_t)\phi_j(X_{t+\tau})\mid X_t\big]\Big]=\E[\phi_i(X_t)\phi_j(X_{t+\tau})],$$ moving $\phi_i(X_t)$ inside the conditional expectation (it is a function of the conditioning variable) and then collapsing the nesting by the tower property. So $\Sigma_\tau$ is the lagged cross-covariance of the features, and $\Sigma_0$ their instantaneous covariance. Replacing each expectation by the time average $\tfrac1S\sum_s$ over a trajectory gives Eq. (23).

By direct integration. With $\K_\tau\phi_j(x)=\int p_\tau(y\mid x)\phi_j(y)\,dy$ (definition of $\K_\tau$ through the transition density) and the stationary inner product written as an integral, $$[\Sigma_\tau]_{ij}=\int\phi_i(x)\Big[\int p_\tau(y\mid x)\phi_j(y)\,dy\Big]\rho_0(x)\,dx=\iint\phi_i(x)\phi_j(y)\,\rho_0(x)p_\tau(y\mid x)\,dx\,dy=\E[\phi_i(X_t)\phi_j(X_{t+\tau})],$$ moving $\phi_i(x)\rho_0(x)$ inside the $y$-integral (it is constant in $y$; Fubini for the order) and recognising $\rho_0(x)p_\tau(y\mid x)$ as the joint density of $(X_t,X_{t+\tau})$ (joint = marginal times conditional). The Galerkin matrix is a lagged moment because the Koopman kernel times the stationary density is the two-time joint density.

The estimate uses: stationarity, so that $X_{t_s}$ and the pairs $(X_{t_s},X_{t_s+\tau})$ have the same law for every $s$ and the averages estimate a single $\Sigma_0$, $\Sigma_\tau$; and ergodicity, so that a time average over one trajectory converges to the ensemble expectation. This is why the paper discusses local stationarity and exponential forgetting immediately after Eq. (23). It is also the reason for choosing the $\rho_0$ inner product rather than, say, Lebesgue: it is the one whose Gram matrices a neuron can estimate from the lagged correlations of its own inputs.

$[\Sigma_0]_{ij}=\E[\mathbf 1_{\{i\}}(X_t)\mathbf 1_{\{j\}}(X_t)]=P(X_t=i,X_t=j)=\pi_i\delta_{ij}$, so $\Sigma_0=\diag(\pi)$. $[\Sigma_\tau]_{ij}=\E[\mathbf 1_{\{i\}}(X_t)\mathbf 1_{\{j\}}(X_{t+1})]=P(X_t=i,X_{t+1}=j)=\pi_iP_{ij}$, so $\Sigma_\tau=\diag(\pi)P$. Hence $K_\tau=\diag(\pi)^{-1}\diag(\pi)P=P$, the Koopman matrix of Problem 2(a). No approximation occurs because the indicators span every function on $N$ states, so $\K_\tau f$ never leaves $F$ and the projection is the identity. Note also that in coefficient coordinates the inner product is $\vec c_f^{\,\top}\Sigma_0\vec c_g=\sum_i\pi_ic_{f,i}c_{g,i}=\ip{f,g}_\pi$: the $\pi$-weighted inner product of Block A is what Supplement Eq. (50) becomes here.

For a continuous process and indicators of a partition $\{B_i\}$, the same computation gives $[\Sigma_0]_{ij}=\rho_0(B_i)\delta_{ij}$ and $[\Sigma_\tau]_{ij}=P(X_t\in B_i,X_{t+\tau}\in B_j)$, so $K_\tau$ is the matrix of transition probabilities between cells, $P(X_{t+\tau}\in B_j\mid X_t\in B_i)$: Ulam's method. Now there is an approximation (functions constant on cells are not invariant under $\K_\tau$), and the numerator of the paper's Eq. (4) is exactly one entry of this matrix.

For OU, $x_{t+\tau}=e^{A\tau}x_t+\eta$ with $\eta$ zero-mean Gaussian independent of $x_t$ (Problem 10a), so $\K_\tau(\vec c\cdot x)=\E[\vec c\cdot(e^{A\tau}x+\eta)]=\vec c\cdot e^{A\tau}x=(e^{A^\top\tau}\vec c)\cdot x$. The image is again linear, so $F=\{\vec c\cdot x\}$ is invariant, the Galerkin condition is satisfied with zero residual, and the representation is exact. Reading off the coefficient map: $K_\tau=e^{A^\top\tau}$.

From Eq. (22): $\Sigma_0=\E[x_tx_t^\top]=\Sigma$ and $\Sigma_\tau=\E[x_tx_{t+\tau}^\top]=\E[x_t(e^{A\tau}x_t+\eta)^\top]=\Sigma e^{A^\top\tau}$, so $K_\tau=\Sigma^{-1}\Sigma e^{A^\top\tau}=e^{A^\top\tau}$. The two agree.

The transpose is not a convention clash. Points move by $x\mapsto e^{A\tau}x$; a linear observable is a covector $\vec c$, and covectors transform by the transpose (contragrediently): $\vec c\cdot(Mx)=(M^\top\vec c)\cdot x$. This is the same transpose that makes the paper's Eq. (18) read $e^{A^\top\tau}\Sigma^{-1}e^{A\tau}\Sigma$ for the $v$'s rather than $e^{A\tau}\cdots$, and it is behind the naming trap noted on the overview page.

By bilinearity, $\ip{f,g}=\sum_{ij}c_{f,i}c_{g,j}\ip{\phi_i,\phi_j}=\vec c_f^{\,\top}\Sigma_0\vec c_g$. So $\R^d$ carries the metric $\Sigma_0$, and the adjoint of a matrix $K$ is defined by $\vec c_g^{\,\top}\Sigma_0(K\vec c_f)=(K^\dagger\vec c_g)^\top\Sigma_0\vec c_f$ for all $\vec c_f,\vec c_g$. The left side is $\vec c_g^{\,\top}\Sigma_0K\vec c_f$ and the right side is $\vec c_g^{\,\top}K^{\dagger\top}\Sigma_0\vec c_f$; equating the matrices gives $\Sigma_0K=K^{\dagger\top}\Sigma_0$, i.e. $K^\dagger=\Sigma_0^{-1}K^\top\Sigma_0$. This is Problem 3(a) with $\diag(\pi)$ replaced by the Gram matrix $\Sigma_0$; (c) shows they coincide for indicator features. The paper writes this as "a direct consequence of Eq. (5)" just above Eq. (24).

With $K_\tau=\Sigma_0^{-1}\Sigma_\tau$ and $\Sigma_\tau^\top=\Sigma_{-\tau}$ (the $(i,j)$ entry $\E[\phi_i(X_t)\phi_j(X_{t+\tau})]$ transposed is $\E[\phi_j(X_t)\phi_i(X_{t+\tau})]=\E[\phi_i(X_{t-\tau})\phi_j(X_t)]$ by stationarity): $K_\tau^\dagger=\Sigma_0^{-1}\Sigma_\tau^\top\Sigma_0^{-\top}\Sigma_0=\Sigma_0^{-1}\Sigma_{-\tau}$, using $\Sigma_0^\top=\Sigma_0$. Then $$F_\tau=K_\tau K_\tau^\dagger=\Sigma_0^{-1}\Sigma_\tau\Sigma_0^{-1}\Sigma_{-\tau},\qquad B_\tau=K_\tau^\dagger K_\tau=\Sigma_0^{-1}\Sigma_{-\tau}\Sigma_0^{-1}\Sigma_\tau,$$ which is Eq. (24). Their eigenvectors are the coefficient vectors of the singular functions (Eq. 25), because singular vectors of $K$ are eigenvectors of $KK^\dagger$ and $K^\dagger K$, exactly as in Problem 4. Problem 14 evaluates these for OU and identifies them with CCA.

$\K_\tau x^2=\E[(e^{-\tau}x+\eta)^2]=e^{-2\tau}x^2+(1-e^{-2\tau})$: a quadratic plus a constant. Also $\K_\tau x=e^{-\tau}x$ and $\K_\tau1=1$. Needed Gaussian moments: $\E[x^2]=1$, $\E[x^3]=0$, $\E[x^4]=3$.

(i) Features $\{x,x^2\}$. $\Sigma_0=\begin{pmatrix}\E[x^2]&\E[x^3]\\ \E[x^3]&\E[x^4]\end{pmatrix}=\diag(1,3)$. For $\Sigma_\tau$ use $[\Sigma_\tau]_{ij}=\E[\phi_i(x)\,\K_\tau\phi_j(x)]$: $\E[x\cdot e^{-\tau}x]=e^{-\tau}$; $\E[x\,(e^{-2\tau}x^2+1-e^{-2\tau})]=0$; $\E[x^2\cdot e^{-\tau}x]=0$; $\E[x^2(e^{-2\tau}x^2+1-e^{-2\tau})]=3e^{-2\tau}+1-e^{-2\tau}=1+2e^{-2\tau}$. So $$K_\tau=\Sigma_0^{-1}\Sigma_\tau=\diag\!\Big(e^{-\tau},\ \tfrac{1+2e^{-2\tau}}{3}\Big).$$ The $x$ entry is exact. The $x^2$ entry is $\tfrac13+\tfrac23e^{-2\tau}$ instead of $e^{-2\tau}$: at $\tau=0.5$, $0.579$ versus $0.368$. The constant $(1-e^{-2\tau})$ in $\K_\tau x^2$ is not in $F$, and it is not orthogonal to $x^2$ either ($\ip{1,x^2}=1\ne0$), so the projection smears part of it onto the $x^2$ direction: coefficient $\ip{x^2,\,1-e^{-2\tau}}/\ip{x^2,x^2}=(1-e^{-2\tau})/3$, which is exactly the excess. A Galerkin matrix is only as good as the invariance of the feature span, and the error lands on whichever features are correlated with what was dropped.

(ii) Features $\{1,x,x^2\}$. Now $\K_\tau$ maps $F$ into itself (polynomials of degree $\le2$ are invariant, Problem 11), so the matrix is exact. Reading coefficients off $\K_\tau1=1$, $\K_\tau x=e^{-\tau}x$, $\K_\tau x^2=(1-e^{-2\tau})\cdot1+e^{-2\tau}x^2$, with column $j$ holding the coefficients of $\K_\tau\phi_j$: $$K_\tau=\begin{pmatrix}1&0&1-e^{-2\tau}\\0&e^{-\tau}&0\\0&0&e^{-2\tau}\end{pmatrix},$$ upper triangular with eigenvalues $1,e^{-\tau},e^{-2\tau}$. Eq. (22) reproduces this (the reader may check: $\Sigma_0=\begin{pmatrix}1&0&1\\0&1&0\\1&0&3\end{pmatrix}$, $\Sigma_\tau=\begin{pmatrix}1&0&1\\0&e^{-\tau}&0\\1&0&1+2e^{-2\tau}\end{pmatrix}$). The eigenvalues are those of the Hermite eigenfunctions $\mathrm{He}_0,\mathrm{He}_1,\mathrm{He}_2$ of Problem 11 (with $a=1$, $\sigma=1$), and the eigenvectors are their coefficient vectors: $(1,0,0)$, $(0,1,0)$ and $(-1,0,1)$, i.e. $x^2-1=\mathrm{He}_2$. Including the constant feature is what makes the trivial pair $u_0=v_0=\mathbf 1$ of Problem 7(b) representable; the paper's Eq. (26) omits it, which is harmless for zero-mean linear features because $\ip{1,x}=0$, but would not be for $x^2$.

(iii) Running the script ($\Delta t=0.01$, $T=2\times10^4$, $\tau=0.5$) gives, for $\{x,x^2\}$, diagonal entries $0.601$ and $0.572$ against the analytic $0.607$ and $0.579$; for $\{1,x,x^2\}$, the $(1,3)$ entry $0.642$ against $0.632$ and eigenvalues $1.000,0.601,0.355$ against $1,0.607,0.368$. The residual gap is Euler–Maruyama bias at this $\Delta t$ (the discrete map has effective decay $1-\Delta t$ per step rather than $e^{-\Delta t}$), as in Problem 14(d).

Maximising correlation is maximising covariance at fixed variances, hence the constrained form. Lagrangian $\mathcal L=a^\top\Sigma_{XY}b-\tfrac\lambda2(a^\top\Sigma_{XX}a-1)-\tfrac\mu2(b^\top\Sigma_{YY}b-1)$. Setting $\partial_a\mathcal L=0$ and $\partial_b\mathcal L=0$: $$\Sigma_{XY}b=\lambda\,\Sigma_{XX}a,\qquad\Sigma_{YX}a=\mu\,\Sigma_{YY}b .$$ Left-multiply the first by $a^\top$ and the second by $b^\top$ and use the constraints: $a^\top\Sigma_{XY}b=\lambda$ and $b^\top\Sigma_{YX}a=\mu$. These are the same scalar (transpose of each other), so $\lambda=\mu=\rho$, the correlation achieved at the optimum. From the second equation $b=\rho^{-1}\Sigma_{YY}^{-1}\Sigma_{YX}a$; substituting into the first, $\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}a=\rho^2\Sigma_{XX}a$, i.e. $\Sigma_{XX}^{-1}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\,a=\rho^2a$, and symmetrically $\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}\,b=\rho^2b$. The top eigenvector gives the first canonical pair; the remaining eigenvectors give the successive pairs, and the constraint that later pairs be uncorrelated with earlier ones is the $\Sigma_{XX}$-orthogonality of eigenvectors of this (self-adjoint in the $\Sigma_{XX}$ metric, cf. Problem 3) matrix.

Whitened form. Put $\tilde a=\Sigma_{XX}^{1/2}a$, $\tilde b=\Sigma_{YY}^{1/2}b$ and $\tilde C=\Sigma_{XX}^{-1/2}\Sigma_{XY}\Sigma_{YY}^{-1/2}$. The constraints become $\|\tilde a\|=\|\tilde b\|=1$, the objective $\tilde a^\top\tilde C\tilde b$, and the stationarity conditions $\tilde C\tilde b=\rho\tilde a$, $\tilde C^\top\tilde a=\rho\tilde b$: the singular value equations of $\tilde C$. So the canonical correlations are the singular values of $\tilde C$ and the canonical directions its singular vectors, unwhitened by $\Sigma^{-1/2}$. This is the same isometry trick as Problem 3(c): whitening turns a weighted problem into a Euclidean one.

Invariance. If $X\mapsto GX$ then $a^\top X=(G^{-\top}a)^\top(GX)$, so the set of achievable projections $\{a^\top X\}$ is unchanged and so are the optima and the $\rho_i$; only the coordinates of $a$ change. CCA therefore depends on $\operatorname{span}\{X_1,\dots,X_p\}$ alone. For Problem 13 this says the Galerkin matrices depend on the feature span $F$, not on the basis chosen for it, exactly as the derivation there (a projection onto $F$) implies; recombining the features by an invertible matrix changes $K_\tau$ by a similarity and leaves the singular functions untouched.

(i) $\E[\vec\psi\vec\psi^\top]=\Sigma_0^{-1/2}\E[\vec\phi\vec\phi^\top]\Sigma_0^{-1/2}=I$, so $\ip{\psi_i,\psi_j}=\delta_{ij}$: an orthonormal basis of $F$. The matrix of $C$ in an orthonormal basis has entries $\ip{\psi_i,C\psi_j}=\ip{\psi_i,\K_\tau\psi_j}$ (the projection drops because $\psi_i\in F$) $=\E[\psi_i(x_t)\psi_j(x_{t+\tau})]$ by the lagged-moment identity of Problem 13(b) (tower property, or direct integration), and $\E[\vec\psi(x_t)\vec\psi(x_{t+\tau})^\top]=\Sigma_0^{-1/2}\Sigma_\tau\Sigma_0^{-1/2}=\tilde C$. So canonical correlations are the singular values of the compressed Koopman operator, and Problem 15(b) is the statement that CCA solves the two-sided restricted problem. Two proofs of $\rho_i\le1$: statistically, a correlation coefficient is at most one (Cauchy–Schwarz); operator-theoretically, $\sigma_i(C)\le\sigma_i(\K_\tau)\le1$ by Problem 15(a) and the contraction bound of Problem 7. That the two arguments agree is the content of the identification.

(ii) Detailed balance says $(x_t,x_{t+\tau})$ and $(x_{t+\tau},x_t)$ have the same joint law (Problem 7d), so $\E[\phi_i(x_t)\phi_j(x_{t+\tau})]=\E[\phi_j(x_t)\phi_i(x_{t+\tau})]$, i.e. $\Sigma_\tau=\Sigma_\tau^\top$ and $\tilde C$ is symmetric. The SVD of a symmetric matrix is its eigendecomposition up to signs ($\tilde C=Q\Lambda Q^\top$ with $\sigma_i=|\lambda_i|$, $\tilde b_i=\operatorname{sign}(\lambda_i)\tilde a_i$), so $a_i=\pm b_i$: a reversible process has no distinction between predictive and retrospective filters, which is Problem 3(c)/9(c) again. (This uses that the features are functions of a Markov state whose time reversal is the same state. For lag vectors of a reversible scalar signal, time reversal also reverses the order of the lags, so $\Sigma_\tau$ is not symmetric; instead $\Sigma_\tau^\top=J\Sigma_\tau J$ with $J$ the lag-reversal permutation, and the future filter is the lag-reversed past filter. See Problem 16(d).) The stationarity condition $\Sigma_{XY}b=\rho\Sigma_{XX}a$ with $a=b$ reads $\Sigma_\tau a=\rho\Sigma_0a$, the TICA generalised eigenproblem, whose solutions are the Galerkin eigenvectors of $K_\tau=\Sigma_0^{-1}\Sigma_\tau$ themselves. For an irreversible process $\Sigma_\tau\ne\Sigma_\tau^\top$, eigenvectors of $K_\tau$ and singular vectors of $\tilde C$ part company, and one must use the latter: this is why the paper (and Problem 4) insists on singular rather than eigen-functions.

(iii) For a Gaussian vector with covariance $\Sigma$, $h=\tfrac12\log\det(2\pi e\Sigma)$, so $I(X;Y)=h(X)+h(Y)-h(X,Y)=\tfrac12\log\dfrac{\det\Sigma_{XX}\det\Sigma_{YY}}{\det\begin{pmatrix}\Sigma_{XX}&\Sigma_{XY}\\\Sigma_{YX}&\Sigma_{YY}\end{pmatrix}}$. Mutual information is invariant under invertible maps of each argument, so whiten: the joint covariance becomes $\begin{pmatrix}I&\tilde C\\\tilde C^\top&I\end{pmatrix}$, whose determinant by the Schur complement is $\det(I-\tilde C^\top\tilde C)=\prod_i(1-\rho_i^2)$, while the marginal determinants are $1$. Hence $I(X;Y)=-\tfrac12\sum_i\log(1-\rho_i^2)$: each canonical pair contributes independently, in decreasing order of $\rho_i$. For the OU process the lag vector is Gaussian, so this is exactly the predictive information between $x_t$ and $x_{t+\tau}$, and the $\lambda_i$ of Eq. (18) are its coordinates. With the analytic values of part (d), $\rho^2=(0.434,0.115)$: $I=-\tfrac12[\log0.566+\log0.885]=0.35$ nats, of which the subdominant pair carries $0.28$. The paper's "coherent set" filter is, for a Gaussian input, the single linear projection carrying the most predictive information; that is the sense in which the coherent-set, CCA and information-bottleneck pictures coincide, and the sense in which they do not: the information-theoretic statement needs Gaussianity, the CCA statement needs only second moments, and the coherent-set statement is about level sets of the singular function rather than the function itself.

Analytic values for reference: $\Sigma=\begin{pmatrix}0.5&0.1667\\0.1667&0.3333\end{pmatrix}$ and $A\Sigma=\begin{pmatrix}-0.5&-0.1667\\0.1667&-0.5\end{pmatrix}$, not symmetric, so the process is irreversible and $F_\tau\ne B_\tau$ is expected. The estimator uses $k=\tau/\Delta t=50$ samples of lag; $\Sigma_0$ is estimated symmetrically from both ends of each lagged pair, $\Sigma_\tau$ from $X_0^\top X_1/n$.

import numpy as np
from scipy.linalg import expm, solve_sylvester
rng = np.random.default_rng(0)
A = np.array([[-1., 0.], [1., -2.]]); D = np.eye(2)
dt, T, tau = 0.01, 4000., 0.5; n = int(T/dt); k = round(tau/dt)
x = np.zeros((n, 2)); noise = rng.standard_normal((n, 2))*np.sqrt(dt)
for t in range(n-1): x[t+1] = x[t] + A @ x[t]*dt + noise[t]
x = x[n//10:]; X0, X1 = x[:-k], x[k:]
S0 = (X0.T@X0 + X1.T@X1)/(2*len(X0)); St = X0.T@X1/len(X0); Smt = St.T
F_hat = np.linalg.solve(S0, St) @ np.linalg.solve(S0, Smt)     # Eq. (24)
B_hat = np.linalg.solve(S0, Smt) @ np.linalg.solve(S0, St)

From Problem 10(a), $\Sigma_\tau=\E[x_tx_{t+\tau}^\top]=\Sigma e^{A^\top\tau}$ and $\Sigma_{-\tau}=e^{A\tau}\Sigma$, $\Sigma_0=\Sigma$. Substituting into Eq. (24): $$F_\tau=\Sigma^{-1}\Sigma e^{A^\top\tau}\,\Sigma^{-1}e^{A\tau}\Sigma=e^{A^\top\tau}\Sigma^{-1}e^{A\tau}\Sigma=M^\top,\qquad B_\tau=\Sigma^{-1}e^{A\tau}\Sigma\,\Sigma^{-1}\Sigma e^{A^\top\tau}=\Sigma^{-1}e^{A\tau}\Sigma e^{A^\top\tau}=N^\top,$$ exactly the two matrices of Eq. (18). (This is the Supplement Sec. 6 claim; note it needs $\Sigma_0=\Sigma$, i.e. a stationary process with a genuine covariance — the same caveat as Problem 12(d).) Numerically, with $\tau=0.5$:

Eigenvalues agree to the second decimal (the residual bias is Euler–Maruyama discretisation at $\Delta t=0.01$, which perturbs the effective $\Sigma$; halving $\Delta t$ halves it), eigenvectors to three decimals. As required, $F_\tau$ and $B_\tau$ share eigenvalues but have different eigenvectors: $v_1=(0.93,0.37)$ versus $u_1=(0.75,0.66)$, an angle of about $17^\circ$. A "predictive" unit and a "retrospective" unit for this input would carry different filters, and the difference is entirely due to $c\ne0$ (Problem 9d).

CCA between $X=x_t$ and $Y=x_{t+\tau}$ with $\Sigma_{XX}=\Sigma_{YY}=\Sigma_0$ and $\Sigma_{XY}=\Sigma_\tau$ seeks $a,b$ maximising $\operatorname{corr}(a^\top X,b^\top Y)$; the stationary conditions derived in (a) give the eigenproblems $\Sigma_{XX}^{-1}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\,a=\rho^2a$ and $\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}\,b=\rho^2b$, which are precisely $F_\tau a=\rho^2a$ and $B_\tau b=\rho^2b$ with $\rho$ the canonical correlations. Equivalently, the $\rho_i$ are the singular values of the whitened cross-covariance $\Sigma_0^{-1/2}\Sigma_\tau\Sigma_0^{-1/2}$, whose right/left singular vectors are $\Sigma_0^{1/2}a$ and $\Sigma_0^{1/2}b$. Numerically: canonical correlations squared $=(0.438,0.123)=$ eigenvalues of $F_\tau$.

So on linear features the whole construction is past–future CCA (Koltai et al. 2018; Wu & Noé 2020, where $\sum\rho_i^2$ is the VAMP-2 score), and $\Sigma_0^{-1/2}\Sigma_\tau\Sigma_0^{-1/2}$ is $\tilde M^\top$ of Problem 10(d): the whitened transition matrix is the canonical-correlation matrix. The paper's contribution beyond CCA is the coherent-set reading of the sign and the identification of the two canonical bases with predictive and retrospective units; the neural claim is that a cell can run this eigenproblem on its own lagged inputs, for which Lipshutz et al. (2021) give a local Hebbian/anti-Hebbian circuit. High for the CCA identity; the biological plausibility of the online estimator is outside this problem set.

For $S\subset F$, $\|Cf\|=\|P_F\K f\|\le\|\K f\|$ because an orthogonal projection is a contraction. Hence $$\sigma_i(C)=\max_{\substack{S\subset F\\\dim S=i}}\;\min_{f\in S}\frac{\|P_F\K f\|}{\|f\|}\;\le\;\max_{\substack{S\subset F\\\dim S=i}}\;\min_{f\in S}\frac{\|\K f\|}{\|f\|}\;\le\;\max_{\dim S=i}\;\min_{f\in S}\frac{\|\K f\|}{\|f\|}=\sigma_i(\K),$$ the first inequality by the contraction, the second by enlarging the set of subspaces maximised over. The same argument on $C^\dagger=P_F\K^\dagger|_F$ (the $\Sigma_0$-adjoint of Problem 13e) gives nothing new since $C$ and $C^\dagger$ share singular values.

Eigenvalues have no min–max characterisation unless the operator is self-adjoint (or normal). All one can say is that eigenvalues of $C$ lie in the numerical range $\{\ip{f,\K f}/\|f\|^2\}$ of $\K$, which for a non-normal operator can be much larger than the convex hull of its spectrum; that is spectral pollution. The 1-D OU of Problem 13(f) is reversible, so $\K$ is self-adjoint and positive, its singular values are its eigenvalues $e^{-n\tau}$, and the compressed value $0.579$ found there does obey (a): it sits below $\sigma_1=1$, and indeed between the true eigenvalues $0.607$ and $0.368$ whose eigenfunctions $x^2$ mixes (Cauchy interlacing). It is a spurious value in the sense that it approximates no eigenvalue, but not an inflated one. Pollution proper needs a non-normal $\K$, e.g. an irreversible chain or the OU of Problem 9(d).

For $f,g\in F$, $\ip{g,\K f}=\ip{g,P_F\K f}=\ip{g,Cf}$ because $g\in F$ makes $\ip{g,(I-P_F)h}=0$ for any $h$. The maximum of $\ip{g,Cf}/(\|f\|\|g\|)$ over $F\times F$ is the top singular value of $C$, attained at its top singular pair (Cauchy–Schwarz in $g$ gives $\|Cf\|/\|f\|$, then maximise over $f$). By Problem 13(e), $F_\tau$ is the matrix of $CC^\dagger$, so its eigenvalues are $\sigma_i(C)^2$ and its eigenvectors the coefficient vectors of the $g$'s; $B_\tau$ likewise gives the $f$'s. In the language of Eq. (4)–(5), the paper's filters are the pair of representable observables whose present and future values are most correlated, which is exactly the CCA problem of Problem 14(a)–(b).

Summing (a) over $i\le d$: $\sum_{i\le d}\sigma_i(C)^2\le\sum_{i\le d}\sigma_i(\K)^2$, the VAMP-2 score bound. Equality holds iff $F$ is spanned by the top $d$ singular functions of $\K$ (then $C$ restricted to them is $\K$ itself, since both the input and the output lie in $F$). This turns feature choice into an optimisation with a known ceiling, and is the precise sense in which "approximating early" is harmless for singular vectors: the compressed problem is a lower bound on a variational objective of the exact operator, and maximising it over feature sets recovers the exact singular functions.

$\|\K^\dagger g\|=\max_{f}\ip{f,\K^\dagger g}/\|f\|=\max_f\ip{\K f,g}/\|f\|$ over all $f$, which is at least the maximum over $f\in F$; so $\max_{g\in F}\|\K^\dagger g\|/\|g\|\ge\sigma_1(C)$. The one-sided problem lets the predictive function leave $F$.

For $F=\operatorname{span}\{x^2\}$ the compressed operator is the scalar $c=\ip{x^2,\K x^2}/\ip{x^2,x^2}$. With $\K x^2=e^{-2\tau}x^2+(1-e^{-2\tau})$: $\ip{x^2,\K x^2}=3e^{-2\tau}+(1-e^{-2\tau})=1+2e^{-2\tau}$, so $c=\tfrac13(1+2e^{-2\tau})$ (as in Problem 13f) and $\sigma_1(C)^2=c^2$. The one-sided value (the process is reversible, $\K^\dagger=\K$) is $\|\K x^2\|^2/3$ with $\|\K x^2\|^2=e^{-4\tau}\cdot3+2e^{-2\tau}(1-e^{-2\tau})\cdot1+(1-e^{-2\tau})^2=1+2e^{-4\tau}$, i.e. $\tfrac13(1+2e^{-4\tau})$. At $\tau=0.5$: $c^2=0.335$ versus $0.424$; the true values $\sigma_i(\K)^2=e^{-2n\tau}$ are $1,\,0.368,\,0.135,\dots$ and both estimates obey their bounds. The gap is Jensen: $x^2=\mathrm{He}_0+\mathrm{He}_2$ has weights $\tfrac13,\tfrac23$ on eigenvalues $1,e^{-2\tau}$, and the one-sided value averages the squares $(1+2e^{-4\tau})/3$ while the two-sided one squares the average $((1+2e^{-2\tau})/3)^2$. The optimal $f$ in the one-sided problem is $\K^\dagger g/\|\K^\dagger g\|\propto\K x^2=e^{-2\tau}x^2+(1-e^{-2\tau})$, which contains a constant and is not in $F$.

The paper solves the two-sided problem: Eq. (24) is $CC^\dagger$, with the compression applied on both sides of the product. That is the right problem for a predictive and a retrospective neuron that each see only the features, and it is what CCA has always computed. But it is a modelling decision about representability, not an approximation to the one-sided Ritz problem, and neither is an approximation to the other. The paper does not flag the choice.

If $\K F\subset F$ and $\K^\dagger F\subset F$ then $F^\perp$ is also invariant under both (for $h\perp F$ and $f\in F$: $\ip{f,\K h}=\ip{\K^\dagger f,h}=0$ since $\K^\dagger f\in F$; likewise for $\K^\dagger$). So $\K$ is block diagonal with respect to $F\oplus F^\perp$, $P_F\K|_F=\K|_F$ exactly, and $C^\dagger C=\K^\dagger\K|_F$; eigenvectors of $C^\dagger C$ are therefore eigenvectors of $\K^\dagger\K$, i.e. singular functions of $\K$, with singular values from the spectrum of $\K$ (a subset of it, the whole top block if $F$ contains the top singular functions).

OU with $\phi(x)=x$: $\K(\vec c\cdot x)=(e^{A^\top\tau}\vec c)\cdot x$ is linear (Problem 13d), and $\K^\dagger$ is the Koopman operator of the reversed process, which by Problem 9 is again OU with drift $A_{\rm rev}=\Sigma A^\top\Sigma^{-1}$, so $\K^\dagger(\vec c\cdot x)=(e^{A_{\rm rev}^\top\tau}\vec c)\cdot x$ is linear too. Both invariances hold, and Problem 14(d) is the resulting equality of Eq. (24) with Eq. (18). (Strictly, $F$ should include the constant so that $\lambda_0=1$ is represented; for zero-mean linear features $1\perp F$ and the block structure still holds.)

For a nonlinear process with lag-vector features $\phi(x)=x$, $\K x$ is a nonlinear function of the lag vector and neither invariance holds. The construction then returns the best linear filter pair: the linear predictive and retrospective observables most correlated across the lag. That is a well-defined and possibly useful object, but a coherent set of the underlying dynamics is a level set of a singular function which may be strongly nonlinear in the lag vector, and there is no theorem that the best linear pair approximates it. The VAMP score certifies only that $\sum\sigma_i(C)^2$ is as large as this feature class allows; without a reference value for $\sum\sigma_i(\K)^2$ it cannot say how far below the ceiling one is. Problem 13(f)(i) is the toy version: the feature set $\{x,x^2\}$ returns a confident $0.579$ where the truth is $0.368$, and nothing in the score flags it.

$r^2=\|\K f\|^2-2\lambda\ip{g,\K f}+\lambda^2\|g\|^2$. For $f,g\in F$ the last two terms are $\vec c_g^{\,\top}\Sigma_\tau\vec c_f$ and $\vec c_g^{\,\top}\Sigma_0\vec c_g$, both from Eq. (23). The first is $\ip{\K f,\K f}$, which is not a bilinear form in $\phi$ evaluated at two times, so no lagged second moment of the features gives it.

Law of total variance: $\E[f(X_{t+\tau})^2\mid X_t]=(\K f)(X_t)^2+\operatorname{Var}(f(X_{t+\tau})\mid X_t)$. Taking expectations, $\|f\|^2=\|\K f\|^2+\E[\operatorname{Var}(f(X_{t+\tau})\mid X_t)]$, using stationarity on the left. So the single-trajectory quantity one can estimate, $\E[f(X_{t+\tau})^2]$, overestimates $\|\K f\|^2$ by the mean conditional variance, and this is the same gap as between the one-sided and two-sided values in (c) (the contraction bound of Problem 7a is this inequality). For the sibling identity, condition on the present: $\E[f(X_{t+\tau})f(X'_{t+\tau})\mid X_t]=(\K f)(X_t)^2$ by conditional independence, and average. So $\|\K f\|^2=\ip{f,\K^\dagger\K f}$ is a two-trajectory quantity, the sibling kernel of Problem 10, and estimating it needs pairs of futures launched from the same present. A single time series offers each present exactly one future. In a deterministic system the conditional variance vanishes, $\|\K f\|^2=\|f\|^2$, and the residual is computable from one trajectory; that is why residual DMD works there and why its stochastic extension needs extra sampling.

Both identities by direct integration. With $\K f(x)=\int p_\tau(y\mid x)f(y)\,dy$, $$\|\K f\|^2=\int\rho_0(x)\Big[\int p_\tau(y\mid x)f(y)\,dy\Big]^2dx=\iiint\rho_0(x)\,p_\tau(y\mid x)\,p_\tau(y'\mid x)\,f(y)f(y')\,dx\,dy\,dy' ,$$ writing the square as a double integral over two dummy futures. The weight $\rho_0(x)p_\tau(y\mid x)p_\tau(y'\mid x)$ is the joint density of a present $X_t\sim\rho_0$ and two conditionally independent futures launched from it (product of the two conditionals given the shared start), so the triple integral is $\E[f(X_{t+\tau})f(X'_{t+\tau})]$: the sibling identity. For the variance identity, $\int p_\tau(y\mid x)f(y)^2\,dy-\big[\int p_\tau(y\mid x)f(y)\,dy\big]^2=\operatorname{Var}(f(X_{t+\tau})\mid X_t=x)$ is the definition of a variance (second moment minus squared mean) under the kernel $p_\tau(\cdot\mid x)$. Multiply by $\rho_0(x)$ and integrate over $x$: the first term becomes $\iint\rho_0(x)p_\tau(y\mid x)f(y)^2\,dx\,dy=\int\rho_0(y)f(y)^2\,dy=\|f\|^2$, doing the $x$-integral first and using the invariance of $\rho_0$, $\int\rho_0(x)p_\tau(y\mid x)\,dx=\rho_0(y)$ (Problem 7); the second term is $\|\K f\|^2$ by the display above. Hence $\|f\|^2-\|\K f\|^2=\int\rho_0(x)\operatorname{Var}(f(X_{t+\tau})\mid X_t=x)\,dx$, the law of total variance in integral form. The single-trajectory quantity $\E[f(X_{t+\tau})^2]$ is the first term; what is missing is the triple integral, which needs two futures per present.

For the pair of (c), $f=g=x^2/\sqrt3$, $\lambda=c$: $r^2=\|\K f\|^2-c^2=\tfrac13(1+2e^{-4\tau})-\tfrac19(1+2e^{-2\tau})^2=0.424-0.335=0.089$ at $\tau=0.5$, against $r^2=0$ for the true pair $(\mathrm{He}_2/\sqrt2,\,e^{-2\tau})$. The residual would expose the compressed pair, but only with $\|\K f\|^2$ in hand. The script estimates it two ways: from sibling futures (200 000 stationary presents, two Euler–Maruyama futures each) it finds $0.428$ against the analytic $0.424$, giving $r^2\approx0.10$; from the single-trajectory proxy $\E[f(X_{t+\tau})^2]/\|f\|^2=1.000$ it would report $r^2\approx0.67$, wrong by the mean conditional variance $0.576$. The neuron has the second kind of data, not the first.

If $t$ and $t+s$ fall on the same plateau, $y(t)y(t+s)=L^2$ with $\E L^2=\sigma^2$; if on different plateaus, the levels are independent with zero mean and the product averages to zero. Hence $C(s)=\sigma^2\Pi(s)$. In a stationary renewal process the plateau containing a typical time is length-biased: the residual length $R$ from $t$ to the next edge has density $P(T>u)/\E T$, so $\Pi(s)=P(R>s)=\frac1{\E T}\int_s^\infty P(T>u)\,du$. For $T\sim\mathrm{Exp}(\lambda)$, $P(T>u)=e^{-\lambda u}$ and $\Pi(s)=e^{-\lambda s}$.

A one-dimensional OU process $dx=-\lambda x\,dt+\sqrt{2\lambda\sigma^2}\,dW$ has the same autocovariance. Linear CCA between lag vectors depends on the process only through $\Sigma_0$ and $\Sigma_\tau$, which are built from $C$ (Problem 13b), so the dead-leaves signal and this OU give identical canonical correlations and directions. For the OU the lag vector is a noisy image of a one-dimensional state, so exactly one canonical correlation is non-zero, equal to $e^{-\lambda\tau}$ (in the adjacent-window construction the two windows meet at $y(t)$ and $y(t+1)$, and the correlation is $e^{-\lambda}$ per sample). The script confirms it: correlations $0.95,\,0.01,\,0.01,\dots$ for exponential plateaus with mean length $20$. A transient neuron would have nothing to predict.

If $\Pi$ is not a single exponential, $C$ is not the autocovariance of any one-dimensional linear process, and the lag vector carries more than one linearly predictable direction (a sum of $m$ exponentials corresponds to an $m$-dimensional OU, Problem 13d generalised). For Pareto lengths the script finds correlations $0.92,\,0.21,\,0.14,\dots$: a clear second pair.

Condition on the present level $L$ and on the age $a$ (time since the last edge). Over the horizon either no edge occurs, and $y(t+\tau)=L$, or one does, and the new level has mean zero. So $\E[y(t+\tau)\mid L,a]=L\,P(\text{no edge in }[t,t+\tau]\mid a)=L\,\Pi_a(\tau)$ with $\Pi_a(\tau)=P(T>a+\tau)/P(T>a)$. For heavy-tailed $T$ this is increasing in $a$ (decreasing hazard: a plateau that has already lasted a long time is likely to last longer), while for exponential $T$ it is constant. Numerically, the probability that the level changes within $\tau$ is $0.84$ for plateaus younger than three samples and $0.16$ for plateaus older than thirty; for exponential lengths it is $0.45$ at every age.

So a recent edge means a young plateau, and a young plateau means the current level is a poor predictor of the future: the level will most likely be redrawn before $t+\tau$, i.e. revert toward the mean. That is a second-moment effect visible to a linear filter because the signed recent jump $y(t)-y(t-k)$ is correlated with the current level (after an up-step the level tends to be high) and, through the age it reveals, with the persistence of that level. The part of the jump uncorrelated with the level, which is what the second canonical variate is, predicts the correction to the level-based forecast.

The future filters are the lag-reversed past filters (part d). Pair 1 is the level: its whitened filter is a decaying low-pass shape like the paper's dashed curve, even though in raw lag coordinates it is essentially $y(t)$ alone; whitening by $\Sigma_0^{1/2}$ smears a delta into the first column of $\Sigma_0^{1/2}$, which for a slowly decorrelating signal is a low-pass filter. Pair 2 is uncorrelated with the level by construction and reads the signed recent jump with negative sign: "just stepped down" gives a positive past variate. Its future partner, oriented so the pair is positively correlated, reads "the level at $t+\tau$ is high, and higher than now". So the coherent pair is

just stepped down $\;\longrightarrow\;$ will be higher by $t+\tau$, $\qquad$ just stepped up $\;\longrightarrow\;$ will be lower by $t+\tau$.

An edge does not predict a flat future; every moment is flat except at edges. It predicts that the current level will not last, and since levels are redrawn independently, that the signal will move back toward the mean. The ON transient unit, $H(+\mathbf w^\top\hat{\mathbf x})$ with $\mathbf w$ the derivative filter, is the membership index of "recent up-step, expect a decrease"; the OFF unit is its complement. The sustained pair reports the level; the transient pair reports, in effect, how much the level should be trusted. The whitened shapes match Fig. 2b, and the sustained/transient split of bipolar cells is read as the two most predictable directions of this input class. High for the mechanism in this model; whether real dead-leaves scans have the required heavy tail at the relevant scale, and whether the second correlation survives the retina's noise, is an empirical question the paper does not test directly.

Reversibility of $y$ means $(y(s))_s$ and $(y(-s))_s$ have the same law. The lag vector at time $t$ is $\hat{\mathbf x}(t)=(y(t),\dots,y(t-d+1))$; under time reversal the sample that was most recent becomes the oldest, so the reversed process's lag vector is $J\hat{\mathbf x}$, not $\hat{\mathbf x}$. Hence $[\Sigma_\tau]_{ij}=\E[y(t-i)\,y(t+\tau-j)]=C(\tau-j+i)$, whereas $[\Sigma_\tau]_{ji}=C(\tau-i+j)$: not equal in general. The script reports $\|\Sigma\text{-whitened }C-C^\top\|_\infty\approx0.4$. What reversibility does give is $[\Sigma_\tau]_{ij}=C(\tau+i-j)=[\Sigma_\tau]_{d-1-j,\,d-1-i}$, i.e. $\Sigma_\tau^\top=J\Sigma_\tau J$; the script finds this identity holds to $10^{-3}$. Since $\Sigma_0$ is Toeplitz, $J\Sigma_0J=\Sigma_0$ as well, so the whitened cross-covariance satisfies $\tilde C^\top=J\tilde CJ$. If $\tilde C\tilde{\mathbf b}=\rho\tilde{\mathbf a}$ and $\tilde C^\top\tilde{\mathbf a}=\rho\tilde{\mathbf b}$, then applying $J$ to the second and using $J\tilde C^\top J=\tilde C$ gives $\tilde C(J\tilde{\mathbf a})=\rho(J\tilde{\mathbf b})$: the pair $(J\tilde{\mathbf b},J\tilde{\mathbf a})$ is also a singular pair with the same $\rho$, so for a non-degenerate $\rho$, $\tilde{\mathbf b}=\pm J\tilde{\mathbf a}$ and likewise $\mathbf b=\pm J\mathbf a$. The table in (c) shows exactly this: each future filter is the past filter read backwards.

The paper's phrasing is this fact in plotting coordinates. Its Fig. 2b plots the past filters against time lag on $[-1,0]$; the future filter of a pair, plotted against time in the future window, is the mirror image ("invert the time axis"). The sign is a convention: with the orientation used here, where each pair is positively correlated, the mirrored derivative filter reads "will go down" when the past one reads "just went up" with no extra sign; if instead one insists that the mirrored filter keep the past filter's sign convention ("late minus early" in both windows), the pair needs a minus sign, which is the paper's "inverting the sign for the second". Either way the content is the same: the transient pair encodes mean reversion, and the sustained pair does not need a sign because $y(t)$ predicts $y(t+\tau)$ with a positive coefficient. The general lesson for Problem 14(b)(ii): "reversible implies $\Sigma_\tau$ symmetric" is a statement about features of a Markov state, and a lag vector is a Markov state only up to the reordering that time reversal induces.