Description Usage Arguments Details Value Author(s) References Examples
Implements the DeBoor (1979) algorithm to compute the product A * y
, where A
is the (large) kronecker product of M matrices mat_1,...,mat_M, each of dimension n_i x n_i,i=1,…,M, and y
is a compatible vector of length N = n_1*n_2*…*n_M
1 | armamax(y,matrices)
|
y |
a numeric vector |
matrices |
A list of square matrices |
This function is useful for approximation high dimensional functional spaces with basis functions. Suppose we want to approximate an M-dimensional object f that maps R^M into R, on the tensor product A
of univariate grids x_i of length n_i,i=1,…,M each. Denoting the size n_i x n_i basis matrix for dimension i by mat_i, this kronecker product is A <- kronecker(mat_M,kronecker(mat_M-1,...,kronecker(mat_2,mat_1))...)
.
Approximating coefficients coefs
may be obtained by solving the system A * coefs = y
, where y_{i,j,...,k} = f(x_{1,i},x_{2,j},…,x_{M,k}) are function values at the grid. Beyond a certain number of dimensions M, or number of data points n_i this is infeasible. Firstly because A
does not fit in memory, and secondly because solving the system is very costly. armakron
implements the efficient deBoor (1979) algorithm to compute coefs
. The kronecker product is never formed, and the calculation involves a series of repeated matrix multiplications. To obtain coefficients, it is most efficient to supply inverse basis matrices. This algorithm is orders of magnitude faster than for example setting up A
as a product of sparse matrices (as is the case with spline basis), and to use solve(A,y)
from a sparse library.
a N x 1 matrix of approximating spline coefficients
<florian.oswald@gmail.com>
C. De Boor. Efficient computer manipulation of tensor products. ACM Transactions on Mathematical Software (TOMS), 5(2):173-182, 1979
M. Miranda and P. Fackler. Applied computational economics and finance. MIT Press, 2002
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | ## Not run:
x1 <- matrix(rnorm(25),nrow=5,ncol=5)
x2 <- matrix(rnorm(9),nrow=3,ncol=3)
x3 <- matrix(rnorm(4),nrow=2,ncol=2)
y <- rnorm(5*3*2)
R.kron <- kronecker(x3,kronecker(x2,x1)) # form R kronecker product. impossible beyond a certain size
R.result <- R.kron %*% y
Arma.result <- armakron(y=y,list(x1,x2,x3)) # armakron() computes result directly without forming kronecker product
all.equal(R.result,Arma.result) # TRUE
# spline example: estimate spline coefficients on kronecker product of 3 univariate spline basis x,y and z.
library(splines)
# degrees
dx <- 3
dy <- 3
dz <- 3
# knots
kx <- c(rep(0,times=dx),seq(0,1,length=17),rep(1,times=dx))
ky <- c(rep(1,times=dy),seq(1,10,length=8),rep(10,times=dy))
kz <- c(rep(-3,times=dz),seq(-3,3,length=11),rep(3,times=dz))
# evaluation points: choose such that # of datapoints is equal # of coefficients
# the algorithm works ONLY with square systems, i.e. number of data points must be equal number of spline coefficients!
x <- seq(0,1,length=length(kx)-dx-1)
y <- seq(1,10,length=length(ky)-dy-1)
z <- seq(-3,3,length=length(kz)-dz-1)
# basis matrices
X <- splineDesign(x=x,knots=kx,ord=dx+1)
Y <- splineDesign(x=y,knots=ky,ord=dy+1)
Z <- splineDesign(x=z,knots=kz,ord=dz+1)
# data generating process
dgp <- function(x) x[1]^2 + x[2]^3 + x[1]*x[2]*x[3]
# create sample of data: note that ordering of data is crucial. grid order and order implied by kronecker product must cooincide.
vals <- apply(expand.grid(x,y,z),1,dgp)
# plot a slice
persp(array(vals,dim=c(length(x),length(y),length(z)))[ , ,5],theta=100)
# estimate spline coefficients to be able to evaluate dgp off the grid
# the problem is:
# kronprod %*% coefs = vals, where kronprod = kronecker(Z,kronecker(Y,X))
# coefs = inv(kronprod) %*% vals
#
# often the forming of kronprod is infeasible because of memory limitations. If not that, then computing the inverse
# is very costly. The algorithm in armakron never forms the kronprod. Instead we give it the inverse matrices one by one:
library(MASS) # adds ginv()
coefs <- armakron(y=vals,matrices=list(ginv(X),ginv(Y),ginv(Z)))
# check against R
kron <- kronecker(Z,kronecker(Y,X))
r.coefs <- solve(kron,vals) # takes couple of seconds already.
all.equal(as.numeric(coefs),r.coefs) # TRUE
# predict the grid.
pred <- armakron(y=coefs,matrices=list(X,Y,Z))
all.equal(as.numeric(pred),vals) # TRUE
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.