Accelerating Conjugate Gradients fitting for small localized kernel (like cubic B-spline)

Question:

Is there some pre-conditioner for Conjugate-Gradient (CG) cheap enough, that it is worth using even if my operator is very local (i.e. already has a low number of non-zero elements), as it is when fitting B-spline functions?

or any other way how to improve the speed of CG convergence?

Background:

I fit function described on 2D/3D grid using localized basis function (e.g. B-splines) using Conjugate Gradient (CG) algorithm.

The least-square fitting reduce to soltion of this matrix problem:

$$(B^T B) {\vec x}= B {\vec b}$$

where $$B$$ is overlap-matrix between the basis function. Iterative method like CG use matrix-vector multiplications:

$$\vec{y} = (B^T B) \vec{x}$$

which in my case reduces into 2D/3D convolution with a stencil formed by tensor-product of 1D Cubic Spline (i.e. 3x3 or 3x3x3 tensor in 2D/3D case)

e.g. for 2D case:

$$[ [1,4,1], [4,16,4], [1,4,1], ]$$

You can see that this matrix is very sparse and relatively cheap (although not so cheap if the grid has ~1 million points).

Is there any preconditioner other than Diagonal Jacobi worth using in such a situation?

Currently CG takes 50 iterations ot achieve $$10^{-6}$$ preccision:

CG[1] err 0.126261
CG[2] err 0.052267
CG[3] err 0.0294687
CG[4] err 0.0193294
...
CG[47] err 4.449e-05
CG[48] err 4.24922e-05
CG[49] err 3.60459e-05
CG[50] err 3.32389e-05


Just in case code can be found here (1,2,3) (but not sure if it is easy to understand)