The industry standard way of doing this is is to calculate a "mask" of sorts, and apply it to each one.
By "mask", I mean, say, a function for the current point based on the points before and after. For example:
$f(x) = \frac{1}{9}[x-2] + \frac{2}{9}[x-1] + \frac{3}{9}[x] + \frac{2}{9}[x+1] + \frac{1}{9}[x+2] $
where $[p]$ is the value of the pixel at $p$.
So to find the new pixel value at, say, pixel $4$, you'd use $ \frac{1}{9}[2] + \frac{2}{9}[3] + \frac{3}{9}[4] + \frac{2}{9}[5] + \frac{1}{9}[6] $ (remembering that, say, $[3]$ is the value of pixel 3)
All that's left to do is then apply your mask to every pixel. (What should you use for pixels that are close to the edge? That's up to you.)
Note that to be a true blurring mask, your coefficients must add up to 1.
A "Gaussian Blur" is just applying a special mask based off of the Gaussian Curve to your pixels.
That is, you make the "coefficient" of each term a number based off of the Gaussian Curve at that point.
For example, wikipedia lists:
$f(x) = 0.00038771[x-3] + 0.01330373[x-2] +
0.11098164[x-1] + 0.22508352[x] + 0.11098164[x+1] +
0.01330373[x+2] + 0.00038771[x+2]$
as a mask with σ = 0.84089642 and a "cut-off" radius of 3, although you can pick your own standard deviation and cut-off radius to your liking.
Your first expression,
$G(x,\sigma_1^2) G(y,\sigma_2^2) = \frac{1}{\sqrt{2\pi}\sigma_1}e^{\frac{-x^2}{2\sigma_1^2}}\frac{1}{\sqrt{2\pi}\sigma_2}e^{\frac{-y^2}{2\sigma_2^2}}$ = $\frac{1}{2\pi\sigma_1\sigma_2}e^{-(\frac{x^2}{2\sigma_1^2}+\frac{y^2}{2\sigma_2^2})}$
is a correct anisotropic 2D Gaussian. The expression you are trying to obtain,
$\frac{1}{2\pi(\sigma_1^2 + \sigma_2^2)}e^{\frac{-(x^2+y^2)}{2\sigma_1^2 + 2\sigma_2^2}}$
is not (it’s an isotropic Gaussian). There is no way you can go from the one to the other. However, one can incorrectly match the two expressions by creating a new coordinate system that is a scaled version of the original one, in such a way that the Gaussian is isotropic in the new coordinate system. The other answer does this. But that is not useful.
It seems that OP is confused about the combination of the variances. The convolution of two 1-dimensional Gaussian functions with variances $\sigma_1^2$ and $\sigma_2^2$ is equal to a 1-dimensional Gaussian function with variance $\sigma_1^2 + \sigma_2^2$. There is no such expectation for the multiplication of Gaussians (in fact, when multiplying them, assuming the same orientation and the same mean, the resulting variance is smaller, not larger.
Best Answer
It can be shown using some simple convolution theory. First, recall that convolution is associative: $$ f\ast(g\ast h) = (f\ast g )\ast h $$
Next recall that Gaussian blurring an image $I$ is simply convolving a Gaussian kernel $G$ to it, where $$ G(x,y|\sigma)=(2\pi\sigma^2)^{-1}\exp\left( -\frac{x^2+y^2}{2\sigma^2} \right) $$ So Gaussian blurring twice is equivalent to convolving twice, so $$ I_B=G_1\ast (G_2\ast I)=(G_1\ast G_2)\ast I = G\ast I $$ where we know that $G$ is a Gaussian kernel because the convolution of two Gaussians is a Gaussian.
Now we just need to show that $$ G(x,y|\sigma) = G\left(x,y\;\middle|\sqrt{\sigma_1^2+\sigma_2^2}\right) = G_1(x,y|\sigma_1)\ast G_2(x,y|\sigma_2) $$ One way to do this is by definition: just compute $$ G(x,y|\sigma)=\iint_{-\infty}^\infty G_1(\tau,\xi|\sigma_1)\ast G_2(x-\tau,y-\xi|\sigma_2)\, d\tau\, d\xi $$ which will eventually be equal to the desired result (see e.g. here to see it done).
However, there is an easier way using some simple probability theory: Recall that the sum of two independent random variables gives a random variable with a density equal to the convolution of the two randomc variables.
So, if $A\sim\mathcal{N}(\mu_A,\sigma_A^2)$, $B\sim\mathcal{N}(\mu_B,\sigma^2_B)$ are independent, then $C = A+B$ has $C\sim\mathcal{N}(\mu_A+\mu_B,\sigma^2_A+\sigma^2_B)$.
The multivariate generalization is also true: Notice that if $X_1\sim \mathcal{N}(0,\sigma_1^2I_2)$, $X_2\sim \mathcal{N}(0,\sigma_2^2I_2)$, then they have density functions $G_1$ and $G_2$ respectively. Thus, the sum $Z = X_1 + X_2$ has a density function given by $G=G_1 \ast G_2$. But we know that $Z\sim \mathcal{N}(0+0, \sigma^2_1I_2+\sigma^2_2I_2)=\mathcal{N}(0, (\sigma^2_1 +\sigma^2_2)I_2)$. Thus, the density function of $Z$ is given by: \begin{align*} p_Z(z) &= \frac{1}{\sqrt{4\pi^2|\Sigma|}}\exp\left( -\frac{1}{2}(z-0)^T\Sigma^{-1}(z-0) \right) \\ &= \frac{1}{2\pi (\sigma^2_1 +\sigma^2_2) }\exp\left( -\frac{1}{2}\frac{z^Tz}{[\sigma_1^2+\sigma_2^2]} \right) \\ &= \frac{1}{2\pi (\sigma^2_1 +\sigma^2_2) }\exp\left( -\frac{x^2+y^2}{2[\sigma_1^2+\sigma_2^2]} \right) \\ &= G\left(x,y\;\middle|\sqrt{\sigma_1^2+\sigma_2^2}\right) \\ &=: G(x,y|\sigma) \end{align*} where $z=(x,y)$ and $|\Sigma|=|(\sigma^2_1 +\sigma^2_2)I_2|=(\sigma^2_1 +\sigma^2_2)^2$.
This shows that $ \sigma = \sqrt{\sigma_1^2+\sigma_2^2} $, as desired.
See also: [1], [2], [3], [4], [5].