Оптимизация гауссовского процесса в Станции
Недавно я столкнулся с гауссовыми моделями процессов и случайно подумал, что они могут стать решением проблемы, над которой я работал в своей лаборатории. У меня есть открытый и связанный вопрос о перекрестной проверке, но я хотел отделить мои вопросы по моделированию / математике от моих вопросов по программированию. Отсюда и второй, связанный пост. Если бы знание о происхождении моей проблемы помогло бы, хотя вот ссылка на мой открытый вопрос 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);
}
Я возился с априорами для параметров, включенных в мои ядра, некоторые параметризации работают немного быстрее (в некоторых случаях до двух раз быстрее, но все же могут создавать относительно медленные цепочки даже в этих случаях).
Я пытаюсь предсказать значения для данных на 3,5 с (выборка с частотой 10 Гц, то есть с 35 точками данных), используя данные за 15 секунд до и после загрязненного участка (с выборкой с частотой 3,33 Гц, то есть 100 точек данных).
Модель, как указано в 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 поскольку у них есть довольно хорошо известные диапазоны в состоянии покоя у здорового взрослого).
Любые предложения более чем приветствуются, чтобы ускорить время выполнения моей программы.
Благодарю.
1 ответ
Вы можете взять ветвь разработки библиотеки Stan Math, которая имеет недавно обновленную версию multi_normal_cholesky который использует аналитические градиенты внутри вместо авторазличия. Для этого вы можете выполнить в 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 диагональной матрицей ставит N1 узлы в квадрате на дерево автоопределения, как и добавление его в Sigma, который делает много выделения памяти и освобождения. Явный цикл по диагонали только ставит N1 узлы на дереве автоопределения, или, возможно, он просто обновляет существующее дерево, если мы умны с += оператор. Та же сделка внутри вашего post_pred_rng функция, хотя это менее критично, потому что сгенерированный блок количеств оценивается с помощью double, а не пользовательского типа Stan для автодифферсии в обратном режиме
Я думаю / надеюсь, что
vector[N2] Ypred = post_pred_rng(...);
немного быстрее чем
vector[N2] Ypred; // preallocates Ypred with NaNs
Ypred = post_pred_rng(...);
избегая шага предварительного выделения; в любом случае, читать приятнее.
Наконец, хотя это не влияет на скорость, вы не обязаны указывать свои параметры в форме журнала, а затем антилогировать их в блоке преобразованных параметров. Вы можете просто объявить вещи с <lower=0> и это приведет к тому же самому, хотя тогда вы будете указывать свои приоритеты на положительно ограниченных вещах, а не на неограниченных вещах. В большинстве случаев это более интуитивно понятно. Те ученики-ученики с 3 степенями свободы очень тяжелые хвосты, что может привести к тому, что Стэн совершит много прыжковых шагов (по умолчанию до 10), по крайней мере, во время прогрева. Количество шагов (скачков) чехарды является основным фактором времени выполнения, поскольку каждая итерация требует 2^s - 1 оценка функции / градиента.