Deriving _rotated_ values for GP samples obtained from gp.conditional()

I’m doing extrapolation from the posterior of a GP, then re-sampling the model as I gather actual data for the extrapolation points. It seems that the starts argument requires samples for both the user-facing GP variables as well as their hidden counterparts with the suffix _rotated, which exist for the values previously inferred via sampling but not for values inferred via gp.conditional(). Has anyone already worked out how to transform the unrotated samples to obtain their rotated counterparts?

I think I worked it out, but would appreciate someone more math-adept double-checking my work:

from pytensor.tensor.linalg import cholesky
from pytensor.tensor.slinalg import cho_solve

with predictive_model:
    predicted_gp_samples = gp.conditional('predicted_gp_samples', Xnew = time_to_predict)
    predicted_gp_samples_rotated_ = pm.Deterministic(
        'predicted_gp_samples_rotated_'
        , cho_solve(
            (cholesky(gp.cov_func(time_to_predict)),True)
            , predicted_gp_samples - gp.mean_func(time_to_predict)
        )
    )
    predict_trace = pm.sample_posterior_predictive(
        trace = trace
        , var_names = [
            'predicted_gp_samples'
            , 'predicted_gp_samples_rotated_'
        ]
        , model = predictive_model
        , idata_kwargs = {'include_transformed': True}
    )