Let’s derive the centered vs non-centered versions of these. @Bob_Carpenter is right about the differences between models 1 and 2.
Let g_{\text{free}} be a typical, unconstrained GP
g_{\text{free}} \sim \mathcal{N}(0, \tau^2 K).
One key element to my Stancon 2026 talk is that we derive the distribution of the sum of g_{\text{free}} where we will call it S_{\text{free}} = \sum g_{\text{free}}. The distribution is
S_{\text{free}} \sim \mathcal{N}(0, \tau^2 \bf{1}^{\top} K \bf{1})
where \bf{1} is a row vector of ones. The covariance between the vector g_{\text{free}} and S_{\text{free}} is
\begin{aligned}
\mathrm{Cov} &= \mathrm{Cov}(g, \mathbf{1}^{\top} g) \\
&=\mathrm{Cov}(g, g) \mathbf{1} \\
&= \tau^2 K \mathbf{1}.
\end{aligned}
We can do two cool things here. The first is that we can regress g on the sum S (dropping the free subscripts for now) to get
a = \frac{\mathrm{Cov}(g, S)}{\mathrm{Var}(S)} = \frac{\tau^2 K \mathbf{1}}{\tau^2 \mathbf{1}^{\top} K \mathbf{1}} = \frac{K \mathbf{1}}{\mathbf{1}^{\top} K \mathbf{1}}.
The residual sum of squares is 0 thus let the residual squares be
r = g - a S
then \sum r = 0. We can derive the covariance matrix of r as
\begin{aligned}
\mathrm{Cov}(r) &= \mathrm{Cov}(g - a S) \\
&= \mathrm{Var}(g) + a \mathrm{Var}(S) a^{\top} - \mathrm{Cov}(g, S) a^{\top} - a \mathrm{Cov}(g, S) \\
&= \tau^2 \underbrace{ \left(K - \frac{K \mathbf{1}\mathbf{1}^{\top} K}{\mathbf{1}^{\top} K \mathbf{1}} \right )}_{K_0}.
\end{aligned}
Now, with those preliminaries out of the way, we get to our sum-to-zero vector by applying the s2z transform: z = H y_{\text{free}} where y is size J - 1 and z is size J. We act on the J - 1 subspace of J where everything sums to 0. The covariance happens to be the same as derived above with the projection to J that lies on the sum-to-zero plane:
\begin{aligned}
C &= H^{\top} \mathrm{Cov}(r) H \\
&= \tau^2 H^{\top} K_0 H
\end{aligned}
The centered parameterization has one issue where if we use the s2z vector we’ve placed an additional Gaussian element that we need to remove which we can below:
data {
int<lower=2> J;
array[J] real x;
}
transformed data {
vector[J] q = rep_vector(inv_sqrt(1.0 * J), J);
vector[J] mu = rep_vector(0.0, J);
}
parameters {
sum_to_zero_vector[J] z;
real<lower=0> tau;
real<lower=0> rho;
}
transformed parameters {
matrix[J, J] K = gp_exp_quad_cov(x, 1.0, rho);
matrix[J, J] L_K = cholesky_decompose(K);
}
model {
tau ~ normal(0, 1);
rho ~ lognormal(0, 1);
z ~ multi_normal_cholesky(mu, tau * L_K);
/*
* Normalize the Gaussian density after conditioning on q' z = 0.
* This matters when tau or rho is estimated.
*/
target -= normal_lpdf(
0.0 | 0.0,
tau * sqrt(dot_product(q, K * q))
);
}
generated quantities {
real sum_z = sum(z);
}
The non-centered
functions {
/**
* Orthonormal basis for the sum-to-zero subspace.
*
* H' * H = I_(J-1)
* H' * 1 = 0
*/
matrix s2z_basis(int J) {
matrix[J, J - 1] H = rep_matrix(0.0, J, J - 1);
for (k in 1:(J - 1)) {
real a = inv_sqrt(k * (k + 1.0));
for (j in 1:k) {
H[j, k] = a;
}
H[k + 1, k] = -k * a;
}
return H;
}
}
data {
int<lower=2> J;
array[J] real x;
}
transformed data {
matrix[J, J - 1] H = s2z_basis(J);
vector[J] q = rep_vector(inv_sqrt(1.0 * J), J);
}
parameters {
vector[J - 1] eta;
real<lower=0> tau;
real<lower=0> rho;
}
transformed parameters {
matrix[J, J] K = gp_exp_quad_cov(x, 1.0, rho);
vector[J] Kq = K * q;
matrix[J, J] K_zero = K - (Kq * Kq') / dot_product(q, Kq);
matrix[J - 1, J - 1] C = cholesky_decompose(quad_form(H, K_zero));
matrix[J - 1, J - 1] L_C;
vector[J - 1] y = tau * (L_C * eta);
vector[J] z = sum_to_zero_constrain(y); //should be same as H * y
}
model {
eta ~ std_normal();
tau ~ std_normal();
rho ~ lognormal(0, 1);
}
generated quantities {
real sum_z = sum(z);
}
I bet there’s a faster way to do quad_form(H, K_zero) since we know H.