lazymatrix 0.4.0

Lazy normalized design matrices: apply column centering/scaling without materializing (X - 1cᵀ)S⁻¹, preserving sparsity.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

//! Stable weighted sums for normalized columns.

use crate::{Scalar, VectorView, VectorViewMut};

pub(crate) fn validate<F, W, O>(
    nrows: usize,
    ncols: usize,
    weights: &W,
    centers: Option<&[F]>,
    scales: Option<&[F]>,
    out: &O,
) where
    F: Scalar,
    W: VectorView<F> + ?Sized,
    O: VectorViewMut<F> + ?Sized,
{
    assert_eq!(
        weights.len(),
        nrows,
        "weighted_column_sums_into: weights length mismatch"
    );
    assert_eq!(
        out.len(),
        ncols,
        "weighted_column_sums_into: output length mismatch"
    );
    crate::gram::validate_normalization(ncols, centers, scales);
}

#[cfg(any(feature = "faer_all", feature = "nalgebra_all", feature = "sprs_all"))]
pub(crate) fn sparse_sums<F, M, W, O>(
    matrix: &M,
    weights: &W,
    centers: Option<&[F]>,
    scales: Option<&[F]>,
    out: &mut O,
) where
    F: Scalar,
    M: crate::SparseColumns<F> + ?Sized,
    W: VectorView<F> + ?Sized,
    O: VectorViewMut<F> + ?Sized,
{
    use crate::gram::{WeightSums, normalize};

    let (n, p) = (matrix.nrows(), matrix.ncols());
    validate(n, p, weights, centers, scales, out);
    if n == 0 || p == 0 {
        for j in 0..p {
            out.set(j, F::zero());
        }
        return;
    }
    let finite_weights = (0..n).all(|i| weights.get(i).is_finite());
    let centered = (0..p).any(|j| normalize(F::zero(), j, centers, scales) != F::zero());
    let ranges = (centered && finite_weights).then(|| WeightSums::new(weights));
    let finite_ranges = ranges.as_ref().is_none_or(WeightSums::is_finite);
    let mut scratch = None;
    for j in 0..p {
        let (rows, values) = matrix.sparse_column(j);
        let background = normalize(F::zero(), j, centers, scales);
        let mut fast = finite_weights
            && finite_ranges
            && background.is_finite()
            && rows.windows(2).all(|pair| pair[0] < pair[1]);
        let mut sum = F::zero();
        if fast {
            let mut next = 0;
            for (&row, &value) in rows.iter().zip(values) {
                if next < row && background != F::zero() {
                    sum = sum + background * ranges.as_ref().unwrap().sum(next, row);
                }
                sum = sum + normalize(value, j, centers, scales) * weights.get(row);
                next = row + 1;
            }
            if next < n && background != F::zero() {
                sum = sum + background * ranges.as_ref().unwrap().sum(next, n);
            }
            fast = sum.is_finite();
        }
        if !fast {
            // Combine duplicate raw entries before normalization, and evaluate
            // implicit zeros individually to preserve nonfinite arithmetic.
            let column = scratch.get_or_insert_with(|| vec![F::zero(); n]);
            column.fill(F::zero());
            for (&row, &value) in rows.iter().zip(values) {
                column[row] = column[row] + value;
            }
            sum = (0..n)
                .map(|i| normalize(column[i], j, centers, scales) * weights.get(i))
                .sum();
        }
        out.set(j, sum);
    }
}