Row Kronecker product/ tensor product smooth construction

Produce model matrices or penalty matrices for a tensor product smooth from the model matrices or penalty matrices for the marginal bases of the smooth (marginals and results can be sparse). The model matrix construction uses row Kronecker products.

tensor.prod.model.matrix(X)
tensor.prod.penalties(S)
a %.% b

Arguments

X

a list of model matrices for the marginal bases of a smooth. Items can be class "matrix" or "dgCMatrix", but not a mixture of the two.

S

a list of penalties for the marginal bases of a smooth.

a

a matrix with the same number of rows as A.

b

a matrix with the same number of rows as B.

Details

If X[[1]], X[[2]] ... X[[m]] are the model matrices of the marginal bases of a tensor product smooth then the ith row of the model matrix for the whole tensor product smooth is given by X[[1]][i,]%x%X[[2]][i,]%x% ... X[[m]][i,], where %x% is the Kronecker product. Of course the routine operates column-wise, not row-wise!

A%.%B is the operator form of this `row Kronecker product'.

If S[[1]], S[[2]] ... S[[m]] are the penalty matrices for the marginal bases, and I[[1]], I[[2]] ... I[[m]] are corresponding identity matrices, each of the same dimension as its corresponding penalty, then the tensor product smooth has m associate penalties of the form:

S[[1]]%x%I[[2]]%x% ... I[[m]],

I[[1]]%x%S[[2]]%x% ... I[[m]]

...

I[[1]]%x%I[[2]]%x% ... S[[m]].

Of course it's important that the model matrices and penalty matrices are presented in the same order when constructing tensor product smooths.

Value

Either a single model matrix for a tensor product smooth (of the same class as the marginals), or a list of penalty terms for a tensor product smooth.

References

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics 62(4):1025-1036

Author

Simon N. Wood simon.wood@r-project.org

See also

te, smooth.construct.tensor.smooth.spec

Examples

require(mgcv)
## Dense row Kronecker product example...
X <- list(matrix(0:3,2,2),matrix(c(5:8,0,0),2,3))
tensor.prod.model.matrix(X)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    0    0    0   10   14    0
#> [2,]    6    8    0   18   24    0
X[[1]]%.%X[[2]]
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    0    0    0   10   14    0
#> [2,]    6    8    0   18   24    0

## sparse equivalent...
Xs <- lapply(X,as,"dgCMatrix")
tensor.prod.model.matrix(Xs)
#> 2 x 6 sparse Matrix of class "dgCMatrix"
#>                   
#> [1,] . . . 10 14 .
#> [2,] 6 8 . 18 24 .
Xs[[1]]%.%Xs[[2]]
#> 2 x 6 sparse Matrix of class "dgCMatrix"
#>                   
#> [1,] . . . 10 14 .
#> [2,] 6 8 . 18 24 .

S <- list(matrix(c(2,1,1,2),2,2),matrix(c(2,1,0,1,2,1,0,1,2),3,3))
tensor.prod.penalties(S)
#> [[1]]
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    2    0    0    1    0    0
#> [2,]    0    2    0    0    1    0
#> [3,]    0    0    2    0    0    1
#> [4,]    1    0    0    2    0    0
#> [5,]    0    1    0    0    2    0
#> [6,]    0    0    1    0    0    2
#> 
#> [[2]]
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    2    1    0    0    0    0
#> [2,]    1    2    1    0    0    0
#> [3,]    0    1    2    0    0    0
#> [4,]    0    0    0    2    1    0
#> [5,]    0    0    0    1    2    1
#> [6,]    0    0    0    0    1    2
#> 
## Sparse equivalent...
Ss <- lapply(S,as,"dgCMatrix")
tensor.prod.penalties(Ss)
#> [[1]]
#> 6 x 6 sparse Matrix of class "dgCMatrix"
#>                 
#> [1,] 2 . . 1 . .
#> [2,] . 2 . . 1 .
#> [3,] . . 2 . . 1
#> [4,] 1 . . 2 . .
#> [5,] . 1 . . 2 .
#> [6,] . . 1 . . 2
#> 
#> [[2]]
#> 6 x 6 sparse Matrix of class "dgCMatrix"
#>                 
#> [1,] 2 1 . . . .
#> [2,] 1 2 1 . . .
#> [3,] . 1 2 . . .
#> [4,] . . . 2 1 .
#> [5,] . . . 1 2 1
#> [6,] . . . . 1 2
#>