4
$\begingroup$

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.

Surface Plot

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – whuber Jan 4 '17 at 16:57
  • $\begingroup$ 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. $\endgroup$ – Eric Jan 4 '17 at 17:23
  • 1
    $\begingroup$ 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! $\endgroup$ – whuber Jan 4 '17 at 17:31
  • $\begingroup$ Good point. I will rewrite this to explain what is going on. $\endgroup$ – Eric Jan 4 '17 at 17:44
0
$\begingroup$

Not vectorizing your response matrix $Y$ is the way to go;

B = ginv(t(C) %*% C) %*% t(C) %*% Y  #OLS
pred = C%*%B  #Predictions
surface3d(x,x, pred, col = "green")  #Plot 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.