## Replication code for ## 'rTensor: An R Package for Multidimensional Array (Tensor) Unfolding, Multiplication, and Decomposition' ## Authors: James Li, Jacob Bien, Martin Wells ## 1. Introduction library("rTensor") set.seed(1234) tnsr <- rand_tensor(rep(20, 6)) Rprof() mtx <- unfold(tnsr, row_idx = c(4, 1, 3), col_idx = c(2, 5, 6)) Rprof(NULL) summaryRprof() ## See code.m for the comparison with MATLAB. ## 2. rTensor basics library("rTensor") indices <- c(10, 20, 30, 40) arr <- array(rnorm(prod(indices)), dim = indices) tnsr <- as.tensor(arr) tnsr fnorm(tnsr) modeSum(tnsr, m = 1, drop = FALSE) modeSum(tnsr, m = 1, drop = TRUE) innerProd(tnsr, tnsr) tnsr[, 1:2, 1, 1] plot_orl(subject = 5, condition = 10) tnsr <- rand_tensor(modes = c(3, 4, 5, 6)) unfold(tnsr, row_idx = c(1, 2), col_idx = c(3, 4)) unfold(tnsr, row_idx = c(2, 3), col_idx = c(1, 4)) tnsr <- rand_tensor(modes = c(3, 4, 5, 6)) unfolded <- unfold(tnsr, row_idx = c(2, 3), col_idx = c(1, 4)) folded_back <- fold(unfolded, row_idx = c(2, 3), col_idx = c(1, 4), modes = c(3, 4, 5, 6)) identical(folded_back, tnsr) tnsr <- rand_tensor(modes = c(2, 3, 4, 5, 6)) k_unfold(tnsr, m = 2) k_unfold(tnsr, m = 4) tnsr <- rand_tensor(modes = c(2, 3, 4, 5, 6)) unfolded <- k_unfold(tnsr, m = 4) k_fold(unfolded, m = 4, modes = c(2, 3, 4, 5, 6)) identical(k_fold(unfolded, m = 4, modes = c(2, 3, 4, 5, 6)), tnsr) tnsr <- rand_tensor(modes = c(4, 6, 8, 10)) mat <- matrix(rnorm(12), ncol = 6) ttm(tnsr = tnsr, mat = mat, m = 2) mat2 <- matrix(rnorm(24), ncol = 8) ttl(tnsr = tnsr, list_mat = list(mat, mat2), ms = c(2, 3)) subject <- faces_tnsr[, , 1, ] new_positions <- c(7, 10, 1, 6, 3, 2, 5, 8, 9, 4) mat <- matrix(0, 10, 10) for (i in 1:10L) { mat[new_positions[i], i] <- 1 } subject_new <- ttm(subject, mat, m = 3) all(identical(subject_new[, , new_positions], subject[, , seq(1, 10)])) tnsr1 <- rand_tensor(modes = c(3, 4, 5)) tnsr2 <- rand_tensor(modes = c(4, 6, 5)) t_mult(tnsr1, tnsr2) ## 3. rTensor decompositions subject14 <- faces_tnsr[, , 14, ] cp1 <- cp(subject14, num_components = 50) cp2 <- cp(subject14, num_components = 10) cp1$norm_percent cp2$norm_percent plot_orl(subject = 5, condition = 10) # Generates Figure 5 greyscale <- grey(seq(0, 1, length = 256)) par(mfrow = c(1, 3), mar = c(0.1, 0.1, 1, 0.1)) image(subject14[, , 1]@data, col = greyscale) image(cp1$est[, , 1]@data, col = greyscale) image(cp2$est[, , 1]@data, col = greyscale) tnsr <- rand_tensor(modes = c(2, 4, 6, 8)) hosvd1 <- hosvd(tnsr) hosvd1$U[[1]] %*% t(hosvd1$U[[1]]) hosvd2 <- hosvd(tnsr, ranks = c(1, 2, 3, 4)) 1 - hosvd2$fnorm_resid/fnorm(tnsr) tucker1 <- tucker(faces_tnsr, ranks = c(46, 56, 35, 8)) tucker1$norm_percent tucker2 <- tucker(faces_tnsr, ranks = c(23, 28, 10, 3)) tucker2$norm_percent # Generates Figure 7 greyscale <- grey(seq(0, 1, length = 256)) par(mfrow = c(1, 3), mar = c(0.1, 0.1, 1, 0.1)) subject11 <- faces_tnsr[, , 11, ] image(subject11[, , 6]@data, col = greyscale) image(tucker1$est[, , 11, 6]@data, col = greyscale) image(tucker2$est[, , 11, 6]@data, col = greyscale) subject21 <- faces_tnsr[, , 21, ] glram1 <- tucker(subject21, ranks = c(46, 56, 10)) glram1$norm_percent glram2 <- tucker(subject21, ranks = c(23, 28, 10)) glram2$norm_percent # Generates Figure 8 greyscale <- grey(seq(0, 1, length = 256)) par(mfrow = c(1, 3), mar = c(0.1, 0.1, 1, 0.1)) image(subject21[, , 2]@data, col = greyscale) image(glram1$est[, , 2]@data, col = greyscale) image(glram2$est[, , 2]@data, col = greyscale) mpca1 <- tucker(faces_tnsr, ranks = c(46, 56, 20, 10)) mpca1$norm_percent mpca2 <- tucker(faces_tnsr, ranks = c(46, 56, 10, 10)) mpca2$norm_percent # Generates Figure 9 greyscale <- grey(seq(0, 1, length = 256)) par(mfrow = c(1, 3), mar = c(0.1, 0.1, 1, 0.1)) subject35 <- faces_tnsr[, , 35, ] image(subject35[, , 8]@data, col = greyscale) image(mpca1$est[, , 35, 8]@data, col = greyscale) image(mpca2$est[, , 35, 8]@data, col = greyscale) subject8 <- faces_tnsr[, , 8, ] pvd1 <- pvd(subject8, uranks = rep(46, 10), wranks = rep(56, 10), a = 46, b = 56) pvd2 <- pvd(subject8, uranks = rep(46, 10), wranks = rep(56, 10), a = 23, b = 28) pvd1$norm_percent pvd2$norm_percent # Generates Figure 11 greyscale <- grey(seq(0, 1, length = 256)) par(mfrow = c(1, 3), mar = c(0.1, 0.1, 1, 0.1)) image(subject8[, , 4]@data, col = greyscale) image(pvd1$est[, , 4]@data, col = greyscale) image(pvd2$est[, , 4]@data, col = greyscale) tnsr <- rand_tensor(c(4, 5, 6)) tnsr t_tnsr <- t(tnsr) t_tnsr identical(t(t_tnsr), tnsr) tnsr <- rand_tensor(c(3, 4, 5)) decomp <- t_svd(tnsr) decomp decomp$S@data recovered <- t_mult(t_mult(decomp$U, decomp$S), t(decomp$V)) mean(recovered@data - tnsr@data)