# Change of variables in a Gaussian integral in matrix form

I have a problem in which I have to compute the following integral: $$\mathop{\idotsint\limits_{\mathbb{R}^k}}_{\sum_{i=1}^{k}y_i=x} e^{-N^2r(\sum_{i=1}^{k}y_i^2-\frac{1}{k}x^2)} dy_1\dots dy_{k},$$ where this notation means that I want to integrate over $$\mathbb{R}^k$$ restricted to the plain where $$\sum_{i=1}^{k}y_i=x$$ (a convolution of gaussians) and $$N$$ and $$r$$ are positive real constants. I have tried two different methods for computing this integral, but they are yielding different results. I would appreciate it very much if someone could take a look and tell me what I'm doing wrong.

Method 1

In method 1 I just wrote it as $$\mathop{\idotsint\limits_{\mathbb{R}^k}}_{\sum_{i=1}^{k}y_i=x} e^{-N^2r(\sum_{i=1}^{k}y_i^2-\frac{1}{k}x^2)} dy_1\dots dy_{k}=\int_{-\infty}^{\infty}\dots\int_{-\infty}^{\infty} e^{-N^2r((x-y_1)^2+\sum_{i=1}^{k-2}(y_i-y_{i+1})^2+y_{k-1}^2-\frac{1}{k}x^2)} dy_1\dots dy_{k-1}=\sqrt{\frac{1} {\pi r^{k-1}k}}\frac{\pi^{k}}{N^{k-1}}$$

I deduced this formula by induction, first integrating in $$y_{k-1}$$, then $$y_{k-2}$$ and so on.

Method 2

In method 2 I tried writting the function in a matrix form $$\mathop{\idotsint\limits_{\mathbb{R}^k}}_{\sum_{i=1}^{k}y_i=x} e^{-N^2r(\sum_{i=1}^{k}y_i^2-\frac{1}{k}x^2)} dy_1\dots dy_{k}=\mathop{\idotsint\limits_{\mathbb{R}^k}}_{\sum_{i=1}^{k}y_i=x} e^{-N^2r(\vec{y},Q\vec{y})} dy_1\dots dy_{k}$$ where $$$$Q:=\left(\begin{array}{cccccccc} (1-\frac{1}{k})& -\frac{1}{k} & -\frac{1}{k} & \cdots & -\frac{1}{k} \\ -\frac{1}{k} & (1-\frac{1}{k}) & -\frac{1}{k} & \cdots & -\frac{1}{k} \\ \vdots & \ddots & & &\vdots \\ -\frac{1}{k} & \dots & &-\frac{1}{k} &(1-\frac{1}{k}) \end{array}\right)$$$$.

This matrix $$Q$$ has eigenvalues $$\lambda_0=0$$, $$\lambda_l=1$$ and corresponding normalized eigenvetors $$$$\vec{\lambda}_l=\frac{1}{\sqrt{k}}\left(\begin{array}{c} 1 \\ e^{\frac{2\pi i}{k}1l} \\ \vdots \\ e^{\frac{2\pi i}{k}(k-1)l} \\ \end{array}\right)$$$$ for $$0\le l\le k-1$$.

As I understand it, the restriction in the integral means that I shouldn't integrate in the $$\lambda_0$$ direction, since in this direction I must have all components equal, and the only place where the components are equal and the bound is satisfied is $$(\frac{x}{k},\dots,\frac{x}{k})$$. So my integration should occour in the orthogonal complement of this vector, which is a hyperplane of dimension $$k-1$$. Everything seems to check to this point, so I diagonalized the matrix $$Q=U\Lambda U^{-1}$$ and so

$$(\vec{y},Q\vec{y})=(\vec{\xi},\Lambda\vec{\xi})=\sum_{i=1}^{k-1}\xi_i^2.$$

The change of variables $$\vec{\xi}=U^{-1}\vec{y}$$ has a Jacobian $$\frac{1}{\sqrt{k^{k-1}}}$$, since $$U^{-1}$$ is the DFT matrix times $$\frac{1}{\sqrt{k^{k-1}}}$$ and the DFT matrix is known to be unitary. So

$$\mathop{\idotsint\limits_{\mathbb{R}^k}}_{\sum_{i=1}^{k}y_i=x} e^{-N^2r(\vec{y},Q\vec{y})} dy_1\dots dy_{k}=\idotsint\limits_{\mathbb{R}^k} e^{-N^2r\sum_{i=1}^{k-1}\xi_i^2} \frac{1}{\sqrt{k^{k-1}}}d\xi_1\dots d\xi_{k-1}= \sqrt{\frac{\pi^{k-1}}{k^{k-1}r^{k-1}}}\frac{1}{N^{k-1}}$$.

These two results are different and I cannot figure out why.

• In your Method 2 the quantity $(\vec{y},Q\vec{y})$ is supposed to be equal to $\sum_i y_i^2-x^2/k$? I don't see this – thedude Jan 14 at 13:10
• It's because $(\vec{y},Q\vec{y})=\sum_{i=1}^{k}y_i^2-\frac{1}{k}(\sum_{i=1}^{k}y_i)^2$ and we are subject to the bound $\sum_{i=1}^{k}y_i=x$ – Rafael Jan 14 at 14:46

$$\newcommand\R{\mathbb R}\newcommand\1{\mathbf1}$$When you say "I want to integrate over $$\mathbb{R}^k$$ restricted to the plain where $$\sum_{i=1}^{k}y_i=x$$", you have to specify the measure over the plane over which you want to integrate.

It appears you want this measure to be induced by the Lebesgue measure on $$\R^k$$. Then the integration can be done as follows. Let $$c:=N^2r\in(0,\infty)$$ and $$t:=x/\sqrt k,$$ the (signed) distance from the the origin to your plane $$\Pi_t:=\{y\in\R^k\colon u\cdot y=t\}=\{y\in\R^k\colon \1\cdot y=x\},$$ where $$\cdot$$ denotes the dot product, $$\1:=(1,\dots,1)\in\R^k$$, and $$u:=\1/\sqrt k$$ is a unit normal vector to the plane $$\Pi_t$$. Thus, instead of the parameter $$x$$, we use the more geometrical parameter $$t$$.

Then the integral in question can be written as $$I_t:=e^{ct^2}J_t,\quad\text{where}\quad J_t:=\int_{\Pi_t}\mu_t(dy)e^{-c|y|^2},$$ $$|y|$$ is the Euclidean norm of $$y$$, and, for each real $$t$$, $$\mu_t$$ is the measure over the plane $$\Pi_t$$ induced by the Lebesgue measure on $$\R^k$$ in the following sense: $$$$\int_a^b dt\, \int_{\Pi_t}\mu_t(dy)g(y) =\int_{\Pi_{a,b}}dy\,g(y) \tag{1}$$$$ for all nonnegative Borel-measurable functions $$g\colon\R^k\to\R$$ and all real $$a$$ and $$b$$ such that $$a, where $$\Pi_{a,b}:=\bigcup_{t\in[a,b]}\Pi_t=\{y\in\R^k\colon a\le u\cdot y\le b\}.$$

Then for such $$a$$ and $$b$$ we have $$\int_a^b dt\, J_t=\int_a^b dt\, \int_{\Pi_t}\mu_t(dy)e^{-c|y|^2} =K_{a,b}:=\int_{\Pi_{a,b}}dy\,e^{-c|y|^2}.$$ (See the remark on this at the end of this answer.)

To compute the integral $$K_{a,b}$$, let us use a substitution of the form $$y=Qz$$, where $$Q$$ is any orthogonal $$k\times k$$ matrix whose first column is the unit vector $$u$$, so that $$y=Qz$$ implies $$z_1=u\cdot y$$, where $$z_j$$ is the $$j$$'s coordinate of the vector $$z$$; such an orthogonal matrix $$Q$$ exists. Then we can write $$K_{a,b}=\int_{\R^k}dy\,e^{-c|y|^2}\,1(a\le u\cdot y\le b) \\ =\int_{\R^k}dz\,e^{-c|z|^2}\,1(a\le z_1\le b) \\ =\int_a^b dz_1\,e^{-cz_1^2}\int_{\R^{k-1}}dw\,e^{-c|w|^2} \\ =\int_a^b dz_1\,e^{-cz_1^2}\,(\pi/c)^{(k-1)/2}.$$ So, $$J_t=\frac d{dt}\,K_{a,t}=e^{-ct^2}\,(\pi/c)^{(k-1)/2}.$$ Thus, the integral in question is $$I_t:=e^{ct^2}J_t=(\pi/c)^{(k-1)/2}=\Big(\frac\pi{N^2r}\Big)^{(k-1)/2}.$$

This differs from both of your answers -- but you never defined the measure over which you integrate.

Remark: Intuitively, think of $$a$$ and $$b$$ as being close to a real number $$t$$, and hence to each other. We approximate the integral $$\int_{\Pi_t}\mu_t(dy)e^{-c|y|^2}$$ over the plane $$\Pi_t$$ by $$\frac1{b-a}\,\int_{\Pi_{a,b}}dy\, e^{-c|y|^2}$$, that is, by the integral over the thin layer $$\Pi_{a,b}$$ between two parallel planes $$\Pi_a$$ and $$\Pi_b$$ (close to the plane $$\Pi_t$$) divided by thickness $$b-a$$ of the layer.

Formally, we are dealing here with disintegration of a measure. That linked theorem deals only with probability measures, but it is trivially extended to finite measures. If we forget about this finiteness condition for a moment, then in that linked theorem we can choose $$X=\R$$, $$Y=\R^k$$, let the map $$\pi\colon Y\to X$$ be the projection map defined by $$\pi(y):=u\cdot y$$ for all $$y\in Y=\R^k$$, and let $$\mu:=\lambda_k$$ and $$\nu:=\lambda_1$$, where $$\lambda_k$$ is the Lebesgue measure over $$\R^k$$. I think the finiteness condition in that linked theorem is inessential, and the proof will hold for any Borel measures $$\mu$$, at least if $$\mu$$ is $$\sigma$$-finite. Alternatively, one can approximate here the Lebesgue measure over $$\R^k$$ by the finite Lebesgue measures over big cubes in $$\R^k$$.

• I'm confused now, because this also differs from the answer form @CarloBeenakker. I indeed want to integrate over a measure induced by the Lebesgue measure, but do you know why all the answers are different? The presence or absence of the factor $\frac{1}{\sqrt{k}}$ is importante for me. – Rafael Jan 14 at 17:18
• @Rafael : As I said, the other answers are different because they deal with integrals over undefined measures. Therefore, those integrals were undefined, and thus you cannot possibly attach any value to them. – Iosif Pinelis Jan 14 at 17:56
• I'm not sure I understand. What would be, for example, the measure from the below example and my first method and how does it differ from the Lebesgue measure? What is then the measure from the second method? And thank you very much for the comment, your answers always help a lot. – Rafael Jan 14 at 18:26
• @Rafael : As I said, I do not see any measures on the plane defined in the other answers. Also, I do not see good ways to introduce measures in those answers. – Iosif Pinelis Jan 14 at 18:35
• @Rafael : I have added a remark that should help one understand how the measures $\mu_t$ are constructed. – Iosif Pinelis Jan 14 at 19:13

The desired integral $$I$$ can be written as $$^\ast$$

$$I=e^{N^2rx^2/k}J(x),\;\;J(x)=\mathop{\idotsint}\delta\left(x-{\textstyle{\sum_{j=1}^{k}}y_j}\right) e^{-N^2r\sum_{j=1}^{k}y_j^2} dy_1\dots dy_{k}.$$ Fourier transform $$J(x)$$, $$F(q)=\int_{-\infty}^\infty e^{iqx}J(x)\,dx=\mathop{\idotsint} e^{-N^2r\sum_{j=1}^{k}y_j^2}e^{\sum_{j=1}^k qy_j} dy_1\dots dy_{k}=N^{-k}(\pi/r)^{k/2}e^{-kq^2/(4N^2r)}.$$ Now Fourier transform back from $$F(q)$$ to $$J(x)$$ and you're done: $$J(x)=(2\pi)^{-1}\int_{-\infty}^\infty e^{-iqx}F(q)\,dq=k^{-1/2} (N^2 r/\pi)^{(1-k)/2} e^{-N^2 r x^2/k}$$ $$\Rightarrow I=k^{-1/2} (N^2 r/\pi)^{(1-k)/2}.$$ This differs from both of the answers in the OP, but answer number 1 will agree if a typo is corrected ($$\pi^k$$ should be $$\pi^{k/2}$$, which is the factor coming from the Gaussian integrals in method 1).

$$^\ast$$ As noted by Iosif Pinelis, the notation $$\idotsint_{\sum_{i=1}^{k}y_i=x}$$ in the OP is ambiguous, I replaced it by the delta function $$\idotsint \,\delta(x-\sum_{i=1}^{k}y_i)$$ to remove the ambiguity.
More generally, one can define $$I_a$$ and $$J_a(x)$$ with the delta function $$\delta(ax-a\sum_{i=1}^{k}y_i)$$ for any $$a>0$$. This amounts to the same Euclidian measure on the hyperplane, but the distance from the hyperplane is rescaled by a factor $$a$$. Since $$\delta(au)=a^{-1}\delta(u)$$, one has $$I_a=a^{-1}I_1(x)=a^{-1}k^{-1/2}(N^2r/\pi)^{(1-k)/2}.$$ Iosif's answer corresponds to the choice $$a=k^{-1/2}$$.

• I'm going to check my iterations for this mistake! Thanks! It is good to know, that method 1 works. Do you have any idea what went wrong with method 2? – Rafael Jan 14 at 14:50
• @Rafael --- you asked "how do you define the delta function measure?" --- you can take a look on Wikipedia --- the only property of the delta function is need is the identity $\int_{-\infty}^\infty f(x)\delta(x-u)dx=f(u)$. – Carlo Beenakker Jan 15 at 8:05
• @CarloBeenakker : I understand that the ordinary integral $\int_{-\infty}^\infty f(x)\delta(x-u)dx$ is defined as $f(u)$. But how is your first $k$-fold integral defined? – Iosif Pinelis Jan 15 at 18:50
• @IosifPinelis --- the same identity applies there as well, you can use it to integrate out one of the $y_j$'s, for example $y_1$: $$\mathop{\idotsint}\delta\left(x-{\textstyle{\sum_{j=1}^{k}}y_j}\right) e^{-N^2r\sum_{j=1}^{k}y_j^2} dy_1\dots dy_{k}=$$ $$\mathop{\idotsint}e^{-N^2r(x-\sum_{j=2}^k y_j)^2}e^{-N^2r\sum_{j=2}^{k}y_j^2}\,dy_2\dots dy_{k}.$$ – Carlo Beenakker Jan 15 at 20:33
• Thank you very much, I guess this is the integral I was looking for – Rafael Jan 16 at 22:40

$$\newcommand\1{\mathbf1}\newcommand{\R}{\mathbb R}\newcommand{\la}{\lambda}$$Here is a more explicit way to define the disintegration of the Lebesgue measure over $$\R^k$$ into the measures $$\mu_t$$ over the planes $$\Pi_t$$ introduced in my other answer on this page, of course with the same result.

Let us recall notations introduced in that answer: $$c:=N^2r\in(0,\infty),\quad t:=x/\sqrt k,$$ $$\Pi_t:=\{y\in\R^k\colon u\cdot y=t\}=\{y\in\R^k\colon \1\cdot y=x\},$$ where $$\cdot$$ denotes the dot product, $$\1:=(1,\dots,1)\in\R^k$$,
$$u:=\1/\sqrt k;$$

The integral in question still is $$\begin{equation*} I_t:=e^{ct^2}J_t,\quad\text{where}\quad J_t:=\int_{\Pi_t}\mu_t(dy)e^{-c|y|^2}, \tag{1} \end{equation*}$$ $$|y|$$ is the Euclidean norm of $$y$$, and (for each real $$t$$) $$\mu_t$$ is the measure over the plane $$\Pi_t$$ now explicitly defined as follows.

Let $$T\colon\R^{k-1}\to\Pi_0$$ be any linear isometry of $$\R^{k-1}$$ onto the $$(k-1)$$-dimensional linear subspace $$\Pi_0$$ of $$\R^k$$. Any such isometry $$T$$ is given by the formula $$Tz=\sum_{j=1}^{k-1}z_jb_j$$ for all $$z=(z_1,\dots,z_{k-1})\in\R^{k-1}$$, where $$(b_1,\dots,b_{j-1})$$ is any orthonormal basis of $$\R^{k-1}$$. For each real $$t$$, we have $$\Pi_t=\Pi_0+tu$$, and so, we can define the affine isometry $$U_t\colon\R^{k-1}\to\Pi_t$$ of $$\R^{k-1}$$ onto the $$(k-1)$$-dimensional affine subspace $$\Pi_t$$ by the formula $$\begin{equation*} U_tz:=Tz+tu \end{equation*}$$ for all $$z\in\R^{k-1}$$. Let now $$\begin{equation*} \mu_t:=\mu^T_t:=\la_{k-1}U_t^{-1}, \tag{2} \end{equation*}$$ the pushforward measure for the Lebesgue measure $$\la_{k-1}$$ over $$\R^{k-1}$$ under the isometry $$U_t$$, so that $$\mu_t(B)=\la_{k-1}(U_t^{-1}(B))=\la_{k-1}(T^{-1}(B-tu))$$ for all Borel subsets $$B$$ of $$\Pi_t$$.

Remark: The measures $$\mu_t=\mu^T_t$$ do not depend on the choice of a linear isometry $$T$$ of $$\R^{k-1}$$ onto $$\Pi_0$$. Indeed, if $$T$$ is any such isometry, then any other such isometry is of the form $$S:=TQ$$, where $$Q$$ is any linear isometry of $$\R^{k-1}$$ (onto $$\R^{k-1}$$). Hence, $$\begin{equation*} \mu^S_t(B)=\la_{k-1}(S^{-1}(B-tu))=\la_{k-1}(Q^{-1}(T^{-1}(B-tu))) =\la_{k-1}(T^{-1}(B-tu))=\mu^T_t(B) \end{equation*}$$ for all Borel subsets $$B$$ of $$\Pi_t$$; the penultimate equality in the above display holds because the Lebesgue measure is invariant with respect to linear isometries. $$\Box$$

It follows from (1) and (2) that $$$$J_t=\int_{\R^{k-1}}dz\,\exp\{-c|U_tz|^2\},$$$$ by the change-of-variables formula for the pushforward measures. Since $$U_tz=Tz+tu$$, $$Tz\perp u$$, $$T$$ is an isometry, and $$|u|=1$$, we have $$|U_tz|^2=|Tz|^2+t^2|u|^2=|z|^2+t^2$$. So, $$$$J_t=e^{-ct^2}\int_{\R^{k-1}}dz\,\exp\{-c|z|^2\} =e^{-ct^2}(\pi/c)^{(k-1)/2}.$$$$ Thus, the integral in question is $$I_t=e^{ct^2}J_t=(\pi/c)^{(k-1)/2}=\Big(\frac\pi{N^2r}\Big)^{(k-1)/2},$$ which is what was obtained a bit differently in my other answer.

$$\newcommand\1{\mathbf1}\newcommand{\R}{\mathbb R}\newcommand{\la}{\lambda}$$Here is yet another, "multivariate calculus" way to treat your integral, with the same result as in my other two answers on this page.

Let $$c:=N^2r\in(0,\infty),$$ $$\Pi_x:=\{y\in\R^k\colon \1_k\cdot y=x\},$$ where $$\cdot$$ denotes the dot product and $$\1_k:=(1,\dots,1)\in\R^k$$.

The integral in question is $$\begin{equation*} I_x:=e^{cx^2/k}J_x,\quad\text{where}\quad J_x:=\int_{\Pi_x}\mu_x(dy)e^{-c|y|_k^2}, \tag{1} \end{equation*}$$ $$|y|_k$$ is the Euclidean norm of $$y\in\R^k$$, and $$\mu_x(dy)$$ is the surface area element on the plane $$\Pi_x$$. This plane is the graph of the function $$f_x\colon\R^{k-1}\to\R$$ given by the formula $$f_x(y):=x-y_1-\cdots-y_{k-1}$$ for $$y=(y_1,\dots,y_{k-1})\in\R^{k-1}$$. So (cf. the case $$k=3$$), \begin{align*} J_x&=\int_{\R^{k-1}}dy_1\cdots dy_{k-1}\,\sqrt{1+\sum_{j=1}^{k-1}\Big(\frac{\partial f_x}{\partial y_j}\Big)^2} \\ & \times \exp\Big\{-c(x-y_1-\cdots-y_{k-1})^2-c\sum_{j=1}^{k-1}y_j^2\Big\} \\ &=\sqrt k\,\int_{\R^{k-1}}dy\, \exp\big\{-c(x-y\cdot\1_{k-1})^2-c|y|_{k-1}^2\big\} \\ &=\sqrt k\,\int_{\R^{k-1}}dz\, \exp\big\{-c(x-z_1\sqrt{k-1})^2-c|z|_{k-1}^2\big\}; \end{align*} for the last displayed equality, we use any substitution of the form $$y=Qz$$, where $$Q$$ is any orthogonal $$(k-1)\times(k-1)$$ matrix whose first column is $$\1_{k-1}/\sqrt{k-1}$$ and $$z=(z_1,\dots,z_{k-1})$$, so that $$z_1=y\cdot\1_{k-1}/\sqrt{k-1}$$.

So, \begin{align*} J_x&=\sqrt k\,\int_\R dz_1\, \exp\big\{-c(x-z_1\sqrt{k-1})^2-cz_1^2\big\} \\ &\times\int_{\R^{k-2}}dz_2\cdots dz_{k-1}\, \exp\Big\{-c\sum_{j=2}^{k-1}z_j^2\Big\} \\ &=\sqrt k\ \;\frac{e^{-cx^2/k}(\pi/c)^{1/2}}{\sqrt k}\ \;(\pi/c)^{(k-2)/2} \\ &=e^{-cx^2/k}\;(\pi/c)^{(k-1)/2}. \end{align*}

Thus, by (1), the integral in question is $$I_x=(\pi/c)^{(k-1)/2}=\Big(\frac\pi{N^2r}\Big)^{(k-1)/2},$$ which is what was obtained a bit differently in my other two answers.