The return from this function is the tensor product of the B-splines
transformations for the given variables. Say we have variables X, Y, and Z
to build the tensor product of. The columns of the returned matrix
correspond to the column products of the three B-splines:
x1y1z1 x2y1z1 x3y1z1 x4y1z1 x1y2z1 x2y2z1 ... x4y4z4
for three fourth order B-splines with no internal knots. The columns of X
cycle the quickest, followed by Y, and then Z. This would be the same result
as
model.matrix( ~ bsplines(X) : bsplines(Y) : bsplines(Z) + 0) .
See vignette(topic = "cnr", package = "cpr") for more details.