Why isn't bilateral filtering shift-invariant?
The bilateral filter can be written as:

$$y[i, j] = \frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k, j+l]\,f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}\tag{1}$$

Shift invariance can be defined as any shift $m,n$ of the input $x$ resulting in the corresponding shift in the output $y$:

$$x[i,j] \to x[i-m,j-n]\quad\Rightarrow\quad y[i,j]\to y[i-m,j-n].\tag{2}$$

By doing the first substitution on all occcurrences of $x$ on the right side of Eq. 1 we get as the output:

$$\begin{gather}\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k-m, j+l-n]\,f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}\\
=\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i-m+k, j-n+l]\,f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)},\end{gather}\tag{3}$$

where we can recognize $i-m$ and $j-n$ as $i$ and $j$ of Eq. 1, so we have:

$$= y[i-m, j-n],\tag{4}$$

showing that the bilateral filter is shift-invariant.

However, these are shifts with integer $m$ and $n$. The non-linearity of the bilateral filter will make it not invariant to fractional shifts carried out by (sinc) interpolation. Perhaps this is what your source was referring to.