Gaussian process - dimensionality reduction

Question

Specific question on Gaussian Processes and dimensionality reduction. I saw a a method for dimensionality reduction for the squared exponential covariance function (not ARD) whereby one uses a GxD projection matrix P (G < D, D = dimension of the inputs) such that the squared distances are calculated in a lower dimensional space. E.g. the kernel function K(x, x') = c * exp(-1/2*(P(x-x'))T*(P(x-x'))) with T denoting transpose. I'm trying to use this in my implementation of gpr, but I'm not sure how to get the gradient with respect to P for optimization. E.g. I need the derivative of the marginal likelihood w.r.t. P to maximize (well minimizing the negative log marginal likelihood in actuality). Has anyone seen this method before or have any suggestions? Thanks.