Realize reducible nonstationary kernels as solution to SDEs and its extensions I am interested in a regression application where my kernel is of the form
\begin{equation}
    k(t,t^{\prime}) = k_s\left(\phi(t),\phi(t^{\prime})\right)= k_s\left(\phi(t)-\phi(t^{\prime})\right),
\end{equation}
where $k_s$ is a stationary kernel. I came across this paper (section 5) which classifies such kernels as reducible kernels. 
I am using such a kernel for time series extrapolation via Gaussian process regression. Assuming I know $\phi(\cdot)$, is there a way to leverage that information to do the regression in the space of the stationary kernel instead of the nonstationary one? 
It is known that for a time series application, Gaussian process regression with a stationary kernel can be done via Kalman smoothing (see here). This is possible via synthesizing a stationary process from white noise.
\begin{eqnarray}
  x(t) &=& w(t) * h(t) \\
  k_x(\tau) &=& h(\tau)*h^\star(-\tau) \\
  S_x(\omega) &=& H(\omega) H^\star(\omega) = \left|H(\omega) \right|^2 
\end{eqnarray}
Therefore, by finding $h(t)$ or equivalently $H(\omega)$, one can express a GP in terms of white noise. Matern kernels accept such decompositions naturally and can be expressed in terms of differential equations of white noise. For other stationary kernels like RBF, aforementioned paper proposes Taylor approximation. 
To give a simple example; Ornstein-Uhlenbeck process with the kernel $k_{OU}(t,t^{\prime}) = \exp\left(-\lambda \left|t-t^\prime \right| \right)=\exp(-\lambda \left| \tau \right|)$ and the spectral density $S_{OU}(\omega) = \frac{2\lambda}{\omega^2 + \lambda^2}$ can expressed as 
\begin{equation}
   \frac{dx}{dt} = - \lambda x(t) + w(t)
\end{equation}
Hence; whole GP regression can be reduced into a standard filtering/smoothing problem. This results in linear computational burden instead of a cubic one, and also naturally fits to the temporal nature of data.
I would like to know whether this is extensible to the reducible nonstationary kernels. For instance if my kernel is $k(t,t^\prime) = k_{OU}(\phi(t),\phi(t^\prime))$ where $\phi(\cdot)$ is some known function, can I somehow utilize the above differential equation?
tldr, are there any benefits if you have a reducible kernel?
 A: $\blacksquare$1.Whether this(the method of SDE solving kernel problem) is extensible to the reducible nonstationary kernels?
Let me restate your question. The method your pointed out in [Hartikainen&Särkkä] is to regard the covariance(OR its spectral form) as a solution to stochastic differential equation. Your question is that if $k(t,t')$ is not a stationary kernel, whether it can still be expressed as a solution to SDE. 
The answer is NO in general. Because kernels can be regarded as inner products on the state space families indexed by the time variable $t\in T$, if you have a nonstationary kernel i.e. $k(t,t')\not\propto |t-t'|$ then you are saying that the families of inner products vary with time. So the best possibility for general kernels is that you get a family of SDEs, whose coefficients varying with time variable $t$, that at each fixed time $t_0$ you get a SDE with $Spectral_{k}(t,t')$ as solution at this specific time point. 
But I guess for special family of kernels $\color{red}{see\, my\, update}$, say the kernel is reducible in sense of [Genton], this result can be extended with coefficient of SDEs varying with the connecting function $\phi$ you pointed out in your OP. 
$\blacksquare$2.Are there any benefits if you have a reducible kernel?
First this is not a "reducible kernel" as I usually used in a representation morphism, it is "stationary reducible kernel" defined as the author in [Genton] said:

We say that a nonstationary kernel $K(x, z)$ is stationary
  reducible if there exist a bijective deformation $\phi$ such that:
  $K(x, z) = K^{*}_{S}(\phi(x) − \phi(z))$

The importance of studying such a kind of kernel traces back to [Sampson&Guttorp], which use a thin-plate spline to recover the missing dimension used to recover stationary-ness. The main idea is that lower dimension nonstationary kernel must be a projection of stationary kernel in higher dimension. There are two basic drawbacks of this idea in practice, one drawback is scale-sensitive; the other is the difficulty of separating the noise and mean.  There are earlier motivations of this idea from mathematics but [Sampson&Guttorp] is the first literature that brings this idea to statisticians.
Reference
[Hartikainen&Särkkä]Hartikainen, Jouni, and Simo Särkkä. "Kalman filtering and smoothing solutions to temporal Gaussian process regression models." 2010 IEEE International Workshop on Machine Learning for Signal Processing. IEEE, 2010.
[Genton]Genton, Marc G. "Classes of kernels for machine learning: a statistics perspective." Journal of machine learning research 2.Dec (2001): 299-312.
[Karlin&Taylor]Karlin, Samuel, and Howard E. Taylor. A second course in stochastic processes. Elsevier, 1981.
[Sampson&Guttorp]Sampson, Paul D., and Peter Guttorp. "Nonparametric estimation of nonstationary spatial covariance structure." Journal of the American Statistical Association 87.417 (1992): 108-119.

Update
My knowledge on nonstationary kernels is mainly from spatial modeling and group representation theory instead of "machine learning", so here is my best knowledge but possibly a bit different from what you want. One highly cited "geological" paper that explained the mathematical difficulty of deriving such a family of SDE indexed by time variable is [Neuman].
Nonstationary covariance function has a special problem of differentiability when we try to obtain a solution family the corresponding SDE family in classical way(It is easily understood because you cannot expect the spectral function to be regular if the underlying manifold defined by the kernel process is not smooth). To visualize this point you can imagine that a stochastic process without smooth path.
One well-known example that the covariance kernel is not only nonstationary but also nonsmooth is given by setting up a Dirichlet process as prior in $M(\mathcal{X})$, the space of a measures on a Polish sample space $\cal{X}$, whose path is almost nowhere smooth. If the (posterior) covariance is determined by a Dirichlet process, say in a Bayesian nonparametric setting like [Gelfand et.al], then the resulting covariance kernel is not smooth let alone stationary. This example will correspond to a local velocity field which is also not smooth. Hence it is expected that such a covariance kernel is nonstationary. If the resulting $K(t,t')$ is not smooth, then the same differentiability problem arise when we try to figure out the connecting function $\phi$ or $\phi_t$(varying with $t$). It is a practice (Disclaimer: I do not know whether it is common or not in statistical community!) that sometimes we want to assume sort of "local differentiability" for $K$ or $\phi$ in order to solve the SDE for a fixed $t_0$. 
Moreover if you assume some smoothness then you are actually assuming some sort of (local) invertibility, you cannot assume non-invertible smooth kernels by implicit function theorem. So your comment to @Dougal seems a bit confusing to me.
One method of overcoming this is to define weak derivatives and try to solve it as PDE by regarding the time variable as a new variate; another method is aforementioned by [Sampson&Guttorp] and later by Perrin as we discussed above.
As for simple example, I'm afraid that there is no a simple/non-pathological example using Gaussian process since its path is smoothly differentiable everywhere.
One example for "special family" I mentioned in my answer to Question 1 above is given in [Solin&Särkkä].
A question for you: what is your motivation of asking such a problem from viewpoint of "machine learning"? Thanks, this is a good discussion!
Update Reference
[Paciorek&Schervish] Paciorek, Christopher J., and Mark J. Schervish. "Spatial modelling using a new class of nonstationary covariance functions." Environmetrics 17.5 (2006): 483-506.
[Neuman]Neuman, Shlomo P. "Eulerian‐Lagrangian theory of transport in space‐time nonstationary velocity fields: Exact nonlocal formalism by conditional moments and weak approximation." Water Resources Research 29.3 (1993): 633-645.
[Solin&Särkkä]Solin, Arno, and Simo Särkkä. "Explicit Link Between Periodic Covariance Functions and State Space Models." AISTATS. 2014.
A: Denote the value of your process $u(t)$. First, note that if you consider $t_1 \neq t_2$ such that $\phi(t_1) = \phi(t_2)$ then you must have $u(t_1) = u(t_2)$. How do we know this? $u$ is Gaussian, so $u(t_1)$ and $u(t_2)$ are jointly Gaussian. But 
$\text{Cov}[ u(t_1),u(t_2) ] = k_s(0) = \text{Var}[u(t_1)] = \text{Var}[u(t_2)]$. Now you have two ways to look at the situation. First, the (Gaussian!) RVs $u(t_1),u(t_2)$ have correlation $1$, so they are identical. Alternatively, you can look at their covariance matrix. All four of its entries are equal, so it is degenerate. This means, again, that $u(t_1) = u(t_2)$. So, in this respect, you really want to do all inference in $x = \phi(t)$ space. Given a new $t^{*}$ for which you want to predict the distribution of $u(t^{*})$, find $x^{*} = \phi(t^{*})$ and your prediction for $t^{*}$ should be that of $x^{*}$. This does not require any invertibility assumptions on $\phi$.
SDE approach with Ito's Lemma is not as straightforward as I thought...
