hpages/Zarr_vs_HDF5_part1_First_benchmarks.md

Zarr vs HDF5: First benchmarks & impact of chunk size

About chunks

Like HDF5, the Zarr format stores the array data in physical chunks:

Chunks are typically small rectangular regions of the array (e.g. 100x100).
They span the full array with no overlaps (tiles).
They are compressed.

Chunks are the read/write units:

Reading values from a chunk typically requires loading and uncompressing the full chunk in memory.
Some compression algorithms allow random read access so a small set of values can be read without loading and uncompressing the full chunk. BUT random read access from compressed data can be slow so only worth if number of values to read is small.
Writing a single array value to a chunk requires loading and uncompressing the full chunk in memory + updating it + recompressing and writing back the full chunk to disk.

Choosing a good chunk geometry is always a trade off:

Small chunks increase granularity so reading a set of elements that are "spread out" tends to generate less I/O when the chunks are small.
However compression ratio degrades when chunks become too small -> array takes more space on disk.
Two types of overheads to pay for each chunk:
1. disk space overhead
2. time overhead (e.g. some book keeping) These are fixed overheads i.e. they don't depend on the size of the chunk.
Because small chunks means more chunks --> small chunks also means more total overhead (space & time).

BIG DIFFERENCE BETWEEN HDF5 AND ZARR:

HDF5: All the chunks are stored in the same file --> .h5 file
Zarr: Uses one file per chunk --> .zarr/ folder

Benchmarks

These 4 objects contain the first 50k cols of the "1.3 Million Brain Cell Dataset" from 10x Genomics (see section below for how these datasets were created):

## HDF5Array objects:
ha100  <- HDF5Array(ha100_path, "counts")   #  100x100  chunks --> 152M on disk
ha1000 <- HDF5Array(ha1000_path, "counts")  # 1000x1000 chunks --> 139M on disk

## ZarrArray objects:
za100  <- ZarrArray(za100_path)             #  100x100  chunks --> 551M on disk
za1000 <- ZarrArray(za1000_path)            # 1000x1000 chunks --> 187M on disk

dim(za1000)
# [1] 27998 50000

Reading a "spread out" selection of array elements is more snappy when the chunks are small:

index <- list(sample(nrow(za100), 100), sample(ncol(za100), 200))  # spread out

system.time(m1 <- extract_array(ha100,  index))
#   user  system elapsed 
#  0.602   0.013   0.615

system.time(m2 <- extract_array(ha1000, index))
#   user  system elapsed 
#  5.105   0.021   5.126

system.time(m3 <- extract_array(za100,  index))
#   user  system elapsed 
#  1.038   0.066   1.104

system.time(m4 <- extract_array(za1000, index))
#   user  system elapsed 
#  3.909   0.935   4.844

stopifnot(identical(m1, m2), identical(m1, m3), identical(m1, m4))

However reading a big block (i.e. "compact selection") is faster with when the chunks are big:

index <- list(4001:9000, 15001:35000)  # compact

system.time(m1 <- extract_array(ha100,  index))
#   user  system elapsed 
#  0.482   0.077   0.559

system.time(m2 <- extract_array(ha1000, index))
#   user  system elapsed 
#  0.430   0.001   0.431

system.time(m3 <- extract_array(za100,  index))
#   user  system elapsed
#  1.713   0.126   1.839

system.time(m4 <- extract_array(za1000, index))
#   user  system elapsed
#  0.397   0.165   0.563

stopifnot(identical(m1, m2), identical(m1, m3), identical(m1, m4))

How the 4 datasets used in the benchmarks above were created

library(HDF5Array)
library(ExperimentHub)
hub <- ExperimentHub()
hub["EH1039"]$description  # 1.3 Million Brain Cell Dataset (sparse)
brain_s_path <- hub[["EH1039"]]
brain_s <- TENxMatrix(brain_s_path, group="mm10")

Using writeHDF5Array() to create the two HDF5 datasets:

ha100_path <- "brain_first_50k_cols_100x100_chunks.h5"
system.time(writeHDF5Array(brain_s[ , 1:50000], ha100_path, "counts", chunkdim=c(100, 100), with.dimnames=FALSE, verbose=TRUE))
#    user  system elapsed
#  44.607   1.727  46.344

ha1000_path <- "brain_first_50k_cols_1000x1000_chunks.h5"
system.time(writeHDF5Array(brain_s[ , 1:50000], ha1000_path, "counts", chunkdim=c(1000, 1000), with.dimnames=FALSE, verbose=TRUE))
#    user  system elapsed 
#  54.420   1.435  55.903

Using writeZarrArray() to create the two Zarr datasets (note that writeZarrArray() is painfully slow when the chunk geometry is set to 100x100):

library(ZarrArray)

za100_path <- "./brain_first_50k_cols_100x100_chunks.zarr"
system.time(writeZarrArray(brain_s[ , 1:50000], za100_path, chunkdim=c(100, 100), verbose=TRUE))
#     user  system  elapsed 
#  773.777   9.710  784.013 

za1000_path <- "./brain_first_50k_cols_1000x1000_chunks.zarr"
system.time(writeZarrArray(brain_s[ , 1:50000], za1000_path, chunkdim=c(1000, 1000), verbose=TRUE))
#     user  system  elapsed
#   26.762   1.631   28.968

sessionInfo():

R version 4.6.0 alpha (2026-04-05 r89793)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS

Matrix products: default
BLAS:   /home/hpages/R/R-4.6.r89793/lib/libRblas.so 
LAPACK: /home/hpages/R/R-4.6.r89793/lib/libRlapack.so;  LAPACK version 3.12.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB              LC_COLLATE=en_US.UTF-8    
 [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       

time zone: America/Los_Angeles
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ZarrArray_0.99.3      ExperimentHub_3.1.0   AnnotationHub_4.1.0  
 [4] BiocFileCache_3.1.0   dbplyr_2.5.2          HDF5Array_1.39.0     
 [7] h5mread_1.3.3         rhdf5_2.55.16         DelayedArray_0.37.1  
[10] SparseArray_1.11.13   S4Arrays_1.11.1       IRanges_2.45.0       
[13] abind_1.4-8           S4Vectors_0.49.1      MatrixGenerics_1.23.0
[16] matrixStats_1.5.0     BiocGenerics_0.57.0   generics_0.1.4       
[19] Matrix_1.7-5         

loaded via a namespace (and not attached):
 [1] KEGGREST_1.51.1      httr2_1.2.2          Biobase_2.71.0      
 [4] lattice_0.22-9       rhdf5filters_1.23.3  vctrs_0.7.2         
 [7] tools_4.6.0          curl_7.0.0           tibble_3.3.1        
[10] AnnotationDbi_1.73.0 RSQLite_2.4.6        paws.common_0.8.9   
[13] blob_1.3.0           pkgconfig_2.0.3      R.oo_1.27.1         
[16] lifecycle_1.0.5      compiler_4.6.0       Biostrings_2.79.5   
[19] Seqinfo_1.1.0        yaml_2.3.12          pillar_1.11.1       
[22] crayon_1.5.3         R.utils_2.13.0       cachem_1.1.0        
[25] digest_0.6.39        tidyselect_1.2.1     paws.storage_0.9.0  
[28] dplyr_1.2.1          BiocVersion_3.23.1   fastmap_1.2.0       
[31] grid_4.6.0           cli_3.6.5            Rarr_1.11.35        
[34] magrittr_2.0.5       filelock_1.0.3       rappdirs_0.3.4      
[37] bit64_4.6.0-1        XVector_0.51.0       httr_1.4.8          
[40] bit_4.6.0            png_0.1-9            R.methodsS3_1.8.2   
[43] memoise_2.0.1        rlang_1.2.0          Rcpp_1.1.1          
[46] glue_1.8.0           DBI_1.3.0            BiocManager_1.30.27 
[49] jsonlite_2.0.0       R6_2.6.1             Rhdf5lib_1.33.6

Zarr vs HDF5: Time of extract_array() as a function of index size

Generate the plot data:

library(HDF5Array)
library(ZarrArray)
ha100_path <- "brain_first_50k_cols_100x100_chunks.h5"
za100_path <- "./brain_first_50k_cols_100x100_chunks.zarr"
ha100  <- HDF5Array(ha100_path, "counts")
za100  <- ZarrArray(za100_path)

set.seed(111)
index_list <- lapply(1:12,
    function(i) {
        n <- 2L ** (i - 1L)
        list(sample(nrow(za100), 12 * n), sample(ncol(za100), 20 * n))
    }
)

index_size <- integer(length(index_list))
T1 <- T3 <- numeric(length(index_list))
for (i in seq_along(index_list)) {
    index <- index_list[[i]]
    index_size[[i]] <- prod(lengths(index))
    cat("index size = ", index_size[[i]], "\n", sep="")
    st1 <- system.time(m1 <- extract_array(ha100,  index))
    print(st1)
    st3 <- system.time(m3 <- extract_array(za100,  index))
    print(st3)
    stopifnot(identical(m1, m3))
    T1[[i]] <- st1[["elapsed"]]
    T3[[i]] <- st3[["elapsed"]]
}

Plot the data:

library(ggplot2)
index_coverage <- 100 * index_size / length(za100)
df1 <- data.frame(index_coverage=index_coverage, t=T1, format="HDF5")
df2 <- data.frame(index_coverage=index_coverage, t=T3, format="Zarr")
df <- rbind(df1, df2)
myplot <- ggplot(df, aes(x=index_coverage, y=t, color=format)) +
    geom_line() +
    scale_x_log10() + scale_y_log10() +
    labs(title="Time of extract_array() as a function of index size (100x100 chunks)",
         x="% of array covered by index", y="time in seconds")
myplot

ggsave("zzz_extract_array_timings_100x100_chunks.pdf", plot=myplot, width=7, height=5, units="in", dpi=150)