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

100 values based on normal density

#--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) 
f.fcenLR

##   [1] -3.272258689 -3.075923167 -2.883594493 -2.695272667 -2.510957688
##   [6] -2.330649556 -2.154348271 -1.982053834 -1.813766245 -1.649485502
##  [11] -1.489211607 -1.332944560 -1.180684359 -1.032431007 -0.888184501
##  [16] -0.747944843 -0.611712032 -0.479486069 -0.351266953 -0.227054685
##  [21] -0.106849263  0.009349311  0.121541037  0.229725916  0.333903948
##  [26]  0.434075132  0.530239469  0.622396959  0.710547601  0.794691396
##  [31]  0.874828343  0.950958443  1.023081696  1.091198101  1.155307659
##  [36]  1.215410370  1.271506233  1.323595249  1.371677418  1.415752739
##  [41]  1.455821212  1.491882839  1.523937618  1.551985549  1.576026634
##  [46]  1.596060871  1.612088260  1.624108802  1.632122497  1.636129344
##  [51]  1.636129344  1.632122497  1.624108802  1.612088260  1.596060871
##  [56]  1.576026634  1.551985549  1.523937618  1.491882839  1.455821212
##  [61]  1.415752739  1.371677418  1.323595249  1.271506233  1.215410370
##  [66]  1.155307659  1.091198101  1.023081696  0.950958443  0.874828343
##  [71]  0.794691396  0.710547601  0.622396959  0.530239469  0.434075132
##  [76]  0.333903948  0.229725916  0.121541037  0.009349311 -0.106849263
##  [81] -0.227054685 -0.351266953 -0.479486069 -0.611712032 -0.747944843
##  [86] -0.888184501 -1.032431007 -1.180684359 -1.332944560 -1.489211607
##  [91] -1.649485502 -1.813766245 -1.982053834 -2.154348271 -2.330649556
##  [96] -2.510957688 -2.695272667 -2.883594493 -3.075923167 -3.272258689

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

# 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
w

##   [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [16] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [31] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [46] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [61] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [76] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##  [91] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01

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
D_mat

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

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
##  [7,]    0    0    0    0    0    0   -1    1    0     0     0     0
##  [8,]    0    0    0    0    0    0    0   -1    1     0     0     0
##  [9,]    0    0    0    0    0    0    0    0   -1     1     0     0
## [10,]    0    0    0    0    0    0    0    0    0    -1     1     0
## [11,]    0    0    0    0    0    0    0    0    0     0    -1     1

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
D_mat

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

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
##  [7,]    0    0    0    0    0    0   -1    1    0     0     0
##  [8,]    0    0    0    0    0    0    0   -1    1     0     0
##  [9,]    0    0    0    0    0    0    0    0   -1     1     0
## [10,]    0    0    0    0    0    0    0    0    0    -1     1

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)
    }
  }
  M

##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
##  [1,] 0.3495868 0.1726354 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
##  [2,] 0.1726354 0.6991736 0.1726354 0.0000000 0.0000000 0.0000000 0.0000000
##  [3,] 0.0000000 0.1726354 0.6991736 0.1726354 0.0000000 0.0000000 0.0000000
##  [4,] 0.0000000 0.0000000 0.1726354 0.6991736 0.1726354 0.0000000 0.0000000
##  [5,] 0.0000000 0.0000000 0.0000000 0.1726354 0.6991736 0.1726354 0.0000000
##  [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.1726354 0.6991736 0.1726354
##  [7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1726354 0.6991736
##  [8,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1726354
##  [9,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [10,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
##            [,8]      [,9]     [,10]
##  [1,] 0.0000000 0.0000000 0.0000000
##  [2,] 0.0000000 0.0000000 0.0000000
##  [3,] 0.0000000 0.0000000 0.0000000
##  [4,] 0.0000000 0.0000000 0.0000000
##  [5,] 0.0000000 0.0000000 0.0000000
##  [6,] 0.0000000 0.0000000 0.0000000
##  [7,] 0.1726354 0.0000000 0.0000000
##  [8,] 0.6991736 0.1726354 0.0000000
##  [9,] 0.1726354 0.6991736 0.1726354
## [10,] 0.0000000 0.1726354 0.3495868

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
D

##           [,1]     [,2]     [,3]      [,4]      [,5]      [,6]      [,7]
##  [1,] 3.829787 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
##  [2,] 0.000000 1.914894 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
##  [3,] 0.000000 0.000000 1.276596 0.0000000 0.0000000 0.0000000 0.0000000
##  [4,] 0.000000 0.000000 0.000000 0.9574468 0.0000000 0.0000000 0.0000000
##  [5,] 0.000000 0.000000 0.000000 0.0000000 0.9574468 0.0000000 0.0000000
##  [6,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.9574468 0.0000000
##  [7,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.9574468
##  [8,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
##  [9,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
## [10,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
## [11,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
## [12,] 0.000000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000
##            [,8]      [,9]    [,10]    [,11]    [,12]
##  [1,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [2,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [3,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [4,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [5,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [6,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [7,] 0.0000000 0.0000000 0.000000 0.000000 0.000000
##  [8,] 0.9574468 0.0000000 0.000000 0.000000 0.000000
##  [9,] 0.0000000 0.9574468 0.000000 0.000000 0.000000
## [10,] 0.0000000 0.0000000 1.276596 0.000000 0.000000
## [11,] 0.0000000 0.0000000 0.000000 1.914894 0.000000
## [12,] 0.0000000 0.0000000 0.000000 0.000000 3.829787

K

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

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

U = D %*% K
U

##            [,1]      [,2]       [,3]       [,4]       [,5]       [,6]
##  [1,]  3.829787  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
##  [2,] -1.914894  1.914894  0.0000000  0.0000000  0.0000000  0.0000000
##  [3,]  0.000000 -1.276596  1.2765957  0.0000000  0.0000000  0.0000000
##  [4,]  0.000000  0.000000 -0.9574468  0.9574468  0.0000000  0.0000000
##  [5,]  0.000000  0.000000  0.0000000 -0.9574468  0.9574468  0.0000000
##  [6,]  0.000000  0.000000  0.0000000  0.0000000 -0.9574468  0.9574468
##  [7,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000 -0.9574468
##  [8,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
##  [9,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
## [10,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
## [11,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
## [12,]  0.000000  0.000000  0.0000000  0.0000000  0.0000000  0.0000000
##             [,7]       [,8]       [,9]     [,10]     [,11]
##  [1,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [2,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [3,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [4,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [5,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [6,]  0.0000000  0.0000000  0.0000000  0.000000  0.000000
##  [7,]  0.9574468  0.0000000  0.0000000  0.000000  0.000000
##  [8,] -0.9574468  0.9574468  0.0000000  0.000000  0.000000
##  [9,]  0.0000000 -0.9574468  0.9574468  0.000000  0.000000
## [10,]  0.0000000  0.0000000 -1.2765957  1.276596  0.000000
## [11,]  0.0000000  0.0000000  0.0000000 -1.914894  1.914894
## [12,]  0.0000000  0.0000000  0.0000000  0.000000 -3.829787

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.1
GG = t(U) %*% ((1 - alpha) * t(S) %*% M %*% S + alpha * t(B) %*%  W %*% B) %*% U
GG

##               [,1]         [,2]        [,3]       [,4]       [,5]       [,6]
##  [1,]  401.7740172 -192.9462243  22.6010279  4.6493415 -0.3590781  0.0000000
##  [2,] -192.9462243  102.2455183 -18.9656547 -0.5386241  0.5901640 -0.1595903
##  [3,]   22.6010279  -18.9656547  10.4509909 -4.2926564  0.8739103  0.2730583
##  [4,]    4.6493415   -0.5386241  -4.2926564  6.1091276 -4.1595600  0.9915281
##  [5,]   -0.3590781    0.5901640   0.8739103 -4.1595600  6.1091276 -4.1595600
##  [6,]    0.0000000   -0.1595903   0.2730583  0.9915281 -4.1595600  6.1091276
##  [7,]    0.0000000    0.0000000  -0.1196927  0.2331608  0.9915281 -4.1595600
##  [8,]    0.0000000    0.0000000   0.0000000 -0.1196927  0.2331608  0.9915281
##  [9,]    0.0000000    0.0000000   0.0000000  0.0000000 -0.1196927  0.2730583
## [10,]    0.0000000    0.0000000   0.0000000  0.0000000  0.0000000 -0.1595903
## [11,]    0.0000000    0.0000000   0.0000000  0.0000000  0.0000000  0.0000000
##             [,7]       [,8]        [,9]        [,10]        [,11]
##  [1,]  0.0000000  0.0000000   0.0000000    0.0000000    0.0000000
##  [2,]  0.0000000  0.0000000   0.0000000    0.0000000    0.0000000
##  [3,] -0.1196927  0.0000000   0.0000000    0.0000000    0.0000000
##  [4,]  0.2331608 -0.1196927   0.0000000    0.0000000    0.0000000
##  [5,]  0.9915281  0.2331608  -0.1196927    0.0000000    0.0000000
##  [6,] -4.1595600  0.9915281   0.2730583   -0.1595903    0.0000000
##  [7,]  6.1091276 -4.1595600   0.8739103    0.5901640   -0.3590781
##  [8,] -4.1595600  6.1091276  -4.2926564   -0.5386241    4.6493415
##  [9,]  0.8739103 -4.2926564  10.4509909  -18.9656547   22.6010279
## [10,]  0.5901640 -0.5386241 -18.9656547  102.2455183 -192.9462243
## [11,] -0.3590781  4.6493415  22.6010279 -192.9462243  401.7740172

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

Error!! Check again each matrix

alpha * t(K) %*% t(D) %*% t(B) %*% W # work!!

##              [,1]         [,2]         [,3]         [,4]          [,5]
##  [1,] 0.003829787 2.423829e-03 1.317838e-03 4.794427e-04 -1.237272e-04
##  [2,] 0.000000000 4.386000e-04 7.235002e-04 8.750819e-04  9.137267e-04
##  [3,] 0.000000000 1.482648e-05 5.530956e-05 1.154547e-04  1.892674e-04
##  [4,] 0.000000000 1.198907e-07 9.591253e-07 3.237048e-06  7.673002e-06
##  [5,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
##  [6,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
##  [7,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
##  [8,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
##  [9,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
## [10,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
## [11,] 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,] -5.240421e-04 -7.538725e-04 -0.0008455888 -8.315616e-04 -7.441613e-04
##  [2,]  8.598158e-04  7.337308e-04  0.0005558531  3.465639e-04  1.262449e-04
##  [3,]  2.707531e-04  3.539172e-04  0.0004327653  5.013028e-04  5.535352e-04
##  [4,]  1.498633e-05  2.589638e-05  0.0000411225  6.138402e-05  8.740029e-05
##  [5,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [6,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [7,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [8,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##               [,11]         [,12]         [,13]         [,14]         [,15]
##  [1,] -6.157584e-04 -0.0004787234 -3.596720e-04 -2.622009e-04 -1.841521e-04
##  [2,] -8.472273e-05 -0.0002659574 -4.013939e-04 -4.922311e-04 -5.439839e-04
##  [3,]  5.834679e-04  0.0005851064  5.543744e-04  4.948687e-04  4.121042e-04
##  [4,]  1.198907e-04  0.0001595745  2.065716e-04  2.586042e-04  3.127947e-04
##  [5,]  0.000000e+00  0.0000000000  1.198907e-07  9.591253e-07  3.237048e-06
##  [6,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
##  [7,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
##  [8,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00
##               [,16]         [,17]         [,18]         [,19]         [,20]
##  [1,] -1.233675e-04 -7.768915e-05 -4.495900e-05 -2.301901e-05 -9.711143e-06
##  [2,] -5.621673e-04 -5.522963e-04 -5.198859e-04 -4.704509e-04 -4.095065e-04
##  [3,]  3.115958e-04  1.988586e-04  7.940758e-05 -4.124239e-05 -1.575763e-04
##  [4,]  3.662660e-04  4.161405e-04  4.595409e-04  4.935898e-04  5.154099e-04
##  [5,]  7.673002e-06  1.498633e-05  2.589638e-05  4.112250e-05  6.138402e-05
##  [6,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [7,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [8,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,21]         [,22]         [,23]         [,24]         [,25]
##  [1,] -2.877376e-06 -3.596720e-07  0.0000000000  0.000000e+00  0.000000e+00
##  [2,] -3.425676e-04 -2.751491e-04 -0.0002127660 -1.598542e-04 -1.165337e-04
##  [3,] -2.640792e-04 -3.552360e-04 -0.0004255319 -4.708905e-04 -4.929904e-04
##  [4,]  5.221238e-04  5.108541e-04  0.0004787234  4.240533e-04  3.499608e-04
##  [5,]  8.740029e-05  1.198907e-04  0.0001595745  2.065716e-04  2.586042e-04
##  [6,]  0.000000e+00  0.000000e+00  0.0000000000  1.198907e-07  9.591253e-07
##  [7,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [8,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##               [,26]         [,27]         [,28]         [,29]         [,30]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,] -8.184536e-05 -5.483000e-05 -3.452851e-05 -1.998178e-05 -1.023067e-05
##  [3,] -4.949486e-04 -4.798823e-04 -4.509088e-04 -4.111450e-04 -3.637083e-04
##  [4,]  2.607622e-04  1.607734e-04  5.431047e-05 -5.431047e-05 -1.607734e-04
##  [5,]  3.127947e-04  3.662660e-04  4.161405e-04  4.595409e-04  4.935898e-04
##  [6,]  3.237048e-06  7.673002e-06  1.498633e-05  2.589638e-05  4.112250e-05
##  [7,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [8,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,31]         [,32]         [,33]         [,34]         [,35]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,] -4.316064e-06 -1.278834e-06 -1.598542e-07 -1.635523e-50  0.000000e+00
##  [3,] -3.117157e-04 -2.582844e-04 -2.065316e-04 -1.595745e-04 -1.198907e-04
##  [4,] -2.607622e-04 -3.499608e-04 -4.240533e-04 -4.787234e-04 -5.108541e-04
##  [5,]  5.154099e-04  5.221238e-04  5.108541e-04  4.787234e-04  4.240533e-04
##  [6,]  6.138402e-05  8.740029e-05  1.198907e-04  1.595745e-04  2.065716e-04
##  [7,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.198907e-07
##  [8,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,36]         [,37]         [,38]         [,39]         [,40]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,] -8.740029e-05 -6.138402e-05 -4.112250e-05 -2.589638e-05 -1.498633e-05
##  [4,] -5.221238e-04 -5.154099e-04 -4.935898e-04 -4.595409e-04 -4.161405e-04
##  [5,]  3.499608e-04  2.607622e-04  1.607734e-04  5.431047e-05 -5.431047e-05
##  [6,]  2.586042e-04  3.127947e-04  3.662660e-04  4.161405e-04  4.595409e-04
##  [7,]  9.591253e-07  3.237048e-06  7.673002e-06  1.498633e-05  2.589638e-05
##  [8,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,41]         [,42]         [,43]         [,44]         [,45]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
##  [3,] -7.673002e-06 -3.237048e-06 -9.591253e-07 -1.198907e-07  0.0000000000
##  [4,] -3.662660e-04 -3.127947e-04 -2.586042e-04 -2.065716e-04 -0.0001595745
##  [5,] -1.607734e-04 -2.607622e-04 -3.499608e-04 -4.240533e-04 -0.0004787234
##  [6,]  4.935898e-04  5.154099e-04  5.221238e-04  5.108541e-04  0.0004787234
##  [7,]  4.112250e-05  6.138402e-05  8.740029e-05  1.198907e-04  0.0001595745
##  [8,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000
##               [,46]         [,47]         [,48]         [,49]         [,50]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,] -1.198907e-04 -8.740029e-05 -6.138402e-05 -4.112250e-05 -2.589638e-05
##  [5,] -5.108541e-04 -5.221238e-04 -5.154099e-04 -4.935898e-04 -4.595409e-04
##  [6,]  4.240533e-04  3.499608e-04  2.607622e-04  1.607734e-04  5.431047e-05
##  [7,]  2.065716e-04  2.586042e-04  3.127947e-04  3.662660e-04  4.161405e-04
##  [8,]  1.198907e-07  9.591253e-07  3.237048e-06  7.673002e-06  1.498633e-05
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,51]         [,52]         [,53]         [,54]         [,55]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,] -1.498633e-05 -7.673002e-06 -3.237048e-06 -9.591253e-07 -1.198907e-07
##  [5,] -4.161405e-04 -3.662660e-04 -3.127947e-04 -2.586042e-04 -2.065716e-04
##  [6,] -5.431047e-05 -1.607734e-04 -2.607622e-04 -3.499608e-04 -4.240533e-04
##  [7,]  4.595409e-04  4.935898e-04  5.154099e-04  5.221238e-04  5.108541e-04
##  [8,]  2.589638e-05  4.112250e-05  6.138402e-05  8.740029e-05  1.198907e-04
##  [9,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,56]         [,57]         [,58]         [,59]         [,60]
##  [1,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,] -0.0001595745 -1.198907e-04 -8.740029e-05 -6.138402e-05 -4.112250e-05
##  [6,] -0.0004787234 -5.108541e-04 -5.221238e-04 -5.154099e-04 -4.935898e-04
##  [7,]  0.0004787234  4.240533e-04  3.499608e-04  2.607622e-04  1.607734e-04
##  [8,]  0.0001595745  2.065716e-04  2.586042e-04  3.127947e-04  3.662660e-04
##  [9,]  0.0000000000  1.198907e-07  9.591253e-07  3.237048e-06  7.673002e-06
## [10,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.0000000000  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,61]         [,62]         [,63]         [,64]         [,65]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,] -2.589638e-05 -1.498633e-05 -7.673002e-06 -3.237048e-06 -9.591253e-07
##  [6,] -4.595409e-04 -4.161405e-04 -3.662660e-04 -3.127947e-04 -2.586042e-04
##  [7,]  5.431047e-05 -5.431047e-05 -1.607734e-04 -2.607622e-04 -3.499608e-04
##  [8,]  4.161405e-04  4.595409e-04  4.935898e-04  5.154099e-04  5.221238e-04
##  [9,]  1.498633e-05  2.589638e-05  4.112250e-05  6.138402e-05  8.740029e-05
## [10,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,66]         [,67]         [,68]         [,69]         [,70]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,] -1.198907e-07 -9.813141e-50  0.000000e+00  0.000000e+00  0.000000e+00
##  [6,] -2.065716e-04 -1.595745e-04 -1.198907e-04 -8.740029e-05 -6.138402e-05
##  [7,] -4.240533e-04 -4.787234e-04 -5.108541e-04 -5.221238e-04 -5.154099e-04
##  [8,]  5.108541e-04  4.787234e-04  4.240533e-04  3.499608e-04  2.607622e-04
##  [9,]  1.198907e-04  1.595745e-04  2.065316e-04  2.582844e-04  3.117157e-04
## [10,]  0.000000e+00  0.000000e+00  1.598542e-07  1.278834e-06  4.316064e-06
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,71]         [,72]         [,73]         [,74]         [,75]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [6,] -4.112250e-05 -2.589638e-05 -1.498633e-05 -7.673002e-06 -3.237048e-06
##  [7,] -4.935898e-04 -4.595409e-04 -4.161405e-04 -3.662660e-04 -3.127947e-04
##  [8,]  1.607734e-04  5.431047e-05 -5.431047e-05 -1.607734e-04 -2.607622e-04
##  [9,]  3.637083e-04  4.111450e-04  4.509088e-04  4.798823e-04  4.949486e-04
## [10,]  1.023067e-05  1.998178e-05  3.452851e-05  5.483000e-05  8.184536e-05
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##               [,76]         [,77]         [,78]         [,79]         [,80]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [6,] -9.591253e-07 -1.198907e-07 -9.813141e-50  0.000000e+00  0.000000e+00
##  [7,] -2.586042e-04 -2.065716e-04 -1.595745e-04 -1.198907e-04 -8.740029e-05
##  [8,] -3.499608e-04 -4.240533e-04 -4.787234e-04 -5.108541e-04 -5.221238e-04
##  [9,]  4.929904e-04  4.708905e-04  4.255319e-04  3.552360e-04  2.640792e-04
## [10,]  1.165337e-04  1.598542e-04  2.127660e-04  2.751491e-04  3.425676e-04
## [11,]  0.000000e+00  0.000000e+00  0.000000e+00  3.596720e-07  2.877376e-06
##               [,81]         [,82]         [,83]         [,84]         [,85]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [5,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [6,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [7,] -6.138402e-05 -4.112250e-05 -2.589638e-05 -1.498633e-05 -7.673002e-06
##  [8,] -5.154099e-04 -4.935898e-04 -4.595409e-04 -4.161405e-04 -3.662660e-04
##  [9,]  1.575763e-04  4.124239e-05 -7.940758e-05 -1.988586e-04 -3.115958e-04
## [10,]  4.095065e-04  4.704509e-04  5.198859e-04  5.522963e-04  5.621673e-04
## [11,]  9.711143e-06  2.301901e-05  4.495900e-05  7.768915e-05  1.233675e-04
##               [,86]         [,87]         [,88]         [,89]         [,90]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [5,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [6,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00
##  [7,] -3.237048e-06 -9.591253e-07 -1.198907e-07  0.0000000000  0.000000e+00
##  [8,] -3.127947e-04 -2.586042e-04 -2.065716e-04 -0.0001595745 -1.198907e-04
##  [9,] -4.121042e-04 -4.948687e-04 -5.543744e-04 -0.0005851064 -5.834679e-04
## [10,]  5.439839e-04  4.922311e-04  4.013939e-04  0.0002659574  8.472273e-05
## [11,]  1.841521e-04  2.622009e-04  3.596720e-04  0.0004787234  6.157584e-04
##               [,91]         [,92]         [,93]         [,94]         [,95]
##  [1,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [2,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [3,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [4,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [5,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [6,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [7,]  0.000000e+00  0.000000e+00  0.0000000000  0.000000e+00  0.000000e+00
##  [8,] -8.740029e-05 -6.138402e-05 -0.0000411225 -2.589638e-05 -1.498633e-05
##  [9,] -5.535352e-04 -5.013028e-04 -0.0004327653 -3.539172e-04 -2.707531e-04
## [10,] -1.262449e-04 -3.465639e-04 -0.0005558531 -7.337308e-04 -8.598158e-04
## [11,]  7.441613e-04  8.315616e-04  0.0008455888  7.538725e-04  5.240421e-04
##               [,96]         [,97]         [,98]         [,99]       [,100]
##  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [2,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [3,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [4,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [5,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [6,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [7,]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000000
##  [8,] -7.673002e-06 -3.237048e-06 -9.591253e-07 -1.198907e-07  0.000000000
##  [9,] -1.892674e-04 -1.154547e-04 -5.530956e-05 -1.482648e-05  0.000000000
## [10,] -9.137267e-04 -8.750819e-04 -7.235002e-04 -4.386000e-04  0.000000000
## [11,]  1.237272e-04 -4.794427e-04 -1.317838e-03 -2.423829e-03 -0.003829787

See clrf, i,e y

clrf

##   [1] -3.272258689 -3.075923167 -2.883594493 -2.695272667 -2.510957688
##   [6] -2.330649556 -2.154348271 -1.982053834 -1.813766245 -1.649485502
##  [11] -1.489211607 -1.332944560 -1.180684359 -1.032431007 -0.888184501
##  [16] -0.747944843 -0.611712032 -0.479486069 -0.351266953 -0.227054685
##  [21] -0.106849263  0.009349311  0.121541037  0.229725916  0.333903948
##  [26]  0.434075132  0.530239469  0.622396959  0.710547601  0.794691396
##  [31]  0.874828343  0.950958443  1.023081696  1.091198101  1.155307659
##  [36]  1.215410370  1.271506233  1.323595249  1.371677418  1.415752739
##  [41]  1.455821212  1.491882839  1.523937618  1.551985549  1.576026634
##  [46]  1.596060871  1.612088260  1.624108802  1.632122497  1.636129344
##  [51]  1.636129344  1.632122497  1.624108802  1.612088260  1.596060871
##  [56]  1.576026634  1.551985549  1.523937618  1.491882839  1.455821212
##  [61]  1.415752739  1.371677418  1.323595249  1.271506233  1.215410370
##  [66]  1.155307659  1.091198101  1.023081696  0.950958443  0.874828343
##  [71]  0.794691396  0.710547601  0.622396959  0.530239469  0.434075132
##  [76]  0.333903948  0.229725916  0.121541037  0.009349311 -0.106849263
##  [81] -0.227054685 -0.351266953 -0.479486069 -0.611712032 -0.747944843
##  [86] -0.888184501 -1.032431007 -1.180684359 -1.332944560 -1.489211607
##  [91] -1.649485502 -1.813766245 -1.982053834 -2.154348271 -2.330649556
##  [96] -2.510957688 -2.695272667 -2.883594493 -3.075923167 -3.272258689

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

## [1]  11 100

length(clrf)

## [1] 100

clrf

##   [1] -3.272258689 -3.075923167 -2.883594493 -2.695272667 -2.510957688
##   [6] -2.330649556 -2.154348271 -1.982053834 -1.813766245 -1.649485502
##  [11] -1.489211607 -1.332944560 -1.180684359 -1.032431007 -0.888184501
##  [16] -0.747944843 -0.611712032 -0.479486069 -0.351266953 -0.227054685
##  [21] -0.106849263  0.009349311  0.121541037  0.229725916  0.333903948
##  [26]  0.434075132  0.530239469  0.622396959  0.710547601  0.794691396
##  [31]  0.874828343  0.950958443  1.023081696  1.091198101  1.155307659
##  [36]  1.215410370  1.271506233  1.323595249  1.371677418  1.415752739
##  [41]  1.455821212  1.491882839  1.523937618  1.551985549  1.576026634
##  [46]  1.596060871  1.612088260  1.624108802  1.632122497  1.636129344
##  [51]  1.636129344  1.632122497  1.624108802  1.612088260  1.596060871
##  [56]  1.576026634  1.551985549  1.523937618  1.491882839  1.455821212
##  [61]  1.415752739  1.371677418  1.323595249  1.271506233  1.215410370
##  [66]  1.155307659  1.091198101  1.023081696  0.950958443  0.874828343
##  [71]  0.794691396  0.710547601  0.622396959  0.530239469  0.434075132
##  [76]  0.333903948  0.229725916  0.121541037  0.009349311 -0.106849263
##  [81] -0.227054685 -0.351266953 -0.479486069 -0.611712032 -0.747944843
##  [86] -0.888184501 -1.032431007 -1.180684359 -1.332944560 -1.489211607
##  [91] -1.649485502 -1.813766245 -1.982053834 -2.154348271 -2.330649556
##  [96] -2.510957688 -2.695272667 -2.883594493 -3.075923167 -3.272258689