Sometimes, by iteration, linear algebra can be used to solve non-linear equations. For example, consider the system $$Ax=a \qquad B(x)y=b(x), $$ where $a$ is a vector with scalar entries, $A$ is a matrix with scalar entries, $b(x)$ is a vector whose entries are functions of $x$, and $B(x)$ is a matrix whose entries are functions of $x$. This system can be solved by first solving $Ax=b$, then substituting the solution into the second equation $By=b$, and then solving the second equation. The system can also be solved by first solving $By=b$ over the ring of functions of $x$, and then solving the first equation.
Similarly, formal$^*$ Gaussian integration techniques can sometimes be used iteratively to compute the exact integrals of non-Gaussian integrands. Here's a 3D example in the variables $a,b,x$; it is easy to raise this example to higher dimensions by replacing scalars with vectors and matrices.
Let $L=\lambda ab+\frac12 q(b)x^2+\alpha a+\beta b+\xi x$, where all the letters represent scalars except for $q(b)$ which is a function of $b$. We wish to compute $I := \int e^Lda\,db\,dx$. This is not a Gaussian integral because the $q(b)x^2$ term is not quadratic in the integration variables.
Yet first computing the $ab$ integral we get $$ I(x) := \int e^Lda\,db = e^{\xi x}e^{q(\partial_\beta)x^2/2}\int e^{\lambda ab+\alpha a+\beta b}da\,db $$ $$ = \frac{2\pi}{\lambda}e^{\xi x}e^{q(\partial_\beta)x^2/2}e^{-\alpha\beta/\lambda} = \frac{2\pi}{\lambda}e^{-\alpha\beta/\lambda + \xi x + q(-\alpha/\lambda)x^2/2}. $$
Thus $I(x)$ is a Gaussian with respect to $x$, so we can (formally) compute $$ I = \int I(x)dx = \frac{(2\pi)^{3/2}}{\lambda\sqrt{q(-\alpha/\lambda)}} e^{-\alpha\beta/\lambda - q(-\alpha/\lambda)^{-1}\xi^2/2}. $$
We could have arrived at the same result by first computing the $x$ integral as a formal Gaussian over the ring of functions of $b$ and then computing the $ab$ integral.
Question. Is there a precedent for this procedure? A name? Is there a place where people routinely iterate Gaussian integration to integrate non-Gaussians?
$^*$Meaning, applying standard formulas without worrying about convergence. Add conditions if you must, or think that we're really imitating some QFT-like context in which convergence is not an issue.
