  Functions on a bounded domain in scientific computing are often approximated
using piecewise polynomial approximations on meshes that adapt to the shape of
the geometry. We study the problem of function approximation using splines on a
regular but oversampled grid that is defined on a bounding box. This approach
allows the use of high order and highly structured splines as a basis for
piecewise polynomials. The methodology is analogous to that of Fourier
extensions, using Fourier series on a bounding box, which leads to spectral
accuracy for smooth functions. However, Fourier extension approximations
involve solving a highly ill-conditioned linear system, and this is an
expensive step. The computational complexity of recent algorithms is $\mathcal
O\left(N\log^2\left(N\right)\right)$ in 1-D and $\mathcal
O\left(N^2\log^2\left(N\right)\right)$ in 2-D. We show that, compared to
Fourier extension, the compact support of B-splines enables improved complexity
for multivariate approximations, namely $\mathcal O(N)$ in 1-D, $\mathcal
O\left(N^{3/2}\right)$ in 2-D and more generally $\mathcal
O\left(N^{3(d-1)/d}\right)$ in $d$-D with $d>1$. By using a direct sparse QR
solver for a related linear system, we also observe that the computational
complexity can be nearly linear in practice. This comes at the cost of
achieving only algebraic rates of convergence. Our statements are corroborated
with numerical experiments and Julia code is available.
