Ajuste de superficies mediante el producto tensorial de B-Splines

Question

Estoy tratando de enseñarme a mí mismo el ajuste de superficies con splines utilizando productos tensoriales. Estoy intentando construir un ejemplo de juguete pero no consigo que mi ejemplo funcione. Intentaré explicarlo lo mejor que pueda. Estoy tratando de ajustar un conjunto de puntos muestreados de Gaussianas bi-variadas usando B-Splines.

Este es el aspecto de la superficie muestreada.

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")

Quiero utilizar un sistema de base B-Spline de 6º orden con 12 funciones de base.

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)

(Nota: la función bsplineS() es del paquete fda)

Dado que estoy interpolando en una cuadrícula de 10 x 10. Las matrices de base B-Spline son las mismas. Por lo que entiendo, de la lectura de la literatura y esta pregunta Quiero hacer un producto Kronecker por filas de las dos matrices de base (ver la pregunta que he enlazado antes). Haciendo eso obtengo un $10 \times (12*12)$ matriz.

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

Aquí es donde estoy confundido, debería ser capaz de resolver el sistema utilizando las ecuaciones normales. Estoy haciendo un ajuste de mínimos cuadrados no ponderados por lo que debería ser capaz de resolver los coeficientes de regresión $\beta$ a través de las ecuaciones normales

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

En código r esto se convierte en

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

Sin embargo, sigo recibiendo el siguiente error:

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

Obviamente, estoy haciendo algo mal. Estoy teniendo problemas para resolver esto, no hay muchos ejemplos de código en la web. Cualquier ayuda, en particular un ejemplo de código sería muy apreciada.

Answer 1

No vectorizar su matriz de respuesta $Y$ es el camino a seguir;

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