Hadamard product

[adinapoli-mndc]

Hey @vbarrielle !

First of all thanks for creating this fantastic library, we are successfully using it at work and it's doing what it advertises.

Lately we were in need of efficiently implementing Hadamard product on two sparse matrices. Theoretically speaking this could be possible by iterating over the rows of the first matrix, taking into account the non-zero indexes & data of both matrixes. Some (non-compiling) Rust code might look like this:

pub fn hadamard_mul<N>(mut lhs: CsMat<N>, rhs: &CsMat<N>) -> CsMat<N>
where
    N: Zero + PartialOrd + Signed + Copy,
{

    for (mut lhs_row_vec, rhs_row_vec) in lhs.outer_iterator_mut().zip(rhs.outer_iterator()) {

        // In order to correctly multiply the two matrixes, we need to consider
        // indices on both matrixes; iterating only over the indices of the LHS
        // would mean ignoring non-zero values for the RHS, and vice-versa.

        let all_ixs = ... // We can use BTreeSet.union to get all the non-zero indices from lhs_row_vec
                                // and rhs_row_vec

        for ix in all_ixs {
            ???
        }
    }

    lhs
}

The question would be what to put inside the for loop. The problem is that all the iterators the library offers (quite rightly) allow us only to mutate the content of a CsVec, not its internal structure, which is something we do want here: if we find a non-zero element belonging to the rhs_row_vec we want to insert it (in the correct position) inside the "internal storage" of lhs_row_vec.

Is there any hope to implement this using normal library combinators or one would have to either open a pull request extending this library or (even worse) resort to unsafe code? 😅

Thanks for your time!

[adinapoli-mndc]

Self-answering: the following Rust code snippet seems to do the trick (albeit by allocating a new matrix):

/// Actual implementation of CSR-CSR Hadamard product.
pub fn csr_hadamard_mul_csr_impl<N, I>(
    lhs: CsMatViewI<N, I>,
    rhs: CsMatViewI<N, I>,
    workspace: &mut [N],
) -> CsMatI<N, I>
where
    N: Num + Copy + Zero,
    I: SpIndex,
{
    let res_rows = lhs.rows();
    let res_cols = rhs.cols();
    if lhs.cols() != rhs.cols() {
        panic!("Dimension mismatch (cols)");
    }
    if lhs.rows() != rhs.rows() {
        panic!("Dimension mismatch (rows)");
    }
    if res_cols != workspace.len() {
        panic!("Bad storage dimension");
    }
    if lhs.storage() != rhs.storage() || !rhs.is_csr() {
        panic!("Storage mismatch");
    }

    let mut res = CsMatI::empty(lhs.storage(), res_cols);
    res.reserve_nnz_exact(lhs.nnz() + rhs.nnz());

    for (lrow_ix, lvec) in lhs.outer_iterator().enumerate() {
        // reset the accumulators
        for wval in workspace.iter_mut() {
            *wval = N::zero();
        }

        // accumulate the resulting row
        for (lcol_ix, &lval) in lvec.iter() {
            // we can't be out of bounds thanks to the checks of dimension
            // compatibility and the structure check of CsMat. Therefore it
            // should be safe to call into an unsafe version of outer_view

            let rvec = rhs.outer_view(lrow_ix).unwrap();
            let rval = match rvec.get(lcol_ix) {
                None => Zero::zero(),
                Some(v) => *v,
            };

            let wval = &mut workspace[lcol_ix];
            let prod = lval * rval;
            *wval = prod;
        }

        // compress the row into the resulting matrix
        res = res.append_outer(&workspace);
    }
    // TODO: shrink res storage? would need methods on CsMat
    assert_eq!(res_rows, res.rows());
    res
}

pub fn hadamard_mul<N>(lhs: &CsMat<N>, rhs: &CsMat<N>) -> CsMat<N>
where
    N: Zero + PartialOrd + Signed + Copy,
{
    let mut ws = workspace_csr(lhs, rhs);
    csr_hadamard_mul_csr_impl(lhs.view(), rhs.view(), &mut ws)
}

I haven't tested this thoroughly, but it's a start I guess.

[vbarrielle]

Hello @adinapoli-mndc, thanks for the kind words, I'm really glad this library suits you.

You should be able to use https://docs.rs/sprs/0.6.5/sprs/binop/fn.mul_mat_same_storage.html to compute the Hadamard product, though I must admit this is an under-tested part of the API. I'll taker some time to read your implementation and compare it to the current one.

[adinapoli-mndc]

You should be able to use https://docs.rs/sprs/0.6.5/sprs/binop/fn.mul_mat_same_storage.html to compute the Hadamard product, though I must admit this is an under-tested part of the API. I'll taker some time to read your implementation and compare it to the current one.

Ah, that's great to know, I love when somebody already did the hard work for me 😉

[vbarrielle]

There's a caveat though: I implemented this function a long time ago, and I'm not sure it's as performant as it could be.