Sparse matrix representation of a spline interpolation

I use spline interpolation within a statistical model, and the transpose of the operator turns up in the gradient of the log-likelihood. Let me set up some notation first. If $x_1 \ldots x_n$ are a set of measurement locations, $f_1 \ldots f_n$ the value of $f(x)$ at these locations, then spline interpolation of $f$ at another set of locations is a linear operator in ${\bf f} = \left[f_1 \ldots f_n\right]$, which I'll call ${\bf A}$. Since splines have compact support ${\bf A}$ is generally quite sparse.

I'm looking for a fast way to compute a representation of ${\bf A}^t$ and ${\bf A}^t{\bf A}$ as sparse matrices. I do have a way to do this in R, but the sparsity pattern is lost so there must be a better way. Here's my current method: let g(x) be the interpolant. We have that $g(x) = \sum{\alpha_i b_i(x)}$, where the $b_i(x)$'s are the B-spline basis functions. We can solve for the weights by imposing that: $g(x_i) = f_i$ for all $i$, and the boundary conditions $g''(x_1) = 0$ and $g''(x_n) = 0$. These are linear constraints on ${\bf \alpha}$, which we can write as ${\bf M}\alpha = {\bf I^+f}$, with ${\bf I^+} = \left[\begin{array}{cc} \mathbf{I} & \mathbf{0}_{2\times1}\\ \mathbf{0}_{1\times2} & \mathbf{0}_{2\times2} \end{array}\right]$.

Putting everything back together we see that ${\bf A}$ can be written as ${\bf B\left(M^{-1}I^+ f\right)}$, where $\bf B$ is the B-spline basis matrix at the interpolated locations. Here's some R code:

interpmat <-  function(x,xi)
{
require(fda)
#Compute the B-spline basis functions at x with breaks at the measurement locations
V <- bsplineS(x,x,returnMat=T)
#The second derivative at the first and last points. We need that to impose boundary counditions (zero second derivative beyond the endpoints)
D <- bsplineS(c(x,x[length(x)]),x,returnMat=T,nderiv=2)
M <- rBind(V,D)
B <- bsplineS(xi,x,returnMat=T) #The B-spline basis at the interpolation locations
(B%*%solve(M))[,1:length(x)]
}

The problem is that the sparsity pattern is lost in the matrix inversion step. Any clever ideas? I know that one way would be to keep track of which basis functions have support on each interpolated point, but that sounds like a lot of painstaking C code.