What is the most efficient way to perform Gaussian blurs for varying sizes of Gaussian kernels in an image?
I'm not really used to thinking in texture fetches as I usually code high level non-hardware-specific signal processing, but low pass filtering also removes information you can couple it with sub-sampling which can give further gains. It is not clear to me if you have already done this.

In practice this means that after the first low-pass filtering, you can disregard $k$ out of $N$ of the samples in the subsequent low-pass filtering.

Simple example is the $2^n \times 2^n$ averaging-filter, first we can do 2x2 mean values and store. Then for the next pass, we only need to use the output of one of each $2 \times 2$ blocks. So we only need consider $1$ in $4^k$ data points at step $k$.

A 8x8 block this would mean:

$$I_1 = \left[\begin{array}{cccccccc}
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1\\
1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1&amp;1
\end{array}\right]$$
First iteration, all are involved in 2x2 filtering. But they are only stored at the pixels which are ones here (so we only need to calculate 1/4 of them):
$$I_2 = \left[\begin{array}{cccccccc}
1&amp;0&amp;1&amp;0&amp;1&amp;0&amp;1&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
1&amp;0&amp;1&amp;0&amp;1&amp;0&amp;1&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
1&amp;0&amp;1&amp;0&amp;1&amp;0&amp;1&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
1&amp;0&amp;1&amp;0&amp;1&amp;0&amp;1&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0
\end{array}\right]$$
Next step only involves the ones in the previous one, now our image is 4x4 instead of 8x8, so calculations are reduced a factor of 4. But we only need to store the results at pixels:
$$I_3 = \left[\begin{array}{cccccccc}
1&amp;0&amp;0&amp;0&amp;1&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
1&amp;0&amp;0&amp;0&amp;1&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0
\end{array}\right]$$
We will end up with only one out-value, in this case upper left. And here's where the interpolations come into play. We can interpolate between these upper left results of neighboring "blocks". The interpolation kernel can be dependent on the low-pass filter, but I assume that this is where the hardware shines with it's strengths, having really fast ways to do this for us?

If you are interested in theory filter-banks, sub-band coding and maybe wavelets are interesting to read up on. The example I brought up here is the implementation low-pass filter of the basic but popular discrete Haar-wavelet. Gaussian pyramids in computer vision may be even more relevant. There's much to read about it so it's difficult to cram it all into just one answer.