Setting up a naive tensor product B-spline example I've simulated data according to
$y = \text{sin}(2\cdot(4x-2))+2\cdot\text{exp}(-(16^2)(x-0.5)^2)+\epsilon$
where $\epsilon \sim \mathbb{N}(0,0.3^2)$
By evaluating the ith B-Spline of degree $k$, i.e., $\text{B}_{i,k}(x)$ for $m$ inner equidistant knot intervals recursively, I`ll get a design-matrix of dimension $n\times (m+k)$ ($q=m+k$) with B-Splines as columns.
By using ordniary least squares approach it is easy to minimize $(\boldsymbol{y-X\beta})^T(\boldsymbol{y-X\beta})$ w.r.t $\boldsymbol{\beta}$.
Now I wanted to set up a similar simple example for a tensor product B-spline.
Here, the model would be $y_i = \sum_{j=0}^q\sum_{l=0}^p \beta_{j,l}\text{B}_{j,k_1}(x_i)\text{B}_{l,k_2}(z_i) + \epsilon_i\hspace{3ex} (1)$.
W.r.t $x$, the B-Spline-Matrix $\boldsymbol{B}_x$ is of dimension $n\times q$ and w.r.t. $z$, i.e., $\boldsymbol{B}_z$ is of dimension $n\times p$.
If I take the tensor product of $\boldsymbol{B}_x$ and $\boldsymbol{B}_z$, i.e., $\boldsymbol{B}_x\otimes\boldsymbol{B}_z$ this matrix will be of dimension $\Big(n^2\times (q\cdot p)\Big)$ but the vector $\boldsymbol{y}$, which will be projected onto  $\boldsymbol{B}_x\otimes\boldsymbol{B}_z$, is only of length $n$?
(1) suggest that, similar to the ordinary interaction terms, I take all possible products of those two B-Splines bases. But those two are column vectors? I would need to do this calculation for every pair, i.e., $q\cdot p$ different elementwise multiplications. But here the wording "tensor" would not make sense. 
Looking at several papers I'm still not sure how this design-matrix constructed from the tensor -product of those B-Spline matrices will look like?
 A: I know this is an old question, but it's a good one, and I thought about it recently myself.
You select the rows of the Kronecker product matrix in which the row indices of each $B_x$ and $B_z$ element are the same. For example, each element of $B_x \otimes B_z$ is (I am going to use $x$ to mean an element of the $B_x$ design matrix and do the same for $z$ because I can't figure out how to do subscripts-on-subscripts)
$x_{ij}z_{rs}$
where
$i=[1:n]$ (rows of $B_x$)
$j=[1:q]$ (columns of $B_x$)
$r=[1:n]$ (rows of $B_z$)
$s=[1:p]$ (columns of $B_z$)
So you only select the rows where $i=r$. The result is an $n$-by-$qp$ matrix, which is what you want, because your tensor product basis has $n$ observations and a basis dimension of $qp$. I recommend looking at the full Kronecker product expansion in the above Wikipedia link, as it will make clear the order in which the row/column indices iterate as you move along the rows/columns of the matrix.
The rationale here is that each row of the design matrices $B_x$ and $B_z$ corresponds to the same actual observation, so the only values of the Kronecker product matrix that have any physical meaning are those in which the row values of each constituent matrix are the same.
