faissR: nearest-neighbor search for Bioconductor workflows

Introduction

faissR provides exact and approximate nearest-neighbor search for dense row-wise matrices. A functional installation uses the FAISS C++ library on CPU; optional CUDA, FAISS GPU, and NVIDIA cuVS providers support GPU execution. The same metric and result conventions are used through the public R API.

This vignette follows one practical workflow:

  1. load a dense single-cell representation;
  2. inspect the installed capabilities;
  3. run exact and HNSW searches;
  4. compare observed HNSW recall with an exact reference;
  5. reuse an index for prediction; and
  6. preflight input and metric edge cases.

All computational chunks are evaluated. CUDA code is evaluated conditionally and prints a skip message when the package was built without CUDA support. Detailed build instructions are available in vignette("installation", package = "faissR").

Example data

The installed example is derived from the public Zeisel mouse-brain single-cell RNA-sequencing data distributed by scRNAseq. It contains 50 PCA coordinates and cell annotations for 3,005 cells, but no expression matrix. The source URLs, complete preprocessing procedure, seed, and package versions are recorded by inst/scripts/prepare_zeisel_brain_pca.R and documented in ?zeisel_brain_pca.

library("faissR")
library("SingleCellExperiment")

example_file <- system.file(
    "extdata", "zeisel_brain_pca.rds", package = "faissR"
)
sce <- readRDS(example_file)
x <- reducedDim(sce, "PCA")
storage.mode(x) <- "double"
groups <- droplevels(factor(colData(sce)$level1class))

dim(x)
#> [1] 3005   50
head(table(groups), 6L)
#> groups
#> astrocytes_ependymal    endothelial-mural         interneurons 
#>                  224                  235                  290 
#>            microglia     oligodendrocytes        pyramidal CA1 
#>                   98                  820                  939

Nearest-neighbor search should normally receive a dense low-dimensional representation such as PCA scores, not a full sparse gene-by-cell count matrix. The following table estimates the principal host allocations for 50 PCs and k = 30 without allocating the larger matrices.

cells <- c(3005, 100000, 1000000)
data.frame(
    cells = cells,
    input_double_MiB = round(8 * cells * 50 / 1024^2, 1),
    native_float32_MiB = round(4 * cells * 50 / 1024^2, 1),
    result_MiB = round((4 + 8) * cells * 30 / 1024^2, 1)
)
#>     cells input_double_MiB native_float32_MiB result_MiB
#> 1    3005              1.1                0.6        1.0
#> 2  100000             38.1               19.1       34.3
#> 3 1000000            381.5              190.7      343.3

The result estimate includes int32 identifiers and R double distances. Exact full self-search also has quadratic computational cost, so large studies will usually require an approximate method, an exact validation subsample, and an explicit memory budget.

Inspect execution capabilities

Hardware visibility and package capability are different questions. The terminology here follows the framework-independent distinction used by the gpuinfo project: a visible GPU or driver does not prove that a consuming R package was compiled with a usable GPU provider.

In faissR:

  • faissR_backend() reads or sets the requested session device;
  • backend_info() reports providers compiled into this installation and whether their runtimes are currently usable;
  • nn_capabilities(runtime = TRUE) checks a method, device, and metric combination; and
  • cuda_available() is package-specific: it requires both a CUDA-enabled faissR build and a runtime-visible device.

The separate gpuinfo project can provide platform-level hardware diagnosis. faissR does not currently import it; package-specific provider checks remain authoritative for deciding whether a faissR call can execute.

faissR_backend()
#> [1] "cpu"

backend_info()[, c(
    "backend", "available", "device"
)]
#>          backend available                device
#> 1            cpu      TRUE CPU (4 logical cores)
#> 2          faiss      TRUE                   CPU
#> 3 faiss_gpu_cuvs     FALSE                  <NA>
#> 4           cuvs     FALSE              CUDA GPU
#> 5           cuda     FALSE                  <NA>

cpu_euclidean <- subset(
    nn_capabilities(runtime = TRUE),
    backend == "cpu" & metric == "euclidean"
)
cpu_euclidean[, c(
    "method", "supported", "runtime_available", "resolved_backend"
)]
#>             method supported runtime_available     resolved_backend
#> 1             auto      TRUE              TRUE             cpu_auto
#> 7            exact      TRUE              TRUE        faiss_flat_l2
#> 13            flat      TRUE              TRUE        faiss_flat_l2
#> 19      bruteforce      TRUE              TRUE        faiss_flat_l2
#> 25            grid      TRUE              TRUE             cpu_grid
#> 31            hnsw      TRUE              TRUE           faiss_hnsw
#> 37             ivf      TRUE              TRUE            faiss_ivf
#> 43           ivfpq      TRUE              TRUE          faiss_ivfpq
#> 49    vamana_style      TRUE              TRUE           cpu_vamana
#> 55       nsg_style      TRUE              TRUE              cpu_nsg
#> 61 nndescent_style      TRUE              TRUE        cpu_nndescent
#> 67  ivfpq_fastscan      TRUE              TRUE faiss_ivfpq_fastscan
#> 73           cagra     FALSE             FALSE                 <NA>

CPU is the default. An explicit function argument has highest precedence, followed by options(faissR.backend), FAISSR_BACKEND, and finally CPU.

old_backend <- faissR_backend("cpu")
stopifnot(identical(faissR_backend(), "cpu"))

An explicit backend = "cuda" never silently falls back to CPU. faissR does not provide a Metal backend.

Approximate search and observed recall

HNSW trades exactness for speed and index storage. tuning = "auto" chooses parameters for the resolved method. target_recall accepts the calibrated tiers 0.90, 0.95, and 0.99; it is not interpolated. For a new application, measure observed recall on representative queries against an exact reference.

hnsw_time <- system.time(
    hnsw <- nn(
        x,
        k = 15,
        backend = "cpu",
        method = "hnsw",
        metric = "euclidean",
        tuning = "auto",
        target_recall = 0.99,
        exclude_self = TRUE,
        n_threads = 2
    )
)[["elapsed"]]

query_recall <- vapply(
    seq_len(nrow(exact$indices)),
    function(i) mean(hnsw$indices[i, ] %in% exact$indices[i, ]),
    numeric(1L)
)

c(
    elapsed_seconds = unname(hnsw_time),
    mean_recall = mean(query_recall),
    fifth_percentile_recall = unname(quantile(query_recall, 0.05)),
    minimum_recall = min(query_recall),
    backend_used = attr(hnsw, "backend_used")
)
#>         elapsed_seconds             mean_recall fifth_percentile_recall 
#>    "0.0629999999999997"     "0.999534109816972"                     "1" 
#>          minimum_recall            backend_used 
#>     "0.866666666666667"            "faiss_hnsw"

Automatic selection

The execution device is selected with backend = "cpu" or backend = "cuda". If backend is omitted, faissR follows the package option, then FAISSR_BACKEND, and finally uses CPU. Two other arguments control different decisions:

Argument Decision
method = "auto" choose a method family
tuning = "auto" choose parameters within that method

method = "auto" is a static, calibration-informed experimental policy. It does not benchmark the current data during an ordinary call. CPU automatic selection has less independent validation than the CUDA policy; use an explicit method when reproducibility, memory, or latency constraints require a fixed route.

auto_result <- nn(
    x,
    k = 8,
    backend = "cpu",
    method = "auto",
    metric = "euclidean",
    exclude_self = TRUE,
    n_threads = 2
)

c(
    requested_method = attr(auto_result, "requested_method"),
    selected_method = attr(auto_result, "auto_selection")$predicted_method,
    backend_used = attr(auto_result, "backend_used"),
    policy_status = attr(auto_result, "auto_policy_status")
)
#>                                   requested_method 
#>                                             "auto" 
#>                                    selected_method 
#>                                            "exact" 
#>                                       backend_used 
#>                                              "cpu" 
#>                                      policy_status 
#> "calibration_informed_not_independently_validated"

Metrics and edge cases

The public metric set is "euclidean", "cosine", and "correlation". Euclidean results are ordinary L2 distances, not squared values. Cosine normalizes rows; correlation centers and then normalizes rows. All routes reject NA, NaN, and infinite input.

Zero vectors have undefined cosine distance, and constant rows have undefined correlation distance. CPU routes apply documented finite conventions; CUDA routes reject those rows because no CPU repair is performed. Use nn_metric_preflight() before an expensive call to identify affected one-based row indices and the backend action.

edge_x <- rbind(c(0, 0), c(1, 0), c(1, 1))
nn_metric_preflight(
    edge_x,
    metric = "cosine",
    backend = "cuda"
)
#> $metric
#> [1] "cosine"
#> 
#> $requested_backend
#> [1] "cuda"
#> 
#> $self_query
#> [1] TRUE
#> 
#> $degenerate_kind
#> [1] "zero_vector"
#> 
#> $data_degenerate_rows
#> [1] 1
#> 
#> $points_degenerate_rows
#> [1] 1
#> 
#> $data_rows
#> [1] 1
#> 
#> $points_rows
#> [1] 1
#> 
#> $data_non_finite_rows
#> integer(0)
#> 
#> $points_non_finite_rows
#> integer(0)
#> 
#> $has_degenerate_rows
#> [1] TRUE
#> 
#> $has_non_finite_rows
#> [1] FALSE
#> 
#> $would_succeed
#> [1] FALSE
#> 
#> $action
#> [1] "error_degenerate_cuda"
#> 
#> $cpu_convention
#> [1] "two_degenerate_rows_distance_0;mixed_rows_distance_1"
#> 
#> $row_index_base
#> [1] 1

Unsupported method, backend, and metric combinations fail explicitly. The capability table can be used to filter a workflow before execution.

tryCatch(
    nn(x, k = 5, backend = "cpu", method = "cagra"),
    error = function(e) conditionMessage(e)
)
#> [1] "`method = \"cagra\"` is only available with `backend = \"cuda\"`."

Dense input contract

faissR searches dense, two-dimensional numeric data only. Accepted inputs are numeric matrices, numeric data frames, and optional float::fl() matrices, with observations in rows. It does not implement sparse nearest- neighbor search. Sparse, delayed, file-backed, and other matrix-like objects are rejected before conversion and are never silently materialized with as.matrix().

sparse_example <- Matrix::Matrix(diag(4), sparse = TRUE)
tryCatch(
    nn(sparse_example, k = 1),
    error = function(e) conditionMessage(e)
)
#> [1] "`data` is a ddiMatrix object. faissR accepts dense numeric matrices, numeric data frames, and optional float::fl() matrices only. Sparse, delayed, and file-backed objects are not converted implicitly because `as.matrix()` can allocate a dense copy. Compute a size-checked dense representation such as PCA scores, or explicitly materialize a manageable block before calling faissR."

For single-cell data, keep the expression or count assay sparse during normalization and dimension reduction, then pass a dense, size-checked cell embedding such as PCA coordinates to faissR. If direct materialization is appropriate for another workflow, it must be performed explicitly by the caller before the faissR call.

Ordinary R double matrices are converted internally to contiguous float32 storage for FAISS and cuVS. This numeric conversion applies only after the input has passed the dense-input check. The original R object remains live during conversion, so the double source, float32 payload, index workspace, and result must all fit in host memory.

Reuse an index for prediction

knn() creates a fitted classification or regression model. Supported FAISS indexes are retained for later predict() calls in the same R session. Native indexes are session-local external pointers: after saveRDS() and readRDS(), the R model remains usable but the native index is rebuilt on the first prediction.

build_time <- system.time(
    fit <- knn(
        Xtrain = x,
        Ytrain = groups,
        k = 7,
        backend = "cpu",
        method = "exact",
        metric = "euclidean",
        n_threads = 2
    )
)[["elapsed"]]

query_time <- system.time(
    prediction <- predict(fit, x[1:12, , drop = FALSE])
)[["elapsed"]]

c(
    build_seconds = unname(build_time),
    query_seconds = unname(query_time),
    query_source = attr(prediction, "faissR_nn")$query_source,
    backend_used = attr(prediction, "faissR_nn")$resolved_backend
)
#>         build_seconds         query_seconds          query_source 
#> "0.00400000000000045" "0.00399999999999956"        "fitted_index" 
#>          backend_used 
#>       "faiss_flat_l2"
head(as.vector(prediction))
#> [1] "interneurons" "interneurons" "interneurons" "interneurons" "interneurons"
#> [6] "interneurons"

The imported Biobase package is supported for conventional ExpressionSet workflows. Samples are stored in columns, so transpose the assay before search.

data("sample.ExpressionSet", package = "Biobase")
biobase_x <- t(Biobase::exprs(sample.ExpressionSet)[seq_len(32L), ])
biobase_result <- nn(
    biobase_x,
    k = 3,
    backend = "cpu",
    method = "exact",
    exclude_self = TRUE
)
dim(biobase_result$indices)
#> [1] 26  3

Optional CUDA execution

Guard CUDA work with the package-specific cuda_available() predicate. The conditional code below is evaluated on every builder; only the GPU operation is skipped when CUDA is unavailable.

if (cuda_available()) {
    cuda_result <- nn(
        x,
        k = 15,
        backend = "cuda",
        method = "exact",
        metric = "euclidean",
        exclude_self = TRUE
    )
    c(
        backend_used = attr(cuda_result, "backend_used"),
        device = attr(cuda_result, "device")
    )
} else {
    cat(
        "CUDA example skipped: this faissR installation has ",
        "no usable CUDA route."
    )
}
#> CUDA example skipped: this faissR installation has  no usable CUDA route.

nn_gpu() serves a different workflow: it returns an owning faissR_gpu_knn object whose neighbor buffers stay on the device until gpu_knn_to_host() is called. This is intended for compiled downstream CUDA consumers. The object must remain alive while non-owning device pointers are in use and cannot be serialized across sessions or devices. See ?nn_gpu and ?faissR_c_api for the ownership and C-callable interface contracts.

candidate_knn() and fast_kmeans() are secondary interfaces that reuse the same metric and backend infrastructure. Their complete contracts and examples are available in ?candidate_knn and ?fast_kmeans.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        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: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#>  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
#>  [5] Seqinfo_1.3.2               IRanges_2.47.5             
#>  [7] S4Vectors_0.51.10           BiocGenerics_0.59.12       
#>  [9] generics_0.1.4              MatrixGenerics_1.25.0      
#> [11] matrixStats_1.5.0           faissR_0.99.47             
#> [13] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-6        jsonlite_2.0.0      compiler_4.6.1     
#>  [4] BiocManager_1.30.27 Rcpp_1.1.2          parallel_4.6.1     
#>  [7] jquerylib_0.1.4     yaml_2.3.12         fastmap_1.2.0      
#> [10] lattice_0.23-1      XVector_0.53.0      R6_2.6.1           
#> [13] S4Arrays_1.13.1     knitr_1.52          DelayedArray_0.39.6
#> [16] maketools_1.3.2     bslib_0.12.0        rlang_1.3.0        
#> [19] cachem_1.1.0        xfun_0.61           sass_0.4.10        
#> [22] sys_3.4.3           otel_0.2.0          SparseArray_1.13.2 
#> [25] cli_3.6.6           digest_0.6.39       grid_4.6.1         
#> [28] lifecycle_1.0.5     evaluate_1.0.5      buildtools_1.0.0   
#> [31] abind_1.4-8         rmarkdown_2.32      tools_4.6.1        
#> [34] htmltools_0.5.9