# Smoothing splines as basis expansion

I am following the discussion on chapter 5 of Elements of statistical learning which discusses basis expansion using splines. The data set I used is the Ozone data which can be found at :http://web.stanford.edu/~hastie/ElemStatLearn/ It has a response variable "ozone",and three predictor variables "radiation","temperature" and "wind".
Essentially what I am trying to do is get the basis matrix for each one of the predictor variables and combine them to train a generalized linear additive model to predict the ozone levels. To be more specific the underlying model is : $$Y = \theta_{0} + \boldsymbol{h}_{1}^{T}\boldsymbol{\theta}_{1} + \boldsymbol{h}_{2}^{T}\boldsymbol{\theta}_{2}+ \boldsymbol{h}_{3}^{T}\boldsymbol{\theta}_{3} = \theta_{0} + f_{1}(X_{1}) +f_{2}(X_{2})+f_{3}(X_{3})$$ , where Where $$\theta_{0}$$ is the intercept term. $$X_{i}$$ is the $$i^{th}$$ feature of data X. And $$\boldsymbol{h}_{1}(X_{1})$$ is of form $$\begin{bmatrix} h_{11}(X_{1})\\ h_{12}(X_{1})\\ h_{13}(X_{1})\\ h_{14}(X_{1})\\ \vdots \end{bmatrix} ,$$. Each entry of $$\boldsymbol{h}_{1}(X_{1})$$ is a spline basis function. That is each variable is expanded by the spline basis (natural cubic spline or smoothing spline etc.) and then do a linear regression on the whole model. The smoothing spline is found by minimizing $$RSS(f,\lambda) = \sum_{i}^{N}\{y_{i}-f(x_{i})\}^{2} + \lambda\int\{f''(t)\}dt$$, and the solution is given by $$\hat{\theta} = (N^{T}N + \lambda\Omega_{N})^{-1}N^{T}\boldsymbol{y}$$, where $$\{\Omega_{N}\}_{jk} = \int N''_{j}(t)N''_{k}(t) dt$$
I am really stuck on how to calculate the $$\Omega$$ matrix here. Currently I am using R and with the built-in function "ns(x,knots = )" from which I can get a basis matrix with specified knot locations. Since smoothing splines place knots at each distinctive x value , the command I use is ns(Ozone$radiation,knots = sort(Ozone$radiation)[2:110]). The data matrix is of size (111,3), since ns function does not allow placing knots at boundary points I only select the values from index 2 to index 110. Not surprisingly the resulting basis matrix is very sparse thus it is not viable to take the inverse and need to add the $$\Omega$$ matrix in order to get the $$\hat{\theta}$$.But I am not sure how to calculate that , what is the range of integration? And is there libraries that I can use to calculate the derivatives of the basis functions?

Aside
I have tried a simpler model with natural cubic splines with degree of freedom = 4. And below is the code I used. The figure plots "radiation v.s. f(radiation)" where f(radiation) is the $$f_{1}(X_{1})$$ in my model and it is equal to $$\boldsymbol{h}_{1}^{T}\boldsymbol{\theta}_{1}$$ and the red lines are the confidence band obtained following the procedures in Section 5.2.2 in ESL.

require(stats); require(graphics)
library(splines)
library(ggplot2)

train = OZ[,c(2,3,4)]

yTrain=OZ$ozone yTrain = scale(yTrain,center = TRUE, scale = TRUE) X=cbind(ns(train$$radiation,4),ns(train$$temperature,4),ns(train$wind,4))
X = scale(X,center = TRUE, scale = TRUE)
X
df_2=data.frame(yTrain,X)

regression = lm(yTrain~.,data = df_2)
B = cbind(1,X)
N = length(yTrain)
p = 4*3+1
#########################
#    Yhat, theta_hat    #
yhat = B%*%(solve(t(B)%*%B))%*%(t(B))%*%(yTrain)
theta_hat = (solve(t(B)%*%B))%*%(t(B))%*%(yTrain)
#########################
#        J              #
J = (solve(t(B)%*%B))%*%(t(B))

#########################

sigma_val = (((t(yTrain - yhat))%*%(yTrain - yhat))/(N-p-1))[1]
V_y = sigma_val * diag(N)
#########################
# Sigma_hat             #
Sigma_hat = (J%*%V_y)%*%(t(J))
#########################

#########################