Understand smoothing splines

1. Compositional spline

This code implements the compositional smoothing splines grounded on the theory of Bayes spaces. The below code is based on the function compositionalSpline() in the robCompositions package.

The code is based on the method in Machalova, J., Talska, R., Hron, K. Gaba, A. Compositional splines for representation of density functions. Comput Stat (2020). https://doi.org/10.1007/s00180-020-01042-7

The target is to estimate the fomular (19) in their formular

\[ \begin{aligned} J_{l}(\mathbf{z})=&(1-\alpha) \mathbf{z}^{\top} \mathbf{K}^{\top} \mathbf{D S}_{l}^{\top} \mathbf{M}_{k l} \mathbf{S}_{l} \mathbf{D} \mathbf{K} \mathbf{z}+\\ &+\alpha\left[\mathbf{y}-\mathbf{B}_{k+1}(\mathbf{x}) \mathbf{D} \mathbf{K}\right]^{\top} \mathbf{W}\left[\mathbf{y}-\mathbf{B}_{k+1}(\mathbf{x}) \mathbf{D K} \mathbf{z}\right] \end{aligned} \]

2. Simulation data

Simulate 10 design points: Here, denoted as knots. The notation in the paper is \(\lambda\) and then \(g = 8\).

Simulate a grid of 100 values based on normal density, here, denoted as t. The notation in the paper is \(x\) and then \(n = 100\)

#--their function
require(robCompositions)
require(splines)
# Example (normal density)
t =  seq(-4.7,4.7, length = 100)
t_step = diff(t[1:2])
mean = 0; sd = 1.5
f = dnorm(t, mean, sd)
f1 = f/trapzc(t_step,f)

Using fcenLR() function with above values

f.fcenLR = fcenLR(t, t_step,  f) 
length(f.fcenLR)

## [1] 100

knots <- seq(-4.7,4.7, length = 10)
knots 

##  [1] -4.7000000 -3.6555556 -2.6111111 -1.5666667 -0.5222222  0.5222222
##  [7]  1.5666667  2.6111111  3.6555556  4.7000000

3. Test each step in compositionalSpline()

Set \(order = 4 = k+1\).

Derivative \(l = der = 2\).

# der : for the penlized

order = 4
der = 2
k = order
r = length(knots)
r

## [1] 10

w = rep( 1/length(t), length(t))
# see w
length(w)

## [1] 100

clrf = f.fcenLR
length(clrf)

## [1] 100

3.1 Create knots

The below code create knots, as below. Then, the knots \(\lambda_{1}\), .., \(\lambda_{r}\) are knots inside the intervals and there are \(2*k\) extra knots. \[ \begin{aligned} &\Delta \lambda:=\lambda_{0}=u<\lambda_{1}<\cdots<\lambda_{10}<v=\dot{\lambda}_{r+1} \\ &\lambda_{-k}=\cdots=\lambda_{-1}=\lambda_{0}, \quad \lambda_{r+1}=\lambda_{r+2}=\cdots=\lambda_{r+k+1} \end{aligned} \]

Count both inside knots, 2 borders and extra knots

Celkova_Delka = 2 * (k - 1) + r  # 
Celkova_Delka

## [1] 16

lambda = c()
for (i in 1:(Celkova_Delka)) {
  if (i <= k - 1) {
    lambda[i] = knots[1]
  }
  if ((i > k - 1) && (i <= r + k - 1)) {
    lambda[i] = knots[i - (k - 1)]
  }
  if (i > r + k - 1) {
    lambda[i] = knots[r]
  }
}
lambda #### 16 knots = 8 insides + 2 borders + 6 extra knots #  6 = 2*(order -1)

##  [1] -4.7000000 -4.7000000 -4.7000000 -4.7000000 -3.6555556 -2.6111111
##  [7] -1.5666667 -0.5222222  0.5222222  1.5666667  2.6111111  3.6555556
## [13]  4.7000000  4.7000000  4.7000000  4.7000000

Then it means

\[ \begin{aligned} &\Delta \lambda:=\lambda_{0}=-4.7=u<\lambda_{1}=-3.66<\cdots<\lambda_{8}=3.66<v=4.7=\dot{\lambda}_{9} \\ &\lambda_{-3}=\lambda_{-2}=\lambda_{-1}=\lambda_{0}=-4.7, \quad \lambda_{9}=\lambda_{10}=\lambda_{11}=\lambda_{12} \end{aligned} \]

3.2. Create a B-splines basic at order \(4\), based on knots \(\lambad\)

splineDesign () is a function to Design Matrix for B-splines basic, given order \(4\), a total of 16 knots (inside, border and extra) and the value t

ord =  k
B = splineDesign(lambda, t, ord, outer.ok = TRUE)
dim(B) 

## [1] 100  12

length(t)

## [1] 100

length(lambda)- ord

## [1] 12

See matrix B

head(B, 12)

##               [,1]      [,2]       [,3]         [,4] [,5] [,6] [,7] [,8] [,9]
##  [1,] 1.0000000000 0.0000000 0.00000000 0.0000000000    0    0    0    0    0
##  [2,] 0.7513148009 0.2368520 0.01170799 0.0001252191    0    0    0    0    0
##  [3,] 0.5477084899 0.4072126 0.04407713 0.0010017531    0    0    0    0    0
##  [4,] 0.3846731781 0.5189707 0.09297521 0.0033809166    0    0    0    0    0
##  [5,] 0.2577009767 0.5800150 0.15426997 0.0080140245    0    0    0    0    0
##  [6,] 0.1622839970 0.5982344 0.22382920 0.0156523917    0    0    0    0    0
##  [7,] 0.0939143501 0.5815177 0.29752066 0.0270473328    0    0    0    0    0
##  [8,] 0.0480841473 0.5377536 0.37121212 0.0429501628    0    0    0    0    0
##  [9,] 0.0202854996 0.4748310 0.44077135 0.0641121963    0    0    0    0    0
## [10,] 0.0060105184 0.4006386 0.50206612 0.0912847483    0    0    0    0    0
## [11,] 0.0007513148 0.3230654 0.55096419 0.1252191335    0    0    0    0    0
## [12,] 0.0000000000 0.2500000 0.58333333 0.1666666667    0    0    0    0    0
##       [,10] [,11] [,12]
##  [1,]     0     0     0
##  [2,]     0     0     0
##  [3,]     0     0     0
##  [4,]     0     0     0
##  [5,]     0     0     0
##  [6,]     0     0     0
##  [7,]     0     0     0
##  [8,]     0     0     0
##  [9,]     0     0     0
## [10,]     0     0     0
## [11,]     0     0     0
## [12,]     0     0     0

tail(B, 12)

##        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]         [,9]      [,10]
##  [89,]    0    0    0    0    0    0    0    0 0.1666666667 0.58333333
##  [90,]    0    0    0    0    0    0    0    0 0.1252191335 0.55096419
##  [91,]    0    0    0    0    0    0    0    0 0.0912847483 0.50206612
##  [92,]    0    0    0    0    0    0    0    0 0.0641121963 0.44077135
##  [93,]    0    0    0    0    0    0    0    0 0.0429501628 0.37121212
##  [94,]    0    0    0    0    0    0    0    0 0.0270473328 0.29752066
##  [95,]    0    0    0    0    0    0    0    0 0.0156523917 0.22382920
##  [96,]    0    0    0    0    0    0    0    0 0.0080140245 0.15426997
##  [97,]    0    0    0    0    0    0    0    0 0.0033809166 0.09297521
##  [98,]    0    0    0    0    0    0    0    0 0.0010017531 0.04407713
##  [99,]    0    0    0    0    0    0    0    0 0.0001252191 0.01170799
## [100,]    0    0    0    0    0    0    0    0 0.0000000000 0.00000000
##            [,11]        [,12]
##  [89,] 0.2500000 0.0000000000
##  [90,] 0.3230654 0.0007513148
##  [91,] 0.4006386 0.0060105184
##  [92,] 0.4748310 0.0202854996
##  [93,] 0.5377536 0.0480841473
##  [94,] 0.5815177 0.0939143501
##  [95,] 0.5982344 0.1622839970
##  [96,] 0.5800150 0.2577009767
##  [97,] 0.5189707 0.3846731781
##  [98,] 0.4072126 0.5477084899
##  [99,] 0.2368520 0.7513148009
## [100,] 0.0000000 1.0000000000

Plot the matrix B, then For the first 4 column, i.e the first 4 knots, it has values at 11 rows. It corresponds the back curve (up to 1). Then, the matrix 4x11 is moved until at the end of matrix B. Each block corresponds to each colorful curves. The last bock correspond to purple curve. (up to 1)

plot(range(t), c(0,1), type = "n", xlab = "x", ylab = "",
     main =  "B-splines at order = 4 - sum to 1 inside inner knots")
#mtext(expression(B[j](x) *"  and "* sum(B[j](x), j == 1, 6)), adj = 0)
abline(v = lambda, lty = 3, col = "light gray")
abline(v = lambda[c(4,length(lambda)-3)], lty = 3, col = "gray10")
#lines(x, rowSums(bb), col = "gray", lwd = 2)
matlines(t, B, ylim = c(0,1), lty = 1)

3.3 Matrix W

  W = diag(w)
dim(W)

## [1] 100 100

3.4. Create all B-splines up to order \(k+1 = 4\)

BB is a matrix to store all B-splines at knots l, evaluated at vector partition, order k.i.e \[ \boldsymbol{B}_{k+1} \]

  parnition = seq(min(lambda), max(lambda), length = 100)
head(  parnition, 20 )

##  [1] -4.700000 -4.605051 -4.510101 -4.415152 -4.320202 -4.225253 -4.130303
##  [8] -4.035354 -3.940404 -3.845455 -3.750505 -3.655556 -3.560606 -3.465657
## [15] -3.370707 -3.275758 -3.180808 -3.085859 -2.990909 -2.895960

  lambda_index = c(0:(r - 1))
  lambda_index

##  [1] 0 1 2 3 4 5 6 7 8 9

The lambda_index indicates the below knots

\[ \lambda_{0}=-4.7=u<\lambda_{1}=-3.66<\cdots<\lambda_{8}=3.66<v=4.7=\dot{\lambda}_{9} \]

g = lambda_index[length(lambda_index) - 1]
g

## [1] 8

k

## [1] 4

  N = g + (k - 1) + 1
  N

## [1] 12

Then \(g\) indicates the last inside knot, here, \(\lambda_8\).

Similar as above, \(N\) indicates the number of B-splines basic at order \(k=4\) given 8 insides knots. We then have a total of 16 knots (inside knots, border knots and extra knots). We have

\[Celkova_Delka = 2 * (k - 1) + r = length(lambda) \] The number of B-spline basis from splineDesign() is \[length(lambda)- ord = length(lambda)- k = 2 * (k - 1) + r - k = k+r -2 \] On the other sides, given \(r\) knots, there are \(r-2\) inside knots, i.e \(g = r-2\). Then, \[ N = g + (k - 1) + 1 = r - 2 +(k - 1) + 1 =r+k-2 \] Below, they create a B-splines, at order \(k\) and given values \(t = parnition = seq(min(lambda), max(lambda), length = 100)\). Check value \(l\) for the next chunk

  l = c()
  for (i in (1:N)) {
    for (j in 1:(k + 1)) {
      l[j] = lambda[i + j - 1]
      print(c(i , j, round(l, 1)))
    }
  }

## [1]  1.0  1.0 -4.7
## [1]  1.0  2.0 -4.7 -4.7
## [1]  1.0  3.0 -4.7 -4.7 -4.7
## [1]  1.0  4.0 -4.7 -4.7 -4.7 -4.7
## [1]  1.0  5.0 -4.7 -4.7 -4.7 -4.7 -3.7
## [1]  2.0  1.0 -4.7 -4.7 -4.7 -4.7 -3.7
## [1]  2.0  2.0 -4.7 -4.7 -4.7 -4.7 -3.7
## [1]  2.0  3.0 -4.7 -4.7 -4.7 -4.7 -3.7
## [1]  2.0  4.0 -4.7 -4.7 -4.7 -3.7 -3.7
## [1]  2.0  5.0 -4.7 -4.7 -4.7 -3.7 -2.6
## [1]  3.0  1.0 -4.7 -4.7 -4.7 -3.7 -2.6
## [1]  3.0  2.0 -4.7 -4.7 -4.7 -3.7 -2.6
## [1]  3.0  3.0 -4.7 -4.7 -3.7 -3.7 -2.6
## [1]  3.0  4.0 -4.7 -4.7 -3.7 -2.6 -2.6
## [1]  3.0  5.0 -4.7 -4.7 -3.7 -2.6 -1.6
## [1]  4.0  1.0 -4.7 -4.7 -3.7 -2.6 -1.6
## [1]  4.0  2.0 -4.7 -3.7 -3.7 -2.6 -1.6
## [1]  4.0  3.0 -4.7 -3.7 -2.6 -2.6 -1.6
## [1]  4.0  4.0 -4.7 -3.7 -2.6 -1.6 -1.6
## [1]  4.0  5.0 -4.7 -3.7 -2.6 -1.6 -0.5
## [1]  5.0  1.0 -3.7 -3.7 -2.6 -1.6 -0.5
## [1]  5.0  2.0 -3.7 -2.6 -2.6 -1.6 -0.5
## [1]  5.0  3.0 -3.7 -2.6 -1.6 -1.6 -0.5
## [1]  5.0  4.0 -3.7 -2.6 -1.6 -0.5 -0.5
## [1]  5.0  5.0 -3.7 -2.6 -1.6 -0.5  0.5
## [1]  6.0  1.0 -2.6 -2.6 -1.6 -0.5  0.5
## [1]  6.0  2.0 -2.6 -1.6 -1.6 -0.5  0.5
## [1]  6.0  3.0 -2.6 -1.6 -0.5 -0.5  0.5
## [1]  6.0  4.0 -2.6 -1.6 -0.5  0.5  0.5
## [1]  6.0  5.0 -2.6 -1.6 -0.5  0.5  1.6
## [1]  7.0  1.0 -1.6 -1.6 -0.5  0.5  1.6
## [1]  7.0  2.0 -1.6 -0.5 -0.5  0.5  1.6
## [1]  7.0  3.0 -1.6 -0.5  0.5  0.5  1.6
## [1]  7.0  4.0 -1.6 -0.5  0.5  1.6  1.6
## [1]  7.0  5.0 -1.6 -0.5  0.5  1.6  2.6
## [1]  8.0  1.0 -0.5 -0.5  0.5  1.6  2.6
## [1]  8.0  2.0 -0.5  0.5  0.5  1.6  2.6
## [1]  8.0  3.0 -0.5  0.5  1.6  1.6  2.6
## [1]  8.0  4.0 -0.5  0.5  1.6  2.6  2.6
## [1]  8.0  5.0 -0.5  0.5  1.6  2.6  3.7
## [1] 9.0 1.0 0.5 0.5 1.6 2.6 3.7
## [1] 9.0 2.0 0.5 1.6 1.6 2.6 3.7
## [1] 9.0 3.0 0.5 1.6 2.6 2.6 3.7
## [1] 9.0 4.0 0.5 1.6 2.6 3.7 3.7
## [1] 9.0 5.0 0.5 1.6 2.6 3.7 4.7
## [1] 10.0  1.0  1.6  1.6  2.6  3.7  4.7
## [1] 10.0  2.0  1.6  2.6  2.6  3.7  4.7
## [1] 10.0  3.0  1.6  2.6  3.7  3.7  4.7
## [1] 10.0  4.0  1.6  2.6  3.7  4.7  4.7
## [1] 10.0  5.0  1.6  2.6  3.7  4.7  4.7
## [1] 11.0  1.0  2.6  2.6  3.7  4.7  4.7
## [1] 11.0  2.0  2.6  3.7  3.7  4.7  4.7
## [1] 11.0  3.0  2.6  3.7  4.7  4.7  4.7
## [1] 11.0  4.0  2.6  3.7  4.7  4.7  4.7
## [1] 11.0  5.0  2.6  3.7  4.7  4.7  4.7
## [1] 12.0  1.0  3.7  3.7  4.7  4.7  4.7
## [1] 12.0  2.0  3.7  4.7  4.7  4.7  4.7
## [1] 12.0  3.0  3.7  4.7  4.7  4.7  4.7
## [1] 12.0  4.0  3.7  4.7  4.7  4.7  4.7
## [1] 12.0  5.0  3.7  4.7  4.7  4.7  4.7

From the above results, given i and j, the vector \(l\) has 1, 2, 3 or 4 knots (the most common). Then, they use the vector \(l\) knots to create B-splines basis, at order 4 and value parnition. Each B-spline is stored in matrix BB.

Question: Where do we apply matrix BB????? Answer: BB is use to plot the ZB-spline basis system, see again at the end of this file.

N 

## [1] 12

BB = array(0, c(length(parnition), N))
  l = c()
  for (i in (1:N)) {
    for (j in 1:(k + 1)) {
      l[j] = lambda[i + j - 1]
    }
    BB[, i] = splineDesign(l, parnition, k, outer.ok = TRUE)
  }
dim(BB)

## [1] 100  12

see again BB, seems to be simlar to B!!!, plot below

head(BB, 12)

##               [,1]      [,2]       [,3]         [,4] [,5] [,6] [,7] [,8] [,9]
##  [1,] 1.0000000000 0.0000000 0.00000000 0.0000000000    0    0    0    0    0
##  [2,] 0.7513148009 0.2368520 0.01170799 0.0001252191    0    0    0    0    0
##  [3,] 0.5477084899 0.4072126 0.04407713 0.0010017531    0    0    0    0    0
##  [4,] 0.3846731781 0.5189707 0.09297521 0.0033809166    0    0    0    0    0
##  [5,] 0.2577009767 0.5800150 0.15426997 0.0080140245    0    0    0    0    0
##  [6,] 0.1622839970 0.5982344 0.22382920 0.0156523917    0    0    0    0    0
##  [7,] 0.0939143501 0.5815177 0.29752066 0.0270473328    0    0    0    0    0
##  [8,] 0.0480841473 0.5377536 0.37121212 0.0429501628    0    0    0    0    0
##  [9,] 0.0202854996 0.4748310 0.44077135 0.0641121963    0    0    0    0    0
## [10,] 0.0060105184 0.4006386 0.50206612 0.0912847483    0    0    0    0    0
## [11,] 0.0007513148 0.3230654 0.55096419 0.1252191335    0    0    0    0    0
## [12,] 0.0000000000 0.2500000 0.58333333 0.1666666667    0    0    0    0    0
##       [,10] [,11] [,12]
##  [1,]     0     0     0
##  [2,]     0     0     0
##  [3,]     0     0     0
##  [4,]     0     0     0
##  [5,]     0     0     0
##  [6,]     0     0     0
##  [7,]     0     0     0
##  [8,]     0     0     0
##  [9,]     0     0     0
## [10,]     0     0     0
## [11,]     0     0     0
## [12,]     0     0     0

3.5. Check the order of derivative, Collocaton matrix.

Recall collocaton matrix \[ \begin{aligned} &\mathbf{B}_{k+1}(\mathbf{x})=\left(B_{i}^{k+1}\left(x_{j}\right)\right)_{j=1, i=-k}^{n, g} \\ &\mathbf{C}_{k+1}(\mathbf{x})=\left(\begin{array}{ccc} B_{-k}^{k+1}\left(x_{1}\right) & \cdots & B_{g}^{k+1}\left(x_{1}\right) \\ \vdots & \ddots & \vdots \\ B_{-k}^{k+1}\left(x_{n}\right) & \cdots & B_{g}^{k+1}\left(x_{n}\right) \end{array}\right) \in \mathbb{R}^{n, g+k+1} \end{aligned} \]

length(t) <= N

## [1] FALSE

qr(B)$rank != N

## [1] FALSE

  if (length(t) <= N) 
    stop("length(t) must be higher then Dimension(space of splines)")
  if (qr(B)$rank != N) 
    stop("Collocaton matrix does not have full column rank.")

3.6. Create matrix the Upper triangular matrix \(S\)

Matrix \(S_l\), with \(l\) is the \(l\)th derivative in the below equation

\[J_{l}\left(s_{k}\right)=(1-\alpha) \int_{a}^{b}\left[s_{k}^{(l)}(x)\right]^{2} \mathrm{~d} x+\alpha \sum_{i=1}^{n} w_{i}\left[y_{i}-s_{k}\left(x_{i}\right)\right]^{2}\]

3.6.1 The general form of \(S_l\)

\[ \mathbf{S}_{l}=\mathbf{D}_{l} \mathbf{L}_{l} \ldots \mathbf{D}_{1} \mathbf{L}_{1} \in \mathbb{R}^{g+k+1-l, g+k+1} \]

where \[ \mathbf{D}_{j}=(k+1-j) \operatorname{diag}\left(d_{-k+j}, \ldots, d_{g}\right) \] with \[ d_{i}=\frac{1}{\lambda_{i+k+1-j}-\lambda_{i}}, \quad i=-k+j, \ldots, g \] and \[ \mathbf{L}_{j}:=\left(\begin{array}{rrrr} -1 & 1 & & \\ & \ddots & \ddots & \\ & & -1 & 1 \end{array}\right) \in \mathbb{R}^{g+k+1-j, g+k+2-j} \]

3.6.2 The form of \(S_l\) when \(l=2\).

\[ \mathbf{S}_{2}=\mathbf{D}_{2} \mathbf{L}_{2} \mathbf{D}_{1} \mathbf{L}_{1} \in \mathbb{R}^{g+k+1-2, g+k+1} = \mathbb{R}^{g+k-1, g+k+1} \]

where

\[ \begin{align} &\mathbf{D}_{1}=(k+1-1) \operatorname{diag}\left(d_{-k+1}, \ldots, d_{g}\right) = k \operatorname{diag}\left(d_{-k+1}, \ldots, d_{g}\right),\quad d_{i}=\frac{1}{\lambda_{i+k+1-1}-\lambda_{i}}= \frac{1}{\lambda_{i+k}-\lambda_{i}}, \quad i=-k+1, \ldots, g\\ & \mathbf{D}_{2}=(k+1-2) \operatorname{diag}\left(d_{-k+2}, \ldots, d_{g}\right) = (k-1) \operatorname{diag}\left(d_{-k+2}, \ldots, d_{g}\right), \quad d_{i}=\frac{1}{\lambda_{i+k+1-2}-\lambda_{i}}= \frac{1}{\lambda_{i+k-1}-\lambda_{i}}, \quad i=-k+2, \ldots, g \end{align} \]

3.6.3 Check the code with \(l = der = 2\)

#S = array(0)
S_pom = diag(1, N, N)
    
for (j in 1:der) {
      D_mat = array(0)
      rozdil = lambda[(1 + k):(N + k - j)] - lambda[(1 +   j):(N)]
      D_mat = (k - j) * diag(1/rozdil)
      L_mat = array(0, c(N - j, N - j + 1))
          for (J in (1:(N - j))) {
              L_mat[J, J] = (-1)
              L_mat[J, J + 1] = 1
            }
      S_pom = D_mat %*% L_mat %*% S_pom
    }
    S = S_pom
    dim(S)

## [1] 10 12

The following code, we see the matrix D and matrix L, at value \(j = 1\)

S_pom = diag(1, N, N)
j = 1
D_mat = array(0)
rozdil = lambda[(1 + k):(N + k - j)] - lambda[(1 +   j):(N)]
D_mat = (k - j) * diag(1/rozdil)
L_mat = array(0, c(N - j, N - j + 1))
          for (J in (1:(N - j))) {
              L_mat[J, J] = (-1)
              L_mat[J, J + 1] = 1
          }
#print matrix
head(D_mat)

##         [,1]    [,2]      [,3]      [,4]      [,5]      [,6] [,7] [,8] [,9]
## [1,] 2.87234 0.00000 0.0000000 0.0000000 0.0000000 0.0000000    0    0    0
## [2,] 0.00000 1.43617 0.0000000 0.0000000 0.0000000 0.0000000    0    0    0
## [3,] 0.00000 0.00000 0.9574468 0.0000000 0.0000000 0.0000000    0    0    0
## [4,] 0.00000 0.00000 0.0000000 0.9574468 0.0000000 0.0000000    0    0    0
## [5,] 0.00000 0.00000 0.0000000 0.0000000 0.9574468 0.0000000    0    0    0
## [6,] 0.00000 0.00000 0.0000000 0.0000000 0.0000000 0.9574468    0    0    0
##      [,10] [,11]
## [1,]     0     0
## [2,]     0     0
## [3,]     0     0
## [4,]     0     0
## [5,]     0     0
## [6,]     0     0

head(L_mat)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]   -1    1    0    0    0    0    0    0    0     0     0     0
## [2,]    0   -1    1    0    0    0    0    0    0     0     0     0
## [3,]    0    0   -1    1    0    0    0    0    0     0     0     0
## [4,]    0    0    0   -1    1    0    0    0    0     0     0     0
## [5,]    0    0    0    0   -1    1    0    0    0     0     0     0
## [6,]    0    0    0    0    0   -1    1    0    0     0     0     0

The following code, we see the matrix D and matrix L, at value \(j = 2\)

S_pom = diag(1, N, N)
j = 2
D_mat = array(0)
rozdil = lambda[(1 + k):(N + k - j)] - lambda[(1 +   j):(N)]
D_mat = (k - j) * diag(1/rozdil)
L_mat = array(0, c(N - j, N - j + 1))
          for (J in (1:(N - j))) {
              L_mat[J, J] = (-1)
              L_mat[J, J + 1] = 1
          }
#print matrix
head(D_mat)

##          [,1]      [,2]      [,3]      [,4]      [,5]      [,6] [,7] [,8] [,9]
## [1,] 1.914894 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000    0    0    0
## [2,] 0.000000 0.9574468 0.0000000 0.0000000 0.0000000 0.0000000    0    0    0
## [3,] 0.000000 0.0000000 0.9574468 0.0000000 0.0000000 0.0000000    0    0    0
## [4,] 0.000000 0.0000000 0.0000000 0.9574468 0.0000000 0.0000000    0    0    0
## [5,] 0.000000 0.0000000 0.0000000 0.0000000 0.9574468 0.0000000    0    0    0
## [6,] 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.9574468    0    0    0
##      [,10]
## [1,]     0
## [2,]     0
## [3,]     0
## [4,]     0
## [5,]     0
## [6,]     0

head(L_mat)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,]   -1    1    0    0    0    0    0    0    0     0     0
## [2,]    0   -1    1    0    0    0    0    0    0     0     0
## [3,]    0    0   -1    1    0    0    0    0    0     0     0
## [4,]    0    0    0   -1    1    0    0    0    0     0     0
## [5,]    0    0    0    0   -1    1    0    0    0     0     0
## [6,]    0    0    0    0    0   -1    1    0    0     0     0

3.7 Prepare matrix $M_{kl} = M_{3,2} $

Recall matrix \(M_{kl}\) is \[ \mathbf{M}_{k l}=\left(m_{i j}^{k l}\right)_{i, j=-k+l}^{g}, \quad \text { with } \quad m_{i j}^{k l}=\int_{a}^{b} B_{i}^{k+1-l}(x) B_{j}^{k+1-l}(x) \mathrm{d} x \] In the below code, kk means \(kk = k+1 - l\). However, as Christine reminding, we usually work with \(l = 2.\)

Then, given order kk, we need a total of celkova_delka knots (inside knots, border knots and extra knots )

kk = k - der
kk

## [1] 2

  celkova_delka = 2 * (kk - 1) + r
   celkova_delka 

## [1] 12

The below code create a total knots for B-spline order \(kk\)!!! They create extra knots based on the given Lambda knots (including inside knots and border knots).

  Lambda = c()
  for (i in 1:celkova_delka) {
    if (i <= (kk - 1)) {
      Lambda[i] = knots[1]
    }
    if ((i > kk - 1) && (i <= r + kk - 1)) {
      Lambda[i] = knots[i - (kk - 1)]
    }
    if (i > (r + (kk - 1))) {
      Lambda[i] = knots[r]
    }
  }
Lambda

##  [1] -4.7000000 -4.7000000 -3.6555556 -2.6111111 -1.5666667 -0.5222222
##  [7]  0.5222222  1.5666667  2.6111111  3.6555556  4.7000000  4.7000000

correspond to

\[ \begin{aligned} &\Delta \lambda:=\lambda_{0}=-4.7=u<\lambda_{1}=-3.66<\cdots<\lambda_{8}=3.66<v=4.7=\dot{\lambda}_{9} \\ &\lambda_{-1}=\lambda_{0}=-4.7, \quad \lambda_{9} = \lambda_{10}=4.7 \end{aligned} \] The following code create a B-spline at order kk based on Lambda knots.

 Parnition = seq(min(Lambda), max(Lambda), length = 100)
 BBB = splineDesign(Lambda, Parnition, kk, outer.ok = TRUE)

 dim(BBB)

## [1] 100  10

plot the BBB basic, it is correct since they are at order \(kk=1.\)

plot(range(Parnition), c(0,1), type = "n", xlab = "x", ylab = "",
     main =  "B-splines at order kk = 1, using knots Lambda")
#mtext(expression(B[j](x) *"  and "* sum(B[j](x), j == 1, 6)), adj = 0)
abline(v =  Lambda, lty = 3, col = "light gray")
abline(v =  Lambda[c(4,length(Lambda)-3)], lty = 3, col = "gray10")
#lines(x, rowSums(bb), col = "gray", lwd = 2)
matlines(Parnition, BBB, ylim = c(0,1), lty = 1)

The following code finishes the calculation of M

\[ m_{i j}^{k l}=\int_{a}^{b} B_{i}^{k+1-l}(x) B_{j}^{k+1-l}(x) \mathrm{d} x \]

     Lambda_index = c(0:(r - 1)) 
  G = Lambda_index[length(Lambda_index) - 1]
  NN = G + (kk - 1) + 1
  step = diff(Parnition[1:2])
  M = array(0, c(NN, NN))
  for (i in 1:NN) {
    for (j in 1:NN) {
      nenulove = c()
      soucin = BBB[, i] * BBB[, j]
      for (m in 1:length(Parnition)) {
        if (soucin[m] != 0) {
          nenulove[m] = soucin[m]
        }
      }
      M[i, j] = trapzc(step, soucin)
    }
  }
 dim(M)

## [1] 10 10

4. Find the optimal function

4.1 Create matrix D, K, U and G

\[ \mathbf{D}_{j}=(k+1-j) \operatorname{diag}\left(d_{-k+j}, \ldots, d_{g}\right) \] with \[ d_{i}=\frac{1}{\lambda_{i+k+1-j}-\lambda_{i}}, \quad i=-k+j, \ldots, g \]

difference = lambda[(1 + k):(r + 2 * (k - 1))] - lambda[(1:(r +  k - 2))]
D = (k) * diag(1/difference)
K = array(0, c(N, N - 1))
K[1, 1] = 1
K[N, N - 1] = -1
for (j in (2:(N - 1))) {
  K[j, j - 1] = (-1)
  K[j, j] = 1
}
# See matrix D and matrix K
dim(D)

## [1] 12 12

dim(K)

## [1] 12 11

Then, Using the notation \[\mathbf{U}:=\mathbf{D K}\]

U = D %*% K
dim(U)

## [1] 12 11

Then, calculate matrix G \[ \mathbf{G}:=\mathbf{U}^{\top}\left[(1-\alpha) \mathbf{S}_{l}^{\top} \mathbf{M}_{k l} \mathbf{S}_{l}+\alpha \mathbf{B}_{k+1}^{\top}(\mathbf{x}) \mathbf{W B}_{k+1}(\mathbf{x})\right] \mathbf{U} \]

alpha = 0.5891077 #1.950193
GG = t(U) %*% ((1 - alpha) * t(S) %*% M %*% S + alpha * t(B) %*%  W %*% B) %*% U
dim(GG)

## [1] 11 11

Then, calculate matrix \(g\) as below

\[ \mathbf{g}:=\alpha \mathbf{K}^{\top} \mathbf{D} \mathbf{B}_{k+1}^{\top}(\mathbf{x}) \mathbf{W} \mathbf{y} \]

gg = alpha * t(K) %*% t(D) %*% t(B) %*% W %*% clrf
dim(gg)

## [1] 11  1

See clrf, i,e y

length(clrf)

## [1] 100

dim(alpha * t(K) %*% t(D) %*% t(B) %*% W)

## [1]  11 100

length(clrf)

## [1] 100

4.2 Solve the optimal

The optimal solution

z = solve(GG) %*% gg
length(z)

## [1] 11

z

##                [,1]
##  [1,] -5.495245e-01
##  [2,] -1.456176e+00
##  [3,] -2.235599e+00
##  [4,] -2.196856e+00
##  [5,] -1.325074e+00
##  [6,]  1.245115e-13
##  [7,]  1.325074e+00
##  [8,]  2.196856e+00
##  [9,]  2.235599e+00
## [10,]  1.456176e+00
## [11,]  5.495245e-01

The optimal ZB-spline, i.e Zbasis and it is the BDK in the following equations:

\[ s_{k}^{*}(x)=\mathbf{B}_{k+1}(x) \mathbf{D K z}^{*} \]

Zbasis = B %*% D %*% K
dim(Zbasis)

## [1] 100  11

Zbasis_finner = BB %*% D %*% K
dim(Zbasis_finner)

## [1] 100  11

5. Plot the final results

5.1. Plot the ZB-spline basis system

#if (basis.plot == TRUE) {
  matplot(parnition, Zbasis_finner, type = "l", lty = 1, 
          las = 1, xlab = "t", ylab = "fcenLR(density)", 
          col = rainbow(dim(Zbasis_finner)[2]), main = "ZB-spline basis")
  abline(v = knots, col = "gray", lty = 2)

#}

5.2 Plot the optimal spline

spline0 = Zbasis_finner %*% z
dim(spline0)

## [1] 100   1

head(spline0)

##           [,1]
## [1,] -2.104562
## [2,] -2.004042
## [3,] -1.903485
## [4,] -1.802961
## [5,] -1.702538
## [6,] -1.602287

parnition is a grid of point spline0 is the optimal spline evaluated at the grid point “parnition”

#if (spline.plot == TRUE) {
  matplot(parnition, spline0, type = "l", las = 1, 
          xlab = "t", ylab = "fcenLR(density)", 
          col = "darkblue", lwd = 2, cex.lab = 1.2, cex.axis = 1.2, 
          ylim = c(min(c(min(clrf), min(spline0))), max(c(max(clrf), 
                                                          max(spline0)))), main = paste("Compositional spline of degree k =",  k - 1))

 # matpoints(t, clrf, pch = 8, col = "darkblue", cex = 1.3)
  #abline(h = 0, col = "red", lty = 2, lwd = 1) # The horizontal line
  #abline(v = knots, col = "gray", lty = 2, lwd = 1)
#}

The line with pch = 8 is the original grid points and the initial \(y = clr(f)\) points.

#if (spline.plot == TRUE) {
  matplot(parnition, spline0, type = "l", las = 1, 
          xlab = "t", ylab = "fcenLR(density)", 
          col = "darkblue", lwd = 2, cex.lab = 1.2, cex.axis = 1.2, 
          ylim = c(min(c(min(clrf), min(spline0))), max(c(max(clrf), 
                                                          max(spline0)))), main = paste("Compositional spline of degree k =",  k - 1))
  matpoints(t, clrf, pch = 8, col = "darkblue", cex = 1.3)

  #abline(h = 0, col = "red", lty = 2, lwd = 1) # The horizontal line
  #abline(v = knots, col = "gray", lty = 2, lwd = 1)
#}

Check again knots

knots

##  [1] -4.7000000 -3.6555556 -2.6111111 -1.5666667 -0.5222222  0.5222222
##  [7]  1.5666667  2.6111111  3.6555556  4.7000000

The below figure are the original code of the Talska team.

#if (spline.plot == TRUE) {
  matplot(parnition, spline0, type = "l", las = 1, 
          xlab = "t", ylab = "fcenLR(density)", 
          col = "darkblue", lwd = 2, cex.lab = 1.2, cex.axis = 1.2, 
          ylim = c(min(c(min(clrf), min(spline0))), max(c(max(clrf), 
                                                          max(spline0)))), main = paste("Compositional spline of degree k =",  k - 1))
  matpoints(t, clrf, pch = 8, col = "darkblue", cex = 1.3)
  abline(h = 0, col = "red", lty = 2, lwd = 1) # The horizontal line
  abline(v = knots, col = "gray", lty = 2, lwd = 1) # The grid vertical lines

#}

It seems done!!

6. Smoothing parameter alpha

clrf = as.matrix(clrf)
Hmat = (B %*% D %*% K) %*% solve(GG) %*% (alpha * t(K) %*% 
                                            t(D) %*% t(B) %*% W)
clrfhat = (B %*% D %*% K) %*% z
reziduals = (clrf - clrfhat)
Hmat_diag = c()
for (i in 1:length(clrf)) Hmat_diag[i] = Hmat[i, i]
Hmat_diag_mean = (1/length(clrf)) * sum(Hmat_diag)
CV = (1/length(clrf)) * sum((reziduals/(rep(1, length(Hmat_diag)) - 
                                          Hmat_diag))^2)
CV

## [1] 0.2173018

GCV = (1/length(clrf)) * (sum((reziduals)^2))/((1 - Hmat_diag_mean)^2)
J = (1 - alpha) * t(z) %*% t(U) %*% t(S) %*% M %*% S %*% 
  U %*% z + alpha * t(clrf - B %*% D %*% K %*% z) %*% W %*% 
  (clrf - B %*% D %*% K %*% z)

GCV

## [1] 0.2079652