Is it possible to perform analytically the following integral? $$ \int d^n a_\alpha \exp \left[ -\frac{1}{2} (z_i -a_{\alpha}g_{\alpha i})C^{-1}_{ij} (z_j -a_{\beta}g_{\beta j}) \right] \exp \left[ -\frac{1}{2} (x_i -a_{\alpha}g_{\alpha i})C^{-1}_{ij} (x_j -a_{\beta}g_{\beta j}) \right] $$ where $C_{ij}$ is the symmetric $N$ X $N$ covariance matrix, $g_{\alpha i}$ is a $n$ X $N$ matrix, repeated indexes are summed over, and $n\le N$.
If $N=n$ and $\det g \neq 0$, then one can change variables and the integral becomes a simple convolution (and one sums the covariance matrixes and means). And if $n<N$?
Not with this notation. If $M^2$ is a positive definite matrix then $$\int_{\mathbb{R}^n} e^{-\pi (\mu- G x)^T M^2 (\mu- G x) } e^{-2 i \pi \langle x, \xi\rangle} d^n x$$ is directly related to the Fourier transform of the Gaussian, thanks to the change of variable $$y = M (\mu -G x), \qquad d^n y = |\det MG|d^n x$$
By analytic continuation, the obtained formula stays true for $\xi \in \mathbb{C}^n$.
Note $e^{-\pi t^2}$ and hence $e^{-\pi \|x\|^2}$ is its own Fourier transform.