# Surface Fit Using Tensor Product of B-Splines

I am trying to teach myself surface fitting with splines using tensor products. I am trying to construct a toy example but I can't seem to get my example to work. I will try to explain the best I can . I am trying to fit a set of a points sampled from bi-variate Gaussian using B-Splines.

This is what the sampled surface looks like.

library(rgl)
library(fda)

binorm = function(x,y)
{
z =  (15/ (2*pi)) * exp(-(1/2) * (x^2 + y^2))
return(z)
}

# Set up data
samples = 10
x = seq(-3, 3, length.out = samples)
grid = cbind(rep(x, times=1, each=samples), rep(x, times=samples)) # cartesian product
Y = matrix(binorm(grid[,1], grid[,2]), nrow = samples, ncol = samples) # evaluate at each point

# Plot the surface
surface3d(x,x, Y, col = "green")


I want to use a 6-th order B-Spline Basis system with 12 basis funcitons.

param.nbasis = 12
param.norder = 6

param.rangeval = c(min(x), max(x))
nbreaks = param.nbasis - param.norder + 2
basis.breaks = seq(param.rangeval[1], param.rangeval[2], length.out = nbreaks)

B.x = bsplineS(x, breaks = basis.breaks, norder = param.norder)
B.y = bsplineS(x, breaks = basis.breaks, norder = param.norder)


(Note: the bsplineS() function is from the fda package)

Given that I am interpolating on a 10 x 10 grid. The B-Spline basis matrices are the same. From what I understand, from reading the literature and this question I want to do a row-wise Kronecker product of the two Basis matrices (see the question I linked to before). Doing that I get a $10 \times (12*12)$ matrix.

C = matrix(nrow = 10, ncol = param.nbasis^2)
for (i in 1:10)
{
C[i,] = kronecker(B.x[i,], B.y[i,])
}


This is where I am confused, I should be able to solve the system using the normal equations. I am doing an unweighted least-squares fit so I should be able to solve for the regression coefficients $\beta$ via the normal equations

$$\beta = (C'C)^{-1} C'y$$

In r code this becomes

y = as.vector(Y)
B = solve(t(C) %*% C) %*% t(C) %*% y


However, I continue to get the following error:

Error in solve.default(t(C) %*% C) :
Lapack routine dgesv: system is exactly singular: U[7,7] = 0


Obviously I am doing something wrong. I am having trouble figuring this out, there are not many code examples on the web. Any help, particularly a code example would be greatly appreciated.

• This is purely a coding problem (not a statistical one), as a minimal amount of investigation will reveal. Make your parameters smaller so you can easily inspect your objects, then take a look at C: you will discover it has columns that are entirely zero. I do suspect there's a statistical question lurking here, though, concerning what the tensor product of splines is really trying to accomplish. – whuber Jan 4 '17 at 16:57
• Thank you, I was aware of there being 0 columns. To your point "there's a statistical question lurking here..." is exactly what I am trying to understand. I would say this is beyond "purely a coding problem". My lack of understanding of how the spline tensor product formalism is supposed to work is the issue, not the code itself. I merely included my code for completeness. – Eric Jan 4 '17 at 17:23
• You might make headway by explaining how your code is intended to work. I find it strange, for instance, that you compute both sets of univariate splines from the $x$ coordinates alone and even stranger that you are creating $144$ variables to describe just $100$ responses! – whuber Jan 4 '17 at 17:31
• Good point. I will rewrite this to explain what is going on. – Eric Jan 4 '17 at 17:44

Not vectorizing your response matrix $Y$ is the way to go;
B = ginv(t(C) %*% C) %*% t(C) %*% Y  #OLS