Matrix operations, algebra and Statistical Methods for Big Data
2025-11-28
Abstract
Description of functions to perform matrix operations, algebra and basic statistical analyses using HDF5 data files.
1 Overview
This package implements several matrix operations using Matrix objects as well as HDF5 data files. Some basic algebra operations that can also be computed that are useful to implement statistical analyses using standard methodologies such as principal component analyses (PCA) or least squares estimation. The package also contains specific statistical methods mainly used in omic data analysis such as lasso regression.
2 Prerequisites
The package requires other packages from CRAN and Bioconductor to be installed.
CRAN:Matrix,RcppEigenandRSpectra.Bioconductor:HDF5array,rhdf5
The user can execute this code to install the required packages:
# Install BiocManager (if not previously installed)
install.packages("BiocManager")
# Install required packages
BiocManager::install(c("Matrix", "RcppEigen", "RSpectra",
"HDF5Array", "rhdf5"))Our package needs to be installed from source code. In such cases, a collection of software (e.g. C, C++, Fortran, …) are required, mainly for Windows users. These programs can be installed using Rtools.
3 Install package
Once required packages and Rtools are installed, BigDataStatMeth package can be installed from our GitHub repository as follows:
# Install devtools and load library (if not previously installed)
install.packages("devtools")
library(devtools)
# Install BigDataStatMeth
install_github("isglobal-brge/BigDataStatMeth")4 Getting started
First, let us start by loading the required packages to describe the main capabilities of the package
library(rhdf5)
library(BigDataStatMeth)5 Previous knowledge
5.1 HDF5 data files
Hierarchical Data Format (HDF) is a set of file formats designed to store and organize large amounts of data. It is supported by The HDF Group, a non-profit corporation whose mission is to ensure continued development of HDF5 technologies and the continued accessibility of data stored in HDF.
HDF5 is a technology suite that makes possible the management of extremely large and complex data collections, can accommodate virtually every kind of data in a single file, sequences, images, SNP matrices, and every other type of data and metadata associated with an experiment.
There is no limit on the number or size of data objects in the collection, giving great flexibility for omic data. Is high-performance I/O with a rich set of integrated performance features that allow for access time and storage space optimizations
HDF5 file structure include only two major types of object:
Datasets, which are multidimensional arrays of a homogeneous type. For example, datasets for omics data could be to genomics, transcriptomics, epigenomics, proteomics and/or metabolomics experimentsGroups, which are container structures which can hold datasets and other groups
This results in a truly hierarchical, filesystem-like data format
Figure 1: HDF5 hierarchical structure
5.2 Basics in hdf5 files using R
Create hdf5 file
We have implemented the bdCreate_hdf5_matrix_file () function to create an hdf5 file with a group and a dataset in a single step. This function allows to create datasets from a standard R matrices from Matrix objects.
library(rhdf5)
set.seed(5234)
n <- 500
m <- 600
A <- matrix(rnorm(n*m,mean=0,sd=1), n,m)
# We also can create a dataset from R matrix object
bdCreate_hdf5_matrix(filename = "robject.hdf5",
object = A,
group = "INPUT",
dataset = "A",
overwriteFile = TRUE)
$fn
[1] "robject.hdf5"
$ds
[1] "INPUT/A"We see 0 in the console, indicating that no errors where found when creating the hdf5 file. Notice that a file called “robject.hdf5” will be created in the working directory
list.files(pattern = "*.hdf5")
[1] "robject.hdf5"Add datasets in hdf5 file
The function bdCreate_hdf5_matrix() also allows to add a dataset in a existing hdf5 file. We can create the dataset in any group, if group doesn’t exists in file, the group is created before append the dataset.
set.seed(5234)
n <- 500
m <- 1000
A <- matrix(rnorm(n*m,mean=0,sd=1), n, m)
set.seed(5234)
n <- 1000
m <- 12000
B <- matrix(rnorm(n*m,mean=3,sd=0.5), n, m)
# Path to HDF5 file
example_fn <- "delayed.hdf5"
# We create another data file (delayed.hdf5) with a matrix A.
# The group is called INPUT, overwriteFile is set to true to
# overwrite a file if exists
bdCreate_hdf5_matrix(filename = example_fn,
object = A,
group = "INPUT",
dataset = "A",
overwriteFile = TRUE)
$fn
[1] "delayed.hdf5"
$ds
[1] "INPUT/A"
# And them, we add another matrix B to the same group
bdCreate_hdf5_matrix(object = B,
filename = example_fn,
group = "INPUT",
dataset = "B")
$fn
[1] "delayed.hdf5"
$ds
[1] "INPUT/B"Open and get hdf5 content file
We can open an existing file show contents and access data using functions from rhdf5 package. rhdf5 is an R interface for HDF5. The file must always be opened before working with it.
With h5ls() function we get the the list of an hdf5 content without the need of open the file.
# Examine hierarchy before open file
h5ls(example_fn)
group name otype dclass dim
0 / INPUT H5I_GROUP
1 /INPUT A H5I_DATASET FLOAT 500 x 1000
2 /INPUT B H5I_DATASET FLOAT 1000 x 12000Once opened, it can be seen just typing its name. In that case, we only get the current level in the hierarchical tree
# Open file
h5fdelay <- H5Fopen(example_fn)
# Show hdf5 hierarchy (groups)
h5fdelay
HDF5 FILE
name /
filename
name otype dclass dim
0 INPUT H5I_GROUP NOTE: We can also use hdfview
Access datasets data
The $ operator can be used to access the next group level, this operator reads the object from disk. We can assign the dataset contents to an R object in order to work with it.
bdata <- h5fdelay$INPUT$B
bdata[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 3.562285 2.343400 2.952424 2.871667 3.569449
[2,] 2.702681 2.455068 3.006128 3.371495 3.256128
[3,] 3.057644 2.993982 3.833475 3.758272 2.448493And finally, we close the opened files.
h5closeAll()6 Convert text file to HDF5
Usually we have files with a large amount of data to be loaded into R. For Big data problem, R can’t deal with this data so we cannot efficiently work with then. In order to resolve this issue, you can work directly from data on disk with datasets stored in an HDF5 data file.
You can easily change your text files to HDF5 using the function bdImportTextFile_hdf5 (), to do that, we only have to define the input file with the filename parameter and the HDF5 destination dataset with parameters outputfile for the output file name, outGroup with the group to store the dataset and the dataset name with the parameter outDataset. It should be considered that:
This function only allows to import numeric data except for rownames and colnames where character are allowed.
Data can be delimited by different characters. By default, we use tabs (“\t”) but it can be changed easily with
sepparameter.
In this example we convert the “colesterol.csv” file and we will save the data to a file ‘colesterol_file.hdf5’ under the group ‘COLESTEROLGROUP’ in the dataset ‘COLESTEROLDATASET’. In this example the text file contains column names so we set the parameter header = TRUE.
import_hdf5 <- bdImportTextFile_hdf5(filename = "colesterol.csv",
sep=',',
outputfile = "colesterol_file.hdf5",
outGroup = "COLESTEROL",
outDataset = "COLESTEROLDATA",
header = TRUE,
overwrite = TRUE)
If we observe the new file colesterol_file.hdf5 and its content
# Show content
h5ls(import_hdf5$fn)
group name otype
0 / COLESTEROL H5I_GROUP
1 /COLESTEROL .COLESTEROLDATA_dimnames H5I_GROUP
2 /COLESTEROL/.COLESTEROLDATA_dimnames 2 H5I_DATASET
3 /COLESTEROL COLESTEROLDATA H5I_DATASET
dclass dim
0
1
2 COMPOUND 10
3 FLOAT 1000 x 10
# We can open the file and have access to the data
res_hdf5 <- h5read(import_hdf5$fn, import_hdf5$ds)
colnames_hdf5 <- h5read(import_hdf5$fn, import_hdf5$ds_cols)[,1]
# Show hdf5 content dataset
res_hdf5[1:5, 1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 223.7348 55.25039 14.90246 0.9026095 10.22067 1.450970
[2,] 248.1820 53.47404 25.12592 0.8710345 17.10225 1.002655
[3,] 180.2071 58.54268 12.95197 0.8882435 11.21530 1.073924
[4,] 200.1522 62.58284 28.27819 0.9361357 12.79399 1.129373
[5,] 234.3281 68.70254 13.21404 0.7775238 11.39689 1.235655
# Show colnames
head(colnames_hdf5)
[1] "TCholesterol" "Age" "Insulin" "Creatinine" "BUN"
[6] "LLDR"
# Show data with colnames
colnames(res_hdf5) <- colnames_hdf5
res_hdf5[1:5, 1:6]
TCholesterol Age Insulin Creatinine BUN LLDR
[1,] 223.7348 55.25039 14.90246 0.9026095 10.22067 1.450970
[2,] 248.1820 53.47404 25.12592 0.8710345 17.10225 1.002655
[3,] 180.2071 58.54268 12.95197 0.8882435 11.21530 1.073924
[4,] 200.1522 62.58284 28.27819 0.9361357 12.79399 1.129373
[5,] 234.3281 68.70254 13.21404 0.7775238 11.39689 1.235655NOTE: We can overwrite an existing hdf5 dataset with new data by setting overwrite = TRUE and file wit overwriteFile = TRUE.
7 Operation with HDF5 files
Once our data are available as hdf5 files, we can operate directly on it without having data on memory (i.e we have access from disk. In BigDataStatMeth we have implemented most of the common matrix operations and algebra functions to help users how are not familiar with this type of implementation. Next section provides several examples
7.1 HDF5 Matrix Multiplication
If we have the matrices stored in hdf5 data file, we can use the function bdblockmult_hdf5(). This function allows to perform a matrix multiplication from two matrices stored in a file.
The bdblockmult_hdf5() function in the BigDataStatMeth package performs efficient matrix multiplication by blocks, directly on datasets stored in an HDF5 file. This allows scalable computation even when the matrices do not fit into memory.
If the outgroup is not specified, the result is written to the "OUTPUT" group by default. Likewise, if the outdataset is not provided, the output dataset name is constructed as "A_x_B", where A and B are the names of the input datasets.
# Perform blockwise matrix multiplication
res <- bdblockmult_hdf5(filename = example_fn, group = "INPUT",
A = "A", B = "B", outgroup = "HDF5_RES")
# Show the content of the HDF5 file
h5ls(res$fn)
group name otype dclass dim
0 / HDF5_RES H5I_GROUP
1 /HDF5_RES A_x_B H5I_DATASET FLOAT 500 x 12000
2 / INPUT H5I_GROUP
3 /INPUT A H5I_DATASET FLOAT 500 x 1000
4 /INPUT B H5I_DATASET FLOAT 1000 x 12000We can then extract a portion of the result to validate correctness:
# Extract the result from HDF5
result_hdf5 <- h5read(res$fn, res$ds)[1:3, 1:5]
result_hdf5
[,1] [,2] [,3] [,4] [,5]
[1,] 146.46071 168.18088 185.33397 213.72110 198.21560
[2,] 165.83147 192.13943 197.74114 181.76050 175.02919
[3,] -29.83376 -48.24969 -67.16883 -53.60578 -46.85563
# Compute the same multiplication in R
result_r <- (A %*% B)[1:3, 1:5]
result_r
[,1] [,2] [,3] [,4] [,5]
[1,] 146.46071 168.18088 185.33397 213.72110 198.21560
[2,] 165.83147 192.13943 197.74114 181.76050 175.02919
[3,] -29.83376 -48.24969 -67.16883 -53.60578 -46.85563
# Compare both results
all.equal((A %*% B), h5read(res$fn, res$ds))
[1] TRUEThis confirms that the result obtained with the C++ HDF5-based backend is numerically identical to the result obtained using R’s in-memory multiplication.
7.2 Crossproduct and Transposed Crossproduct in HDF5
The functions bdCrossprod_hdf5() and bdtCrossprod_hdf5() from the BigDataStatMeth package implement efficient versions of the crossproduct operations using blockwise computation on HDF5 datasets. These functions allow working with large matrices that may not fit in memory, by performing the computation in chunks.
Both functions can operate on a single matrix or on two matrices:
- If only matrix
Ais provided:bdCrossprod_hdf5()computest(A) %*% AbdtCrossprod_hdf5()computesA %*% t(A)
- If both matrices
AandBare provided:bdCrossprod_hdf5()computest(A) %*% BbdtCrossprod_hdf5()computesA %*% t(B)
By default, the output is stored in the "OUTPUT" group. If outdataset is not specified, the default dataset names are:
- "Crossprod_A_x_B" for bdCrossprod_hdf5()
- "tCrossprod_A_x_B" for bdtCrossprod_hdf5()
where A and B refer to the names of the input datasets provided via the A and B arguments, respectively. If only matrix A is specified (i.e., B is missing), then the name becomes “Crossprod_A” or “tCrossprod_A”. This naming convention ensures traceability of the operation and the involved datasets within the HDF5 file structure.
Crossproduct of a single matrix
# Create example matrices
set.seed(123)
A <- matrix(rnorm(1000 * 200), nrow = 1000, ncol = 200)
B <- matrix(rnorm(1000 * 150), nrow = 1000, ncol = 150)
C <- matrix(rnorm(800 * 200), nrow = 800, ncol = 200)
# Save matrices to HDF5 file using BigDataStatMeth
example_fn <- "delayed.hdf5"
bdCreate_hdf5_matrix(filename = example_fn,
object = A,
group = "INPUT",
dataset = "A",
overwriteFile = TRUE)
$fn
[1] "delayed.hdf5"
$ds
[1] "INPUT/A"
bdCreate_hdf5_matrix(filename = example_fn,
object = B,
group = "INPUT",
dataset = "B",
overwriteFile = FALSE)
$fn
[1] "delayed.hdf5"
$ds
[1] "INPUT/B"
bdCreate_hdf5_matrix(filename = example_fn,
object = C,
group = "INPUT",
dataset = "C",
overwriteFile = FALSE)
$fn
[1] "delayed.hdf5"
$ds
[1] "INPUT/C"
# Compute t(A) %*% A
res_cross <- bdCrossprod_hdf5(filename = example_fn,
group = "INPUT",
A = "A")
# Show where the result is stored
h5ls(res_cross$fn)
group name otype dclass dim
0 / INPUT H5I_GROUP
1 /INPUT A H5I_DATASET FLOAT 1000 x 200
2 /INPUT B H5I_DATASET FLOAT 1000 x 150
3 /INPUT C H5I_DATASET FLOAT 800 x 200
4 / OUTPUT H5I_GROUP
5 /OUTPUT CrossProd_A_x_A H5I_DATASET FLOAT 200 x 200
# Compare with R's crossprod
res_hdf5 <- h5read(res_cross$fn, res_cross$ds)
res_r <- crossprod(A)
all.equal(res_r, res_hdf5)
[1] TRUE- For
bdtCrossprod_hdf5(), it is"tcrossprod_A".
Crossproduct of two matrices
# Compute t(A) %*% B
res_cross2 <- bdCrossprod_hdf5(filename = example_fn,
group = "INPUT",
A = "A",
B = "B")
# Compare with R
res_hdf5 <- h5read(res_cross2$fn, res_cross2$ds)
res_r <- crossprod(A, B)
all.equal(res_r, res_hdf5)
[1] TRUETransposed Crossproduct of a single matrix
# Compute A %*% t(A)
res_tcross <- bdtCrossprod_hdf5(filename = example_fn,
group = "INPUT",
A = "A")
h5ls(res_tcross$fn)
group name otype dclass dim
0 / INPUT H5I_GROUP
1 /INPUT A H5I_DATASET FLOAT 1000 x 200
2 /INPUT B H5I_DATASET FLOAT 1000 x 150
3 /INPUT C H5I_DATASET FLOAT 800 x 200
4 / OUTPUT H5I_GROUP
5 /OUTPUT CrossProd_A_x_A H5I_DATASET FLOAT 200 x 200
6 /OUTPUT CrossProd_A_x_B H5I_DATASET FLOAT 200 x 150
7 /OUTPUT tCrossProd_A_x_A H5I_DATASET FLOAT 1000 x 1000
res_hdf5_t <- h5read(res_tcross$fn, res_tcross$ds)
res_r_t <- tcrossprod(A)
all.equal(res_r_t, res_hdf5_t)
[1] TRUETransposed Crossproduct of two matrices
# Compute A %*% t(B)
res_tcross2 <- bdtCrossprod_hdf5(filename = example_fn,
group = "INPUT",
A = "A",
B = "C")
res_hdf5_t <- h5read(res_tcross2$fn, res_tcross2$ds)
res_r_t <- tcrossprod(A, C)
all.equal(res_r_t, res_hdf5_t)
[1] TRUEThese examples confirm that both crossproduct and transposed crossproduct computations done with the HDF5-based C++ backend are equivalent to the results produced by R’s crossprod() and tcrossprod() functions, while allowing better scalability and performance for large datasets.
Crossproduct Combinations Supported
The following table summarizes the types of crossproduct operations supported by the functions bdCrossprod_hdf5() and bdtCrossprod_hdf5() in the BigDataStatMeth package.
| Function | Matrix A | Matrix B (optional) | Operation |
|---|---|---|---|
bdCrossprod_hdf5() |
A | — | t(A) %*% A |
bdCrossprod_hdf5() |
A | B | t(A) %*% B |
bdtCrossprod_hdf5() |
A | — | A %*% t(A) |
bdtCrossprod_hdf5() |
A | B | A %*% t(B) |
These operations are executed efficiently using C++ backends and blockwise computation on HDF5 datasets, making them suitable for large-scale matrix operations.
7.3 Blockwise Matrix Subtraction in HDF5
The bdblockSubstract_hdf5() function performs blockwise element-wise subtraction of two matrices stored in an HDF5 file. It computes \(A - B\), writing the result back into the file, without requiring either matrix to be fully loaded into memory.
# Create two matrices of the same dimensions
set.seed(42)
A_sub <- matrix(rnorm(1000 * 300), nrow = 1000, ncol = 300)
B_sub <- matrix(rnorm(1000 * 300), nrow = 1000, ncol = 300)
# Save them to HDF5
fn_sub <- "subtraction_example.hdf5"
bdCreate_hdf5_matrix(filename = fn_sub, object = A_sub,
group = "INPUT", dataset = "A_sub",
overwriteFile = TRUE)
$fn
[1] "subtraction_example.hdf5"
$ds
[1] "INPUT/A_sub"
bdCreate_hdf5_matrix(filename = fn_sub, object = B_sub,
group = "INPUT", dataset = "B_sub",
overwriteFile = FALSE)
$fn
[1] "subtraction_example.hdf5"
$ds
[1] "INPUT/B_sub"
# Perform subtraction: A - B
res_sub <- bdblockSubstract_hdf5(filename = fn_sub,
group = "INPUT",
A = "A_sub", B = "B_sub")
# Compare a subset with R
result_hdf5 <- h5read(res_sub$fn, res_sub$ds)
result_r <- A_sub - B_sub
all.equal(result_r, result_hdf5)
[1] TRUEThe result stored in HDF5 is numerically identical to the in-memory computation of \(A - B\) in R, validating the blockwise approach.
7.4 Blockwise Matrix Addition in HDF5
The bdblockSum_hdf5() function performs blockwise element-wise addition of two matrices stored in HDF5 format. This is useful for efficiently computing \(A + B\) on large matrices.
# Create two compatible matrices
set.seed(99)
A_add <- matrix(rnorm(800 * 250), nrow = 800, ncol = 250)
B_add <- matrix(rnorm(800 * 250), nrow = 800, ncol = 250)
# Save them to HDF5
fn_add <- "addition_example.hdf5"
bdCreate_hdf5_matrix(filename = fn_add, object = A_add,
group = "INPUT", dataset = "A",
overwriteFile = TRUE)
$fn
[1] "addition_example.hdf5"
$ds
[1] "INPUT/A"
bdCreate_hdf5_matrix(filename = fn_add, object = B_add,
group = "INPUT", dataset = "B",
overwriteFile = FALSE)
$fn
[1] "addition_example.hdf5"
$ds
[1] "INPUT/B"
# Perform addition: A + B
res_add <- bdblockSum_hdf5(filename = fn_add,
group = "INPUT",
A = "A", B = "B")
# Compare result with R
result_hdf5 <- h5read(res_add$fn, res_add$ds)
result_r <- A_add + B_add
all.equal(result_r, result_hdf5)
[1] TRUELike subtraction, the addition operation yields identical results to standard R arithmetic but can scale to matrices stored on disk using the HDF5 format.
8 Matrix Decompositions and Inverses
8.1 Singular Value Decomposition (SVD)
The SVD of an \(m \times n\) real or complex matrix \(A\) is a factorization of the form:
\[U\Sigma { V }^{ T }\]
where :
- \(U\) is a \(m \times m\) real or complex unitary matrix
- \(\Sigma\) is a \(m \times n\) rectangular diagonal matrix with non-negative real numbers on the diagonal
- \(V\) is a \(n \times n\) real or complex unitary matrix.
Notice that:
- The diagonal entries \(\sigma_i\) of \(\Sigma\) are known as the singular values of \(A\).
- The columns of \(U\) are called the left-singular vectors of \(A\).
- The columns of \(V\) are called the right-singular vectors of \(A\).
8.1.1 Block Singular Values Decomposition
This method was developed by M. A. Iwen and B. W. Ong (2016). The authors introduced a distributed and incremental SVD algorithm that is useful for agglomerative data analysis on large networks. The algorithm calculates the singular values and left singular vectors of a matrix A, by first, partitioning it by columns. This creates a set of submatrices of A with the same number of rows, but only some of its columns. After that, the SVD of each of the submatrices is computed. The final step consists of combining the results obtained by merging them again and computing the SVD of the resulting matrix.
Figure 2: Flowchart for a two-level hierarchical Block SVD algorithm
This method is implemented in bdSVD_hdf5 function, this function works directly on hdf5 data format, loading in memory only the data to perform calculations and saving the results again in the hdf5 file for later use.
We have to indicate the file to work with, the dataset name and the group where the dataset is located :
# Create dataframe data with 'odata' matrix in delayed hdf5 file at OMIC group
set.seed(5234)
n <- 100
m <- 15000
omicdata <- matrix(rnorm(n*m, mean=0, sd=1), n,m)
bdCreate_hdf5_matrix(filename = example_fn,
object = omicdata,
group = "OMICS",
dataset = "data",
overwriteDataset = TRUE)
$fn
[1] "delayed.hdf5"
$ds
[1] "OMICS/data"
# Direct from hdf5 data file
svdh5 <- bdSVD_hdf5( filename = example_fn,
group = "OMICS",
dataset = "data",
overwrite = TRUE)
# get results svd (d) from hdf5 data file
svd_hdf5_d <- h5read(svdh5$fn, svdh5$ds_d)
# Results in hdf5 file for d
svd_hdf5_d[1:7]
[1] 131.8042 131.1643 130.6787 130.6179 130.5750 130.3000 130.1092
svd <- svd(scale(omicdata))
svd$d[1:7]
[1] 131.8042 131.1643 130.6787 130.6179 130.5750 130.3000 130.1092Like in Simple Singular Value Decomposition we can normalize, center or scale data before proceed with SVD decomposition with bscale and bcenter parameters, by default this parameter are TRUE, data is normalized before SVD decomposition. To proceed with SVD without normalization :
# Direct from hdf5 data file (using only one thread, serial execution)
svdh5 <- bdSVD_hdf5( filename = example_fn,
group = "OMICS",
dataset = "data",
bcenter = FALSE,
bscale = FALSE,
overwrite = TRUE)# get results svd (d)
svd_hdf5_d <- h5read(svdh5$fn, svdh5$ds_d)[1:7]
# SVD (d) from file - data not normalized
svd_hdf5_d
[1] 132.0820 131.4080 130.9686 130.8555 130.6781 130.4269 130.3248
# with R implementation from data in memory
svd <- svd(omicdata)
svd$d[1:7]
[1] 132.0820 131.4080 130.9686 130.8555 130.6781 130.4269 130.3248In the SVD decomposition by blocks we can indicate the number of decomposition levels and number of local SVDs to concatenate at each level with parameters q and k respectively, by default q = 1 one level with k=2.
# Block decomposition with 1 level and 4 local SVDs at each level using
# two threads (as maximum)
svdh5 <- bdSVD_hdf5( filename = example_fn,
group = "OMICS",
dataset = "data",
q = 1,
k = 4,
threads = 2,
overwrite = TRUE)
# get results svd (d)
svd_hdf5_d <- h5read(svdh5$fn, svdh5$ds_d)[1:7]
# SVD (d) from file - data not normalized
svd_hdf5_d
[1] 131.8042 131.1643 130.6787 130.6179 130.5750 130.3000 130.1092
# with R implementation from data in memory
svd <- svd(scale(omicdata))
svd$d[1:7]
[1] 131.8042 131.1643 130.6787 130.6179 130.5750 130.3000 130.10928.2 Cholesky Decomposition
The Cholesky decomposition is a matrix factorization technique applicable to symmetric, positive-definite matrices. It decomposes a matrix \(A\) into the product of a lower triangular matrix \(L\) and its transpose, such that: \[ A = L \cdot L^\top \] This method is particularly useful for solving systems of linear equations, inverting matrices, and performing efficient numerical computations.
The BigDataStatMeth package provides the function bdInvCholesky_hdf5() to compute the inverse of a symmetric, positive-definite matrix stored in HDF5 format using the Cholesky decomposition. This function operates directly on disk-stored data, minimizing memory usage while maintaining computational efficiency. If the input matrix meets the necessary conditions (Hermitian and positive-definite), the function performs the factorization and inversion in a blockwise manner, making it suitable for large-scale data.
Here’s how to use bdInvCholesky_hdf5():
N <- 100
set.seed(5234)
Y <- matrix(rnorm(N*N), N, N)
Ycp <- crossprod(Y)
bdCreate_hdf5_matrix(filename = example_fn,
object = Ycp,
group = "chol",
dataset = "data",
transp = FALSE,
overwriteFile = TRUE, overwriteDataset = TRUE,
unlimited = FALSE)
$fn
[1] "delayed.hdf5"
$ds
[1] "chol/data"
cholh5 <- bdCholesky_hdf5(filename = example_fn,
group = "chol",
dataset = "data",
outdataset = "matrixDec",
outgroup = "Cholesky_Dec",
overwrite = TRUE)
Warning: 'memory.size()' is Windows-specific
choldesc_hdf5 <- h5read(cholh5$fn, cholh5$ds)
choldesc_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.726702 2.161678 1.3276203 -0.8431624 0.6066569
[2,] 0.000000 10.339840 -0.1125209 -0.3722274 0.2753345
[3,] 0.000000 0.000000 9.9295256 0.9489195 -0.3208095
choldesc_r <- chol(Ycp)
choldesc_r[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.726702 2.161678 1.3276203 -0.8431624 0.6066569
[2,] 0.000000 10.339840 -0.1125209 -0.3722274 0.2753345
[3,] 0.000000 0.000000 9.9295256 0.9489195 -0.3208095
all.equal(choldesc_hdf5, choldesc_r)
[1] TRUEThe computation can be further optimized by specifying the elementsBlock parameter, which determines the approximate number of elements to include in each processing block. For example, setting elementsBlock = 50 instructs the function to process the matrix in chunks of approximately 50 elements, balancing speed and memory use.
cholh5 <- bdCholesky_hdf5(filename = example_fn,
group = "chol",
dataset = "data",
outdataset = "matrixDec_50",
outgroup = "Cholesky_Dec",
elementsBlock = 50,
overwrite = TRUE)
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
# Result dataset
cholh5$ds
[1] "Cholesky_Dec/matrixDec_50"
choldesc_hdf5 <- h5read(cholh5$fn, cholh5$ds)
choldesc_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.726702 2.161678 1.3276203 -0.8431624 0.6066569
[2,] 0.000000 10.339840 -0.1125209 -0.3722274 0.2753345
[3,] 0.000000 0.000000 9.9295256 0.9489195 -0.3208095
choldesc_r <- chol(Ycp)
choldesc_r[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.726702 2.161678 1.3276203 -0.8431624 0.6066569
[2,] 0.000000 10.339840 -0.1125209 -0.3722274 0.2753345
[3,] 0.000000 0.000000 9.9295256 0.9489195 -0.3208095
all.equal(choldesc_hdf5, choldesc_r)
[1] TRUE8.3 QR Decomposition
The QR decomposition is a matrix factorization technique in which a matrix \(A\) is decomposed into the product of an orthogonal matrix \(Q\) and an upper triangular matrix \(R\), such that:
\[ A = Q \cdot R \]
This decomposition is widely used in numerical linear algebra for solving least squares problems, eigenvalue computations, and matrix inversion. In the context of large-scale data, computing the QR decomposition efficiently and in a memory-safe manner becomes essential.
The BigDataStatMeth package implements two functions to compute the inverse of a matrix using QR decomposition:
bdQR()for in-memory matrices.bdQR_hdf5()for matrices stored in HDF5 format.
The computation is done in blocks, allowing it to scale to datasets that do not fit in memory. The result can be used in further downstream statistical procedures, such as CCA or regression.
In this example, we apply the QR decomposition to the same matrix data inside the group chol that was created and stored in HDF5 format in the previous section.
QRh5 <- bdQR_hdf5(filename = example_fn,
group = "chol",
dataset = "data",
outgroup = "QR_Dec",thin = TRUE,
overwrite = TRUE)
# Result dataset
QR_Q_hdf5 <- h5read(QRh5$fn, QRh5$ds_Q)
QR_R_hdf5 <- h5read(QRh5$fn, QRh5$ds_R)
# Q matrix
QR_Q_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -0.68551057 0.04546917 0.01531785 0.0008194793 0.034269378
[2,] -0.15234901 -0.74590341 0.02851012 0.0182487534 0.005718736
[3,] -0.09356695 0.01443328 -0.70676468 0.0584085661 -0.032816534
# R matrix
QR_R_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -138.0121 -40.05205 -22.151218 11.2513535 -14.627475
[2,] 0.0000 -141.41677 -3.081843 0.6397727 -3.862321
[3,] 0.0000 0.00000 -139.144722 -25.7224753 8.435610
# Q matrix in R
QR_Q_r <- qr.Q(qr(Ycp))
QR_Q_r[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -0.68551057 0.04546917 0.01531785 0.0008194793 0.034269378
[2,] -0.15234901 -0.74590341 0.02851012 0.0182487534 0.005718736
[3,] -0.09356695 0.01443328 -0.70676468 0.0584085661 -0.032816534
all.equal(QR_Q_hdf5, QR_Q_r)
[1] TRUEIn this example, we show how to control the block size used during the computation by setting the block_size parameter. Adjusting this value can help improve performance depending on the available memory and system architecture.
QRh5 <- bdQR_hdf5(filename = example_fn,
group = "chol",
dataset = "data",
outgroup = "QR_Dec",
block_size = 256,
overwrite = TRUE)
# Result dataset
QR_Q_hdf5 <- h5read(QRh5$fn, QRh5$ds_Q)
QR_R_hdf5 <- h5read(QRh5$fn, QRh5$ds_R)
# Q matrix
QR_Q_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -0.68551057 0.04546917 0.01531785 0.0008194793 0.034269378
[2,] -0.15234901 -0.74590341 0.02851012 0.0182487534 0.005718736
[3,] -0.09356695 0.01443328 -0.70676468 0.0584085661 -0.032816534
# R matrix
QR_R_hdf5[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -138.0121 -40.05205 -22.151218 11.2513535 -14.627475
[2,] 0.0000 -141.41677 -3.081843 0.6397727 -3.862321
[3,] 0.0000 0.00000 -139.144722 -25.7224753 8.435610
# Q matrix in R
QR_Q_r <- qr.Q(qr(Ycp))
QR_Q_r[1:3,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] -0.68551057 0.04546917 0.01531785 0.0008194793 0.034269378
[2,] -0.15234901 -0.74590341 0.02851012 0.0182487534 0.005718736
[3,] -0.09356695 0.01443328 -0.70676468 0.0584085661 -0.032816534
all.equal(QR_Q_hdf5, QR_Q_r)
[1] TRUE8.4 Cholesky-based Matrix Inversion
The Cholesky decomposition is a matrix factorization technique applicable to symmetric, positive-definite matrices. It expresses a matrix \(A\) as the product of a lower triangular matrix \(L\) and its transpose:
\[ A = L \cdot L^\top \]
This decomposition is particularly useful for solving linear systems and computing matrix inverses in a numerically stable and efficient way.
The BigDataStatMeth package implements two functions to compute the inverse of a matrix using Cholesky decomposition:
bdInvCholesky()for in-memory matrices.bdInvCholesky_hdf5()for matrices stored in HDF5 format.
These functions are especially useful when dealing with large covariance matrices in genomics or other high-dimensional applications.
invCholh5 <- bdInvCholesky_hdf5( filename = example_fn,
group = "chol",
dataset = "data",
outdataset = "invmatrix",
outgroup = "InvCholesky",
fullMatrix = FALSE,
overwrite = TRUE)
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
# Result dataset
invChol_hdf5 <- h5read(invCholh5$fn, invCholh5$ds)
invChol_hdf5[1:5,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 70.55445 -59.13169 5.832314 -20.482406 -57.516370
[2,] 0.00000 50.61402 -4.833963 17.358636 49.337577
[3,] 0.00000 0.00000 0.915268 -1.606462 -4.476971
[4,] 0.00000 0.00000 0.000000 6.222298 16.988061
[5,] 0.00000 0.00000 0.000000 0.000000 48.897665The argument fullMatrix is a logical flag that controls how the inverse is stored:
- If fullMatrix = TRUE, the entire inverse matrix is stored in the output file.
- If fullMatrix = FALSE (default), only the lower triangular part is stored, which reduces disk usage and is sufficient when the matrix is symmetric.
invCholh5 <- bdInvCholesky_hdf5( filename = example_fn,
group = "chol",
dataset = "data",
outdataset = "invmatrix",
outgroup = "InvCholesky",
fullMatrix = TRUE,
overwrite = TRUE)
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
Warning: 'memory.size()' is Windows-specific
# Result dataset
invChol_hdf5 <- h5read(invCholh5$fn, invCholh5$ds)
invChol_hdf5[1:5,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 70.554448 -59.131685 5.832314 -20.482406 -57.516370
[2,] -59.131685 50.614021 -4.833963 17.358636 49.337577
[3,] 5.832314 -4.833963 0.915268 -1.606462 -4.476971
[4,] -20.482406 17.358636 -1.606462 6.222298 16.988061
[5,] -57.516370 49.337577 -4.476971 16.988061 48.897665
# Inverse Cholesky matrix in R
invChol_r <- solve(Ycp)
invChol_r[1:5,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 70.554448 -59.131685 5.832314 -20.482406 -57.516370
[2,] -59.131685 50.614021 -4.833963 17.358636 49.337577
[3,] 5.832314 -4.833963 0.915268 -1.606462 -4.476971
[4,] -20.482406 17.358636 -1.606462 6.222298 16.988061
[5,] -57.516370 49.337577 -4.476971 16.988061 48.897665
all.equal(invChol_hdf5, invChol_r)
[1] TRUE9 Session information
sessionInfo()
R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Madrid
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] BigDataStatMeth_1.0.2 rhdf5_2.52.1 knitr_1.50
[4] BiocStyle_2.36.0
loaded via a namespace (and not attached):
[1] cli_3.6.5 rlang_1.1.6 xfun_0.52
[4] data.table_1.17.8 jsonlite_2.0.0 RCurl_1.98-1.17
[7] htmltools_0.5.8.1 sass_0.4.10 rmarkdown_2.29
[10] evaluate_1.0.4 jquerylib_0.1.4 bitops_1.0-9
[13] fastmap_1.2.0 yaml_2.3.10 lifecycle_1.0.4
[16] Rhdf5lib_1.30.0 bookdown_0.43 BiocManager_1.30.26
[19] compiler_4.5.1 Rcpp_1.1.0 rhdf5filters_1.20.0
[22] rstudioapi_0.17.1 digest_0.6.37 R6_2.6.1
[25] bslib_0.9.0 tools_4.5.1 cachem_1.1.0