r - 在 Stan/rstan 中优化高斯过程

Question

我最近遇到了高斯过程模型，并且碰巧认为它们可能是我实验室一直在研究的问题的解决方案。我有一个关于交叉验证的开放和相关问题，但我想将我的建模/数学问题与我的编程问题分开。因此，这是第二个相关帖子。如果更多地了解我的问题的背景会有所帮助，尽管这里是我打开的链接 CV question .

这是我的 stan 代码，对应于我的 CV 帖子中提供的更新的协方差函数:

functions{
    //covariance function for main portion of the model
    matrix main_GP(
        int Nx,
        vector x,
        int Ny,
        vector y, 
        real alpha1,
        real alpha2,
        real alpha3,
        real rho1,
        real rho2,
        real rho3,
        real rho4,
        real rho5,
        real HR_f,
        real R_f){
                    matrix[Nx, Ny] K1;
                    matrix[Nx, Ny] K2;
                    matrix[Nx, Ny] K3;
                    matrix[Nx, Ny] Sigma;

                    //specifying random Gaussian process that governs covariance matrix
                    for(i in 1:Nx){
                        for (j in 1:Ny){
                            K1[i,j] = alpha1*exp(-square(x[i]-y[j])/2/square(rho1));
                        }
                    }

                    //specifying random Gaussian process incorporates heart rate
                    for(i in 1:Nx){
                        for(j in 1:Ny){
                            K2[i, j] = alpha2*exp(-2*square(sin(pi()*fabs(x[i]-y[j])*HR_f))/square(rho2))*
                            exp(-square(x[i]-y[j])/2/square(rho3));
                        }
                    }

                    //specifying random Gaussian process incorporates heart rate as a function of respiration
                    for(i in 1:Nx){
                        for(j in 1:Ny){
                            K3[i, j] = alpha3*exp(-2*square(sin(pi()*fabs(x[i]-y[j])*HR_f))/square(rho4))*
                            exp(-2*square(sin(pi()*fabs(x[i]-y[j])*R_f))/square(rho5));
                        }
                    }

                    Sigma = K1+K2+K3;
                    return Sigma;
                }
    //function for posterior calculations
    vector post_pred_rng(
        real a1,
        real a2,
        real a3,
        real r1, 
        real r2,
        real r3,
        real r4,
        real r5,
        real HR,
        real R,
        real sn,
        int No,
        vector xo,
        int Np, 
        vector xp,
        vector yobs){
                matrix[No,No] Ko;
                matrix[Np,Np] Kp;
                matrix[No,Np] Kop;
                matrix[Np,No] Ko_inv_t;
                vector[Np] mu_p;
                matrix[Np,Np] Tau;
                matrix[Np,Np] L2;
                vector[Np] yp;

    //--------------------------------------------------------------------
    //Kernel Multiple GPs for observed data
    Ko = main_GP(No, xo, No, xo, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
    Ko = Ko + diag_matrix(rep_vector(1, No))*sn;

    //--------------------------------------------------------------------
    //kernel for predicted data
    Kp = main_GP(Np, xp, Np, xp, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
    Kp = Kp + diag_matrix(rep_vector(1, Np))*sn;

    //--------------------------------------------------------------------
    //kernel for observed and predicted cross 
    Kop = main_GP(No, xo, Np, xp, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);

    //--------------------------------------------------------------------
    //Algorithm 2.1 of Rassmussen and Williams... 
    Ko_inv_t = Kop'/Ko;
    mu_p = Ko_inv_t*yobs;
    Tau=Kp-Ko_inv_t*Kop;
    L2 = cholesky_decompose(Tau);
    yp = mu_p + L2*rep_vector(normal_rng(0,1), Np);
    return yp;
    }
}

data { 
    int<lower=1> N1;
    int<lower=1> N2;
    vector[N1] X; 
    vector[N1] Y;
    vector[N2] Xp;
    real<lower=0> mu_HR;
    real<lower=0> mu_R;
}

transformed data { 
    vector[N1] mu;
    for(n in 1:N1) mu[n] = 0;
}

parameters {
    real loga1;
    real loga2;
    real loga3;
    real logr1;
    real logr2;
    real logr3;
    real logr4;
    real logr5;
    real<lower=.5, upper=3.5> HR;
    real<lower=.05, upper=.75> R;
    real logsigma_sq;
}

transformed parameters {
    real a1 = exp(loga1);
    real a2 = exp(loga2);
    real a3 = exp(loga3);
    real r1 = exp(logr1);
    real r2 = exp(logr2);
    real r3 = exp(logr3);
    real r4 = exp(logr4);
    real r5 = exp(logr5);
    real sigma_sq = exp(logsigma_sq);
}

model{ 
    matrix[N1,N1] Sigma;
    matrix[N1,N1] L_S;

    //using GP function from above 
    Sigma = main_GP(N1, X, N1, X, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
    Sigma = Sigma + diag_matrix(rep_vector(1, N1))*sigma_sq;

    L_S = cholesky_decompose(Sigma);
    Y ~ multi_normal_cholesky(mu, L_S);

    //priors for parameters
    loga1 ~ student_t(3,0,1);
    loga2 ~ student_t(3,0,1);
    loga3 ~ student_t(3,0,1);
    logr1 ~ student_t(3,0,1);
    logr2 ~ student_t(3,0,1);
    logr3 ~ student_t(3,0,1);
    logr4 ~ student_t(3,0,1);
    logr5 ~ student_t(3,0,1);
    logsigma_sq ~ student_t(3,0,1);
    HR ~ normal(mu_HR,.25);
    R ~ normal(mu_R, .03);
}

generated quantities {
    vector[N2] Ypred;
    Ypred = post_pred_rng(a1, a2, a3, r1, r2, r3, r4, r5, HR, R, sigma_sq, N1, X, N2, Xp, Y);
}

我已经修改了内核中包含的参数的先验知识，有些参数化要快一些(在某些情况下最多快两倍，但即使在这些情况下仍然可以产生相对较慢的链)。

我正在尝试使用受污染部分前后 15 秒的数据(以 3.33 Hz 采样，总共 100 个数据点)来预测 3.5 秒的数据值(以 10 Hz 采样 - 因此有 35 个数据点)。

R中指定的模型如下:

 fit.pred2 <- stan(file = 'Fast_GP6_all.stan',
                 data = dat, 
                 warmup = 1000,
                 iter = 1500,
                 refresh=5,
                 chains = 3,
                 pars= pars.to.monitor
                 )

老实说，我不知道我是否需要那么多热身迭代。我想估计缓慢的部分原因是相当无信息先验的结果(心率和呼吸除外 HR 和 R 因为它们在健康成年人的休息时具有相当广为人知的范围)。

任何建议都非常受欢迎，以加快我的程序的运行时间。

谢谢。

Answer 1

你可以获取 Stan Math Library 的开发分支，它有一个最近更新的版本 multi_normal_cholesky在内部使用解析梯度而不是 autodiff。为此，您可以在 R 中执行
source("https://raw.githubusercontent.com/stan-dev/rstan/develop/install_StanHeaders.R")
但你需要有 CXXFLAGS+=-std=c++11在您的 ~/.R/Makevars 文件中，可能需要重新安装 rstan 之后打包。

两者 multi_normal_cholesky和您的 main_GP是 O(N^3)，所以你的 Stan 程序永远不会特别快，但这两个函数的增量优化将产生最大的不同。

除此之外，还有一些小事，比如
Sigma = Sigma + diag_matrix(rep_vector(1, N1))*sigma_sq;
应该改为
for (n in 1:N1) Sigma[n,n] += sigma_sq;
原因是乘以sigma_sq由对角矩阵 puts N1平方节点到 autodiff 树上，将其添加到 Sigma ，它做了大量的内存分配和释放。沿对角线的显式循环仅放置 N1节点到 autodiff 树上，或者如果我们聪明地使用 +=，它可能只是更新现有树。运算符(operator)。您的 post_pred_rng 中的相同交易函数，尽管这不太重要，因为生成的数量块是用 double 计算的，而不是用于反向模式 autodiff 的自定义 Stan 类型。

我认为/希望
vector[N2] Ypred = post_pred_rng(...);
略快于
vector[N2] Ypred; // preallocates Ypred with NaNs
Ypred = post_pred_rng(...);
通过避免预分配步骤；无论如何，读起来更好。

最后，虽然它不会影响速度，但您没有义务以对数形式指定参数，然后在转换后的参数块中对它们进行反对数。你可以用 <lower=0> 声明东西这将导致同样的事情，尽管那样你会在积极约束的事物而不是不受约束的事物上指定你的先验。在大多数情况下，这更直观。那些具有 3 个自由度的学生 t 先验是非常重尾的，这可能导致 Stan 至少在热身期间采取很多跳跃式步骤(默认情况下最多为 10 步)。由于每次迭代需要2^s - 1，所以跳跃步数(s)是运行时间的主要贡献者。函数/梯度评估。