tensor_predictors/tensorPredictors/R/mtvk.R

64 lines
2.1 KiB
R

#' Matrix Times Vectorized Kronecker product
#'
#' \deqn{C = A vec(B_1\otimes ... \otimes B_r)}{%
#' C = A vec(B_1 %x% ... %x% B_r)}
#'
#' This function is equivalent to `c(A %*% c(Reduce("%x%", Bs)))`.
#'
#' @param A matrix of dimensions `n` by `pp * qq`
#' @param Bs list of matrices such that the product there Kronecker product has
#' dimensions `pp` by `qq`.
#'
#' @returns vector of length `n`
#'
#' @examples
#' n <- 17
#' p <- c(2, 5, 11)
#' q <- c(3, 7, 13)
#'
#' A <- matrix(rnorm(n * prod(p * q)), n)
#' Bs <- Map(matrix, Map(rnorm, p * q), p)
#'
#' stopifnot(all.equal(
#' c(A %*% c(Reduce(`%x%`, Bs))),
#' mtvk(A, Bs)
#' ))
#'
#' @note May be slower than `c(A %*% c(Reduce("%x%", Bs)))`.
#' @TODO C++ version using Rcpp is much faster than plain C using `R`s C API!
#'
#' @export
mtvk <- function(A, Bs) {
c(A %*% c(Reduce("%x%", Bs)))
# storage.mode(A) <- "double"
# Bs <- Map(`storage.mode<-`, Bs, list("double"))
# .Call("C_mtvk", A, Bs)
}
# n <- 17
# p <- rev(c(11, 7, 20))
# q <- rev(c(13, 5, 30))
# r <- length(p)
# A <- matrix(rnorm(n * prod(p * q)), n)
# Bs <- Map(matrix, Map(rnorm, p * q), p)
# microbenchmark::microbenchmark(
# A %*% c(Reduce("%x%", Bs)),
# mtvk(A, Bs)
# )
# gcc -I"/usr/share/R/include" -DNDEBUG -fpic -g -O2 -fdebug-prefix-map=/build/r-base-zYgbYq/r-base-4.2.1=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c mtvk.c -o mtvk.o
# g++ -I"/usr/share/R/include" -DNDEBUG -fpic -g -O2 -fdebug-prefix-map=/build/r-base-zYgbYq/r-base-4.2.1=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c mtvk.cpp -o mtvk.o \
# -std=gnu++14 -I"/usr/local/lib/R/site-library/Rcpp/include" -I"/home/loki/Work/tensorPredictors/wip"
# g++ -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions -Wl,-z,relro -o sourceCpp_2.so mtvk.o -L/usr/lib/R/lib -lR \
# -std=gnu++14
# gcc -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions -Wl,-z,relro -o tensorPredictors.so -L/usr/lib/R/lib -lR \
# init.o mcrossprod.o mtvk.o poi.o ttm.o -lblas -lgfortran -lm -lquadmath