Sparse Statistics

Median, quantile, variance and standard-deviation reductions optimized for sparse vectors and matrices.

A quantile of a sparse vector with nnz non-zero values out of n total entries can be located in O(nnz) expected time instead of the dense O(n * log(n)) time spent on a full sort. The key observation is that the values, if conceptually sorted, lay out as [sorted negatives... | zeros... | sorted positives...] , where any explicit zero entries that may be stored in nzval are folded together with the structural zeros. Once we know how many negatives and how many positives are stored in nzval , mapping a quantile position to a value just requires running an in-place quickselect on (a copy of) nzval to extract the one or two values needed. The implementation is allocation-free when the caller provides a scratch buffer of the appropriate size.

The variance is computed in O(nnz) by running Welford's algorithm over nzval and merging the resulting stream stats (n, mean, M2) with the all-zeros stream (n_zero, 0, 0) via the standard pairwise variance combiner; this is exact and numerically stable, with no scratch buffer needed. The standard deviation is just the square root of the variance.

If the data is known to be non-negative (e.g., UMI counts), passing positive = true to the median/quantile functions skips the negativity scan entirely and yields a small additional speedup.

The functions exposed here mirror the calling convention of Statistics.var and Statistics.std : they accept either an AbstractVector (returning a single value) or an AbstractMatrix together with a dims keyword argument (returning a vector with one of the dimensions reduced to a single entry, or a scalar when dims is omitted).

The same [negatives... | zeros... | positives...] layout also speeds up the one-vs-rest AUC (area under the ROC curve) of a feature across labeled samples: only the stored non-zero values need to be ranked, the many zero samples sharing a single tied rank. See auc_per_group! and auc_per_category .

auc_per_category(
    feature_per_sample::AbstractVector{<:Real},
    category_per_sample::AbstractVector;
    positive::Bool = false,
)::Tuple{Vector, Vector{Float64}}

Compute the one-vs-rest AUC of the feature_per_sample for each distinct value of category_per_sample (an integer or string label per sample). Returns the distinct categories (in order of first appearance) aligned with their AUCs. See auc_per_group! for the definition and the sparse-feature optimization; use it directly to avoid re-deriving the categories when scoring many features against the same categories.

using SparseArrays

feature = sparse([0, 0, 0, 1, 2, 3])
category = [1, 2, 3, 1, 2, 3]
categories, aucs = auc_per_category(feature, category)
println(categories)
println(aucs)

# output

[1, 2, 3]
[0.625, 0.5, 0.625]

auc_per_group!(;
    auc_per_group::AbstractVector{Float64},
    feature_per_sample::AbstractVector{<:Real},
    group_index_per_sample::AbstractVector{<:Integer},
    n_samples_per_group::AbstractVector{<:Integer},
    positive::Bool = false,
)::AbstractVector{Float64}

Fill auc_per_group with the one-vs-rest AUC of the feature_per_sample , distinguishing each group (the samples whose group_index_per_sample holds that group's 1 -based index) from the rest of the included samples. Samples with a group index of 0 are excluded from the population. n_samples_per_group[group] is the number of samples in each group, supplied by the caller so that scoring many features against the same grouping computes it just once.

The AUC measures how predictive the feature is of membership in the group, regardless of direction: it is max(a, 1 - a) of the Mann-Whitney statistic a = (rank_sum - n * (n + 1) / 2) / (n * (N - n)) (tie-averaged ranks), where n is the group size and N the number of included samples. It is thus always >= 0.5 , so a group covering none or all of the included samples - where the AUC is undefined - is reported as 0 . Being rank based, any monotonic transform of the feature (e.g. a logarithm) yields the same result.

When feature_per_sample is sparse this costs O(nnz * log(nnz)) rather than the dense O(N * log(N)) : the many zero samples share the single tied rank of the [negatives... | zeros... | positives...] layout, so only the stored non-zero values are ranked. Pass positive = true when the feature is known to be non-negative (e.g. UMI counts) to skip the negative-value handling.

sparse_median(
    vector::AbstractVector{<:Real};
    positive::Bool = false,
    scratch::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
)::Float64

sparse_median(
    matrix::AbstractMatrix{<:Real};
    dims::Maybe{Integer} = nothing,
    positive::Bool = false,
    result::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    scratch::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    progress::Maybe{AbstractProgress} = nothing,
    progress_chunk::Maybe{Integer} = nothing,
)::Union{Float64, AbstractVector{<:AbstractFloat}}

Compute the median of the values of a vector , or of each column ( dims = 1 / Rows ) or each row ( dims = 2 / Columns ) of a matrix , or of all the elements of a matrix when dims is omitted. This is equivalent to sparse_quantile with p = 0.5 and shares its calling convention, optimizations, positive flag, optional result / scratch buffers, and progress / progress_chunk reporting.

using SparseArrays

vector = sparse([0, 0, 1, 2, 3])
println(sparse_median(vector))
println(sparse_median(vector; positive = true))

# output

1.0
1.0

using SparseArrays

matrix = sparse([0 1 0; 2 0 3; 0 4 0; 5 0 6])
sparse_median(matrix; dims = Rows)

# output

3-element Vector{Float64}:
 1.0
 0.5
 1.5

sparse_quantile(
    vector::AbstractVector{<:Real},
    p::Real;
    positive::Bool = false,
    scratch::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
)::Float64

sparse_quantile(
    matrix::AbstractMatrix{<:Real},
    p::Real;
    dims::Maybe{Integer} = nothing,
    positive::Bool = false,
    result::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    scratch::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    progress::Maybe{AbstractProgress} = nothing,
    progress_chunk::Maybe{Integer} = nothing,
)::Union{Float64, AbstractVector{<:AbstractFloat}}

Compute the p -quantile ( 0 <= p <= 1 ) of the values of a vector , or of each column ( dims = 1 / Rows ) or each row ( dims = 2 / Columns ) of a matrix , or of all the elements of a matrix when dims is omitted. With dims given, the matrix variant returns a vector of length n_columns ( dims = Rows ) or n_rows ( dims = Columns ), holding one quantile value per slice; without dims , it returns a scalar Float64 .

The quantile is computed via the same linear interpolation rule as Statistics.quantile (corresponding to R 's default "type 7"), so the result matches quantile(Vector(input), p) exactly. For a sparse vector (or each sparse column slice of a column-major sparse matrix , or each sparse row slice of a row-major sparse matrix ), the quantile is located in O(nnz) expected time using an in-place quickselect over a copy of the stored non-zero values.

If positive is true , all input values are assumed to be non-negative; this lets the implementation skip the scan that locates the boundary between negative and non-negative values. No validation is performed on this assumption.

If scratch is given, no allocation is performed. The required scratch length is nnz(input) for sparse data and length(input) for dense vectors (or n_rows * n_columns for dense matrices). If result is given for the matrix variant (length n_columns for dims = Rows , or n_rows for dims = Columns ), no result allocation is performed either.

When operating on a matrix , the iteration is parallelized via parallel_loop_wo_rng . Operating against the major_axis of a sparse matrix will trigger the GLOBAL_INEFFICIENT_ACTION_HANDLER and fall back to allocating per-iteration views.

If progress is given (matrix variant only), it is advanced once per processed slice; that is, by n_columns ticks when dims = Rows , and by n_rows ticks when dims = Columns . The progress_chunk is passed through to parallel_loop_wo_rng to throttle the rate of progress updates.

println(sparse_quantile([0, 0, 1, 2, 3], 0.5))
println(sparse_quantile([0, 0, 1, 2, 3], 0.5; positive = true))
println(sparse_quantile([-2, 0, 0, 1, 3], 0.25))

# output

1.0
1.0
0.0

using SparseArrays

matrix = sparse([0 1 0; 2 0 3; 0 4 0; 5 0 6])
sparse_quantile(matrix, 0.5; dims = Rows)

# output

3-element Vector{Float64}:
 1.0
 0.5
 1.5

sparse_quantile([1.0 4.0; 2.0 5.0; 3.0 6.0], 0.5; dims = Rows)

# output

2-element Vector{Float64}:
 2.0
 5.0

sparse_std(
    vector::AbstractVector{<:Real};
    corrected::Bool = true,
)::Float64

sparse_std(
    matrix::AbstractMatrix{<:Real};
    dims::Maybe{Integer} = nothing,
    corrected::Bool = true,
    result::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    progress::Maybe{AbstractProgress} = nothing,
    progress_chunk::Maybe{Integer} = nothing,
)::Union{Float64, AbstractVector{<:AbstractFloat}}

Compute the standard deviation of the values; equivalent to taking sqrt of sparse_var . Shares the calling convention, corrected flag, and optional result / progress / progress_chunk parameters of sparse_var .

using SparseArrays

matrix = sparse([0 1 0; 2 0 3; 0 4 0; 5 0 6])
sparse_std(matrix; dims = Rows)

# output

3-element Vector{Float64}:
 2.362907813126304
 1.8929694486000912
 2.8722813232690143

sparse_var(
    vector::AbstractVector{<:Real};
    corrected::Bool = true,
)::Float64

sparse_var(
    matrix::AbstractMatrix{<:Real};
    dims::Maybe{Integer} = nothing,
    corrected::Bool = true,
    result::Maybe{AbstractVector{<:AbstractFloat}} = nothing,
    progress::Maybe{AbstractProgress} = nothing,
    progress_chunk::Maybe{Integer} = nothing,
)::Union{Float64, AbstractVector{<:AbstractFloat}}

Compute the variance of the values of a vector , or of each column ( dims = 1 / Rows ) or each row ( dims = 2 / Columns ) of a matrix , or of all the elements of a matrix when dims is omitted. The shape conventions match Statistics.var : with dims the matrix variant returns a vector of length n_columns ( dims = Rows ) or n_rows ( dims = Columns ); without dims it returns a scalar Float64 .

If corrected is true (the default), the bias-corrected sample variance is returned (sum of squared deviations divided by n - 1 ); otherwise the population variance is returned (divided by n ).

The implementation runs Welford's algorithm over the stored non-zero values and merges the result with the all-zeros stream (n_zero, 0, 0) using the standard pairwise variance combiner. This is O(nnz) for sparse inputs, numerically stable, and allocation-free apart from the result vector. No scratch buffer is needed.

When operating on a matrix , the iteration is parallelized via parallel_loop_wo_rng . Operating against the major_axis of a sparse matrix will trigger the GLOBAL_INEFFICIENT_ACTION_HANDLER .

If progress is given (matrix variant only), it is advanced once per processed slice ( n_columns ticks for dims = Rows , n_rows for dims = Columns ). The progress_chunk is passed through to parallel_loop_wo_rng .

println(sparse_var([1.0, 2.0, 3.0, 4.0]))
println(sparse_var([1.0, 2.0, 3.0, 4.0]; corrected = false))

# output

1.6666666666666667
1.25

using SparseArrays

matrix = sparse([0 1 0; 2 0 3; 0 4 0; 5 0 6])
sparse_var(matrix; dims = Rows)

# output

3-element Vector{Float64}:
 5.583333333333333
 3.5833333333333335
 8.25