2. Tensor arithmetic by DelayedTensor
26 October 2021
Package
DelayedTensor 1.0.0
Authors: Koki Tsuyuzaki [aut, cre]
Last modified: 2021-10-26 15:33:21
Compiled: Tue Oct 26 17:15:40 2021
1 Setting
suppressPackageStartupMessages(library("DelayedTensor"))
suppressPackageStartupMessages(library("DelayedArray"))
suppressPackageStartupMessages(library("HDF5Array"))
suppressPackageStartupMessages(library("DelayedRandomArray"))
darr1 <- RandomUnifArray(c(2,3,4))
darr2 <- RandomUnifArray(c(2,3,4))There are several settings in DelayedTensor.
First, the sparsity of the intermediate DelayedArray objects
calculated inside DelayedTensor is set by setSparse.
Note that the sparse mode is experimental.
Whether it contributes to higher speed and lower memory is quite dependent on the sparsity of the DelayedArray, and the current implementation does not recognize the block size, which may cause out-of-memory errors, when the data is extremely huge.
Here, we specify as.sparse as FALSE
(this is also the default value for now).
DelayedTensor::setSparse(as.sparse=FALSE)Next, the verbose message is suppressed by setVerbose.
This is useful when we want to monitor the calculation process.
Here we specify as.verbose as FALSE
(this is also the default value for now).
DelayedTensor::setVerbose(as.verbose=FALSE)The block size of block processing is specified by setAutoBlockSize.
When the sparse mode is off, all the functions of DelayedTensor
are performed as block processing,
in which each block vector/matrix/tensor is expanded to memory space
from on-disk file incrementally so as not to exceed the specified size.
Here, we specify the block size as 1E+8.
setAutoBlockSize(size=1E+8)## automatic block size set to 1e+08 bytes (was 1e+08)Finally, the temporal directory to store the intermediate HDF5 files during running DelayedTensor is specified by setHDF5DumpDir.
Note that in many systems the /var directory has the storage limitation, so if there is no enough space, user should specify the other directory.
# tmpdir <- paste(sample(c(letters,1:9), 10), collapse="")
# dir.create(tmpdir, recursive=TRUE))
tmpdir <- tempdir()
setHDF5DumpDir(tmpdir)These specified values are also extracted by each getter function.
DelayedTensor::getSparse()## $delayedtensor.sparse
## [1] FALSEDelayedTensor::getVerbose()## $delayedtensor.verbose
## [1] FALSEgetAutoBlockSize()## [1] 1e+08getHDF5DumpDir()## [1] "/tmp/RtmpTsz5op"2 Tensor Arithmetic Operations
2.1 Unfold/Fold Operations
Unfold (a.k.a. matricizing) operations are used to reshape a tensor into a matrix.
Figure 1: Unfold/Fold Operasions
In unfold, row_idx and col_idx are specified to set which modes are used
as the row/column.
dmat1 <- DelayedTensor::unfold(darr1, row_idx=c(1,2), col_idx=3)
dmat1## <6 x 4> matrix of class HDF5Matrix and type "double":
## [,1] [,2] [,3] [,4]
## [1,] 0.9403318 0.7022740 0.2367886 0.7467352
## [2,] 0.5682368 0.7879356 0.3669046 0.9762840
## [3,] 0.3605710 0.6461899 0.8715041 0.5052938
## [4,] 0.4407712 0.4724037 0.7793043 0.8901269
## [5,] 0.9857139 0.6971698 0.1015803 0.7391881
## [6,] 0.9169803 0.9595190 0.7256086 0.9145331fold is the inverse operation of unfold, which is used to reshape
a matrix into a tensor.
In fold, row_idx/col_idx are specified to set which modes correspond
the row/column of the output tensor and modes
is specified to set the mode of the output tensor.
dmat1_to_darr1 <- DelayedTensor::fold(dmat1,
row_idx=c(1,2), col_idx=3, modes=dim(darr1))
dmat1_to_darr1## <2 x 3 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 0.9403318 0.3605710 0.9857139
## [2,] 0.5682368 0.4407712 0.9169803
##
## ,,2
## [,1] [,2] [,3]
## [1,] 0.7022740 0.6461899 0.6971698
## [2,] 0.7879356 0.4724037 0.9595190
##
## ,,3
## [,1] [,2] [,3]
## [1,] 0.2367886 0.8715041 0.1015803
## [2,] 0.3669046 0.7793043 0.7256086
##
## ,,4
## [,1] [,2] [,3]
## [1,] 0.7467352 0.5052938 0.7391881
## [2,] 0.9762840 0.8901269 0.9145331identical(as.array(darr1), as.array(dmat1_to_darr1))## [1] TRUEThere are some wrapper functions of unfold and fold.
For example, in k_unfold, mode m is used as the row, and the other modes
are is used as the column.
k_fold is the inverse operation of k_unfold.
dmat2 <- DelayedTensor::k_unfold(darr1, m=1)
dmat2_to_darr1 <- k_fold(dmat2, m=1, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat2_to_darr1))## [1] TRUEdmat3 <- DelayedTensor::k_unfold(darr1, m=2)
dmat3_to_darr1 <- k_fold(dmat3, m=2, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat3_to_darr1))## [1] TRUEdmat4 <- DelayedTensor::k_unfold(darr1, m=3)
dmat4_to_darr1 <- k_fold(dmat4, m=3, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat4_to_darr1))## [1] TRUEIn rs_unfold, mode m is used as the row, and the other modes
are is used as the column.
rs_fold and rs_unfold also perform the same operations.
On the other hand, cs_unfold specifies the mode m as the column
and the other modes are specified as the column.
cs_fold is the inverse operation of cs_unfold.
dmat8 <- DelayedTensor::cs_unfold(darr1, m=1)
dmat8_to_darr1 <- DelayedTensor::cs_fold(dmat8, m=1, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat8_to_darr1))## [1] TRUEdmat9 <- DelayedTensor::cs_unfold(darr1, m=2)
dmat9_to_darr1 <- DelayedTensor::cs_fold(dmat9, m=2, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat9_to_darr1))## [1] TRUEdmat10 <- DelayedTensor::cs_unfold(darr1, m=3)
dmat10_to_darr1 <- DelayedTensor::cs_fold(dmat10, m=3, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat10_to_darr1))## [1] TRUEIn matvec, m=2 is specified as unfold.
unmatvec is the inverse operation of matvec.
dmat11 <- DelayedTensor::matvec(darr1)
dmat11_darr1 <- DelayedTensor::unmatvec(dmat11, modes=dim(darr1))
identical(as.array(darr1), as.array(dmat11_darr1))## [1] TRUEttm multiplies a tensor by a matrix.
m specifies in which mode the matrix will be multiplied.
dmatZ <- RandomUnifArray(c(10,4))
DelayedTensor::ttm(darr1, dmatZ, m=3)## <2 x 3 x 10> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 1.0149876 0.7891889 0.9992852
## [2,] 1.2213071 1.0469729 1.2977946
##
## ,,2
## [,1] [,2] [,3]
## [1,] 0.8374310 0.5572095 0.8407914
## [2,] 0.6759429 0.5841716 0.9790184
##
## ,,3
## [,1] [,2] [,3]
## [1,] 1.879543 1.742376 1.806140
## [2,] 2.090041 1.906731 2.557365
##
## ...
##
## ,,8
## [,1] [,2] [,3]
## [1,] 1.0774054 1.1412568 0.9926901
## [2,] 1.0360721 1.2516714 1.5304638
##
## ,,9
## [,1] [,2] [,3]
## [1,] 1.797303 1.635262 1.725379
## [2,] 1.719593 1.692060 2.395502
##
## ,,10
## [,1] [,2] [,3]
## [1,] 0.9896005 1.0710785 0.9202772
## [2,] 1.0682445 1.1023058 1.4353843ttl multiplies a tensor by multiple matrices.
ms specifies in which mode these matrices will be multiplied.
dmatX <- RandomUnifArray(c(10,2))
dmatY <- RandomUnifArray(c(10,3))
dlizt <- list(dmatX = dmatX, dmatY = dmatY)
DelayedTensor::ttl(darr1, dlizt, ms=c(1,2))## <10 x 10 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] 1.4334490 1.6424464 1.3071755 . 0.9305834 1.5458981
## [2,] 0.8677958 1.0182677 0.8434610 . 0.6166978 0.9519311
## ... . . . . . .
## [9,] 2.1645734 2.4991337 2.0151627 . 1.4474592 2.3470950
## [10,] 0.6971417 0.8073639 0.6543977 . 0.4716842 0.7575832
##
## ...
##
## ,,4
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] 1.408210 1.819961 1.521677 . 0.8359982 1.4969348
## [2,] 1.079933 1.439732 1.222937 . 0.6379620 1.1473655
## ... . . . . . .
## [9,] 2.3065625 3.0158583 2.5367535 . 1.3668185 2.4514057
## [10,] 0.7663364 1.0061826 0.8481411 . 0.4538144 0.81440132.2 Vectorization
vec collapses a DelayedArray into
a 1D DelayedArray (vector).
Figure 2: Vectorization
DelayedTensor::vec(darr1)## <24> array of class HDF5Array and type "double":
## [1] [2] [3] . [23] [24]
## 0.9403318 0.5682368 0.3605710 . 0.7391881 0.91453312.3 Norm Operations
fnorm calculates the Frobenius norm of a DelayedArray.
Figure 3: Norm Operations
DelayedTensor::fnorm(darr1)## [1] 3.541311innerProd calculates the inner product value of two
DelayedArray.
DelayedTensor::innerProd(darr1, darr2)## [1] 8.8041012.4 Outer Product
Inner product multiplies two tensors and collapses to 0D tensor (norm). On the other hand, the outer product is an operation that leaves all subscripts intact.
Figure 4: Outer Product
DelayedTensor::outerProd(darr1[,,1], darr2[,,1])## <2 x 3 x 2 x 3> array of class HDF5Array and type "double":
## ,,1,1
## [,1] [,2] [,3]
## [1,] 0.5283294 0.2025884 0.5538276
## [2,] 0.3192663 0.2476492 0.5152093
##
## ,,2,1
## [,1] [,2] [,3]
## [1,] 0.8012317 0.3072329 0.8399006
## [2,] 0.4841794 0.3755694 0.7813345
##
## ,,1,2
## [,1] [,2] [,3]
## [1,] 0.7846660 0.3008808 0.8225354
## [2,] 0.4741689 0.3678044 0.7651802
##
## ,,2,2
## [,1] [,2] [,3]
## [1,] 0.06217333 0.02384042 0.06517393
## [2,] 0.03757097 0.02914313 0.06062937
##
## ,,1,3
## [,1] [,2] [,3]
## [1,] 0.2863585 0.1098044 0.3001786
## [2,] 0.1730447 0.1342277 0.2792473
##
## ,,2,3
## [,1] [,2] [,3]
## [1,] 0.3756170 0.1440306 0.3937449
## [2,] 0.2269830 0.1760667 0.36628922.5 Diagonal Operations
Using DelayedDiagonalArray, we can originally create a diagonal
DelayedArray by specifying the dimensions (modes) and the values.
Figure 5: Diagonal Operations
dgdarr <- DelayedTensor::DelayedDiagonalArray(c(5,6,7), 1:5)
dgdarr## <5 x 6 x 7> sparse array of class DelayedArray and type "integer":
## ,,1
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
##
## ...
##
## ,,7
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0Similar to the diag of the base package,
the diag of DelayedTensor is used to extract
and assign values to DelayedArray.
DelayedTensor::diag(dgdarr)## <5> array of class DelayedArray and type "integer":
## [1] [2] [3] [4] [5]
## 1 2 3 4 5DelayedTensor::diag(dgdarr) <- c(1111, 2222, 3333, 4444, 5555)
DelayedTensor::diag(dgdarr)## <5> array of class DelayedArray and type "double":
## [1] [2] [3] [4] [5]
## 1111 2222 3333 4444 55552.6 Mode-wise Operations
modeSum calculates the summation for a given mode m of
a DelayedArray.
The mode specified as m is collapsed into 1D as follows.
Figure 6: Mode-wise Operations
DelayedTensor::modeSum(darr1, m=1)## <1 x 3 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 1.5085685 0.8013422 1.9026943
##
## ,,2
## [,1] [,2] [,3]
## [1,] 1.490210 1.118594 1.656689
##
## ,,3
## [,1] [,2] [,3]
## [1,] 0.6036933 1.6508084 0.8271889
##
## ,,4
## [,1] [,2] [,3]
## [1,] 1.723019 1.395421 1.653721DelayedTensor::modeSum(darr1, m=2)## <2 x 1 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1]
## [1,] 2.286617
## [2,] 1.925988
##
## ,,2
## [,1]
## [1,] 2.045634
## [2,] 2.219858
##
## ,,3
## [,1]
## [1,] 1.209873
## [2,] 1.871817
##
## ,,4
## [,1]
## [1,] 1.991217
## [2,] 2.780944DelayedTensor::modeSum(darr1, m=3)## <2 x 3 x 1> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 2.626130 2.383559 2.523652
## [2,] 2.699361 2.582606 3.516641Similar to modeSum, modeMean calculates the average value
for a given mode m of a DelayedArray.
DelayedTensor::modeMean(darr1, m=1)## <1 x 3 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 0.7542843 0.4006711 0.9513471
##
## ,,2
## [,1] [,2] [,3]
## [1,] 0.7451048 0.5592968 0.8283444
##
## ,,3
## [,1] [,2] [,3]
## [1,] 0.3018466 0.8254042 0.4135944
##
## ,,4
## [,1] [,2] [,3]
## [1,] 0.8615096 0.6977103 0.8268606DelayedTensor::modeMean(darr1, m=2)## <2 x 1 x 4> array of class DelayedArray and type "double":
## ,,1
## [,1]
## [1,] 0.7622056
## [2,] 0.6419961
##
## ,,2
## [,1]
## [1,] 0.6818779
## [2,] 0.7399528
##
## ,,3
## [,1]
## [1,] 0.4032910
## [2,] 0.6239392
##
## ,,4
## [,1]
## [1,] 0.6637390
## [2,] 0.9269813DelayedTensor::modeMean(darr1, m=3)## <2 x 3 x 1> array of class DelayedArray and type "double":
## ,,1
## [,1] [,2] [,3]
## [1,] 0.6565324 0.5958897 0.6309130
## [2,] 0.6748403 0.6456515 0.87916032.7 Tensor Product Operations
There are some tensor specific product such as Hadamard product, Kronecker product, and Khatri-Rao product.
2.7.1 Hadamard Product
Suppose a tensor \(A \in \Re ^{I \times J}\) and a tensor \(B \in \Re ^{I \times J}\).
Hadamard product is defined as the element-wise product of \(A\) and \(B\).
Figure 7: Hadamard Product
Hadamard product can be extended to higher-order tensors.
\[ A \circ B = \begin{bmatrix} a_{11}b_{11} & a_{12}b_{12} & \cdots & a_{1J}b_{1J} \\ a_{21}b_{21} & a_{22}b_{22} & \cdots & a_{2J}b_{2J} \\ \vdots & \vdots & \ddots & \vdots \\ a_{I1}b_{I1} & a_{I2}b_{I2} & \cdots & a_{IJ}b_{IJ} \\ \end{bmatrix} \]
hadamard calculates Hadamard product of two DelayedArray
objects.
prod_h <- DelayedTensor::hadamard(darr1, darr2)
dim(prod_h)## [1] 2 3 4hadamard_list calculates Hadamard product of multiple
DelayedArray objects.
prod_hl <- DelayedTensor::hadamard_list(list(darr1, darr2))
dim(prod_hl)## [1] 2 3 42.7.2 Kronecker Product
Suppose a tensor \(A \in \Re ^{I \times J}\) and a tensor \(B \in \Re ^{K \times L}\).
Kronecker product is defined as all the possible combination of element-wise product and the dimensions of output tensor are \({IK \times JL}\).
Figure 8: Kronecker Product
Kronecker product can be extended to higher-order tensors.
\[ A \otimes B = \begin{bmatrix} a_{11}B & a_{12}B & \cdots & a_{1J}B \\ a_{21}B & a_{22}B & \cdots & a_{2J}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{I1}B & a_{I2}B & \cdots & a_{IJ}B \\ \end{bmatrix} \]
kronecker calculates Kronecker product of two DelayedArray
objects.
prod_kron <- DelayedTensor::kronecker(darr1, darr2)
dim(prod_kron)## [1] 4 9 16kronecker_list calculates Kronecker product of multiple
DelayedArray objects.
prod_kronl <- DelayedTensor::kronecker_list(list(darr1, darr2))
dim(prod_kronl)## [1] 4 9 162.7.3 Khatri-Rao Product
Suppose a tensor \(A \in \Re ^{I \times J}\) and a tensor \(B \in \Re ^{K \times J}\).
Khatri-Rao product is defined as the column-wise Kronecker product and the dimensions of output tensor is \({IK \times J}\).
\[ A \odot B = \begin{bmatrix} a_{1} \otimes a_{1} & a_{2} \otimes a_{2} & \cdots & a_{J} \otimes a_{J} \\ \end{bmatrix} \]
Figure 9: Khatri-Rao Product
Khatri-Rao product can only be used for 2D tensors (matrices).
khatri_rao calculates Khatri-Rao product of two DelayedArray
objects.
prod_kr <- DelayedTensor::khatri_rao(darr1[,,1], darr2[,,1])
dim(prod_kr)## [1] 4 3khatri_rao_list calculates Khatri-Rao product of multiple
DelayedArray objects.
prod_krl <- DelayedTensor::khatri_rao_list(list(darr1[,,1], darr2[,,1]))
dim(prod_krl)## [1] 4 32.8 Utilities Functions
list_rep replicates an arbitrary number of any R object.
str(DelayedTensor::list_rep(darr1, 3))## List of 3
## $ :Formal class 'RandomUnifArray' [package "DelayedRandomArray"] with 1 slot
## .. ..@ seed:Formal class 'RandomUnifArraySeed' [package "DelayedRandomArray"] with 6 slots
## .. .. .. ..@ min : num 0
## .. .. .. ..@ max : num 1
## .. .. .. ..@ dim : int [1:3] 2 3 4
## .. .. .. ..@ chunkdim: int [1:3] 2 3 4
## .. .. .. ..@ seeds :List of 1
## .. .. .. .. ..$ : int [1:2] -179811032 1063580519
## .. .. .. ..@ sparse : logi FALSE
## $ :Formal class 'RandomUnifArray' [package "DelayedRandomArray"] with 1 slot
## .. ..@ seed:Formal class 'RandomUnifArraySeed' [package "DelayedRandomArray"] with 6 slots
## .. .. .. ..@ min : num 0
## .. .. .. ..@ max : num 1
## .. .. .. ..@ dim : int [1:3] 2 3 4
## .. .. .. ..@ chunkdim: int [1:3] 2 3 4
## .. .. .. ..@ seeds :List of 1
## .. .. .. .. ..$ : int [1:2] -179811032 1063580519
## .. .. .. ..@ sparse : logi FALSE
## $ :Formal class 'RandomUnifArray' [package "DelayedRandomArray"] with 1 slot
## .. ..@ seed:Formal class 'RandomUnifArraySeed' [package "DelayedRandomArray"] with 6 slots
## .. .. .. ..@ min : num 0
## .. .. .. ..@ max : num 1
## .. .. .. ..@ dim : int [1:3] 2 3 4
## .. .. .. ..@ chunkdim: int [1:3] 2 3 4
## .. .. .. ..@ seeds :List of 1
## .. .. .. .. ..$ : int [1:2] -179811032 1063580519
## .. .. .. ..@ sparse : logi FALSE2.8.1 Bind Operations
modebind_list collapses multiple DelayedArray objects
into single DelayedArray object.
m specifies the collapsed dimension.
Figure 10: Bind Operations
dim(DelayedTensor::modebind_list(list(darr1, darr2), m=1))## [1] 4 3 4dim(DelayedTensor::modebind_list(list(darr1, darr2), m=2))## [1] 2 6 4dim(DelayedTensor::modebind_list(list(darr1, darr2), m=3))## [1] 2 3 8rbind_list is the row-wise modebind_list and
collapses multiple 2D DelayedArray objects
into single DelayedArray object.
dim(DelayedTensor::rbind_list(list(darr1[,,1], darr2[,,1])))## [1] 4 3cbind_list is the column-wise modebind_list and
collapses multiple 2D DelayedArray objects
into single DelayedArray object.
dim(DelayedTensor::cbind_list(list(darr1[,,1], darr2[,,1])))## [1] 2 6Session information
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DelayedRandomArray_1.2.0 HDF5Array_1.22.0 rhdf5_2.38.0
## [4] DelayedArray_0.20.0 IRanges_2.28.0 S4Vectors_0.32.0
## [7] MatrixGenerics_1.6.0 matrixStats_0.61.0 BiocGenerics_0.40.0
## [10] Matrix_1.3-4 DelayedTensor_1.0.0 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 rTensor_1.4.8 bslib_0.3.1
## [4] compiler_4.1.1 BiocManager_1.30.16 jquerylib_0.1.4
## [7] rhdf5filters_1.6.0 tools_4.1.1 digest_0.6.28
## [10] jsonlite_1.7.2 evaluate_0.14 lattice_0.20-45
## [13] rlang_0.4.12 parallel_4.1.1 yaml_2.2.1
## [16] xfun_0.27 fastmap_1.1.0 stringr_1.4.0
## [19] knitr_1.36 sass_0.4.0 grid_4.1.1
## [22] R6_2.5.1 BiocParallel_1.28.0 rmarkdown_2.11
## [25] bookdown_0.24 irlba_2.3.3 Rhdf5lib_1.16.0
## [28] magrittr_2.0.1 BiocSingular_1.10.0 htmltools_0.5.2
## [31] rsvd_1.0.5 beachmat_2.10.0 dqrng_0.3.0
## [34] ScaledMatrix_1.2.0 stringi_1.7.5 einsum_0.1.0