This is my first time answering on stack exchange.

**Proof of** : $\int f(\vec x) \, \exp\left( - \frac 1 2 \sum_{i,j=1}^{n}A_{ij} x_i x_j \right) d^nx
=
\sqrt{(2\pi)^n\over \det A} \, \left. \exp\left({1\over 2}\sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right)f(\vec{x})\right|_{\vec{x}=0}$

I found a silly simple proof that actually gives a generalization of the above result.

To make it all the more amusing I will give a Gauss style proof where, as Niels Abel cleverly said, like the fox, I will efface the tracks in the sand with my tail.

Consider a probability density $\mu$ that takes as input $\textbf{x}=(x_1,x_2,...,x_n)$ and define:

$$F(\textbf{y})=\int{\mu(\textbf{x})f(\textbf{x}+\textbf{y})d^n\textbf{x}}=\langle f(\textbf{x}+\textbf{y})\rangle,$$

where $\langle ... \rangle$ represents the average with respect to $\mu$. In the original problem, $\mu$ is Gaussian.

Next, we use the fact that the generator of translations $\textbf{x}\rightarrow\textbf{x}+\textbf{a}$ is given by

$$ P_\textbf{a}=\textbf{a}\cdot\nabla $$

where I used $P$ like in quantum mechanics for that extra mystique. By the definiton of what a generator is, this implies that

$$f(\textbf{x}+\textbf{y})=e^{\textbf{x}\cdot\nabla_y}f(y)$$

which is just a fancy way of writing an all order Taylor expansion.

Now the magical part,

$$\int{\mu(\textbf{x})f(\textbf{x}+\textbf{y})d^n\textbf{x}}=\int{\mu(\textbf{x})e^{\textbf{x}\cdot\nabla_y}f(\textbf{y})d^n\textbf{x}}=\langle e^{\textbf{x}\cdot\nabla_y}\rangle f(\textbf{y})=K(\nabla_\textbf{y})f(\textbf{y}),$$

where K is by definition the moment generating function.

Using the formula for the moment generating function of a Gaussian and setting $\textbf{y}=0$ gives the formula above.

**Proof of** : $$
\frac{1}{2^N N!} \left(\partial_x^\mathrm{T} A^{-1}\partial_x\right)^N x_{k_1}\cdots x_{k_{2N}} = \frac{1}{2^N N!} \sum_{\sigma \in S_{2N}}(A^{-1})_{k_{\sigma(1)}k_{\sigma(2)}} \cdots (A^{-1})_{k_{\sigma(2N-1)}k_{\sigma(2N)}} $$

The following is independent from the proof above and simply provides extra steps in the derivation of @user26872.

We start by writing explicitly the left hand side as a product of sums:

$$\left(\partial_x^\mathrm{T} A^{-1}\partial_x\right)^N= \left(\sum_{i,j}A^{-1}_{i,j}\partial_{x_i}\partial_{x_j}\right)^N$$

Next, we consider the following notation:

**The set of integers from 1 to N**: $\{1,...,N\}=[N]$
**The set of all functions from set $A$ to set $B$**: $A\rightarrow B$

The formula for a product of sums may then be written as (explanation of formula below):

$$\left(b_{1,1}+b_{1,2}+\ldots + b_{1,M}\right)\ldots\left(b_{N,1}+\ldots + b_{N,M}\right)=\prod^N_{k=1} \sum^M_{r=1} b_{k,r}=\sum_{f \, \in \, [N]\rightarrow[M]}\prod^N_k b_{k,f(k)}$$

In words, for each factor in the product (so for each index $k$), pick one term in that sum at random (so one random $\hat r$ which we call $\hat r=f(k)$) and then multiply all of the collected terms (so $\prod_k b_{k,f(k)}$). Finally, to get all of the terms of the expansion, sum over all configurations (so over all $f$).

In the case of a double sum we have:

$$\prod^N_{k=1} \sum^{M_1}_{r=1}\sum^{M_2}_{s=1} b_{k,r,s}=\sum_{\begin{gather} f_1 \, \in \, [N]\rightarrow [M_1] \\ f_2 \, \in \, [N]\rightarrow [M_2] \end{gather}}\prod^N_k b_{k,f_1(k),f_2(k)}$$

if $[M_1]=[M_2]=[M]$, then we may consider the change of notation $f_1(k)=f(1,k)$ and $f_2(k)=f(2,k)$ which leads to :

$$\prod^N_{k=1} \sum^{M}_{r=1}\sum^{M}_{s=1} b_{k,r,s}=\sum_{f \, \in \, [2]\mathrm{x}[N]\rightarrow [M]}\prod^N_k b_{k,f(1,\,k),f(2,\,k)}$$

To check that this last statement is true we show that both sums are over the same set.

Denoting Card as the cardinality of a set, recall that the cardinality of the set of all functions from set $A$ to set $B$ is $\textrm{Card}\left(B\right)^{\textrm{Card}\left(A\right)}$ and so

$$\mathrm{Card}\left[\left([N]\rightarrow[M]\right)\mathrm{x}\left([N]\rightarrow[M]\right)\right]=M^N M^N=\left(M^2\right)^N=M^{2 N}$$

Furthermore, $\mathrm{Card}\left([2]\textrm{x}[N]\right)=2 N$ such that

$$\mathrm{Card}\left[\left([2]\mathrm{x}[N]\right)\rightarrow[M]\right]=M^{2N}$$

And so both sets have equal cardinality and both sums are equal.

Notice that the cardinality $M^{2N}$ leads to a third formulation of the sum as this is also the cardinality of the set $[2N]\rightarrow[M]$. We thus stack both copies of $[N]$, that is $\left(1,[N]\right)$ and $\left(2,[N]\right)$ and consider all functions from this new stacked set to the set $[M]$. We then arrive at the crucial formula of this post:

$$\prod^N_{k=1} \sum^M_{r=1}\sum^M_{s=1} b_{k,r,s}=\sum_{f \, \in \, [2N]\rightarrow [M]}\prod^N_{k=1} b_{k,f(k),f(k+N)}$$

For the original problem we have $b_{k,r,s}=A^{-1}_{r,s}\partial_{x_r}\partial_{x_s}$ which leads to:

$$\left(\sum_{i,j}A^{-1}_{i,j}\partial_{x_i}\partial_{x_j}\right)^N=\sum_{f \, \in \, [2N]\rightarrow [M]}\prod^N_{k=1} A^{-1}_{f(k),f(k+N)}\partial_{x_{f(k)}}\partial_{x_{f(k+N)}}$$

Reordering the partial derivatives and applying this operator on the product $ x_{k_1} \cdots x_{k_{2N}} $ then leads to

$$\left(\sum_{i,j}A^{-1}_{i,j}\partial_{x_i}\partial_{x_j}\right)^N\, x_{k_1} \cdots x_{k_{2N}}=\sum_{f \, \in \, [2N]\rightarrow [M]}\prod^N_{k=1} A^{-1}_{f(k),f(k+N)}\prod^{2N}_{k=1}\partial_{x_{f(k)}}\; x_{k_1} \cdots x_{k_{2N}}$$

As $\partial_{x_i} x_j=\delta_{i,j}$, the only non zero terms of the sum are those for which

$$\left\{ f(1),...,f(2N) \right\}=\left\{ k_1,..., k_{2N} \right\}$$

so that $f$ is a bijection from the set $[2N]$ to $\left\{k_1,..., k_{2N} \right\}$. Hence, for every $f$ there exist a unique $\sigma\in S_{2N}$ such that $f(i)=k_{\sigma(i)}$. The sum over $f$ may then be replaced by a sum over $\sigma$ which leads to

$$\left(\partial_x^\mathrm{T} A^{-1}\partial_x\right)^N x_{k_1}\cdots x_{k_{2N}} = \sum_{\sigma \in S_{2N}}(A^{-1})_{k_{\sigma(1)}k_{\sigma(N+1)}} \cdots (A^{-1})_{k_{\sigma(N)}k_{\sigma(2N)}}$$

Finally, for any $\alpha\in S_{2N}$, $\sum_{\sigma}=\sum_{\sigma\circ\alpha}$, where $\circ$ represents composition, as this is simply a rearrangement of terms. Taking $\alpha$ such that

$$\left[\alpha(1),\alpha(N+1),\alpha(2),\alpha(N+2)...,\alpha(N),\alpha(2N)\right]=[1,2,3,4...,2N-1,2N]$$

leads to the desired result.