Sparse mat / dense vector product needs cleanup

[vbarrielle]

Right now this product is either available as a special case of the sparse/dense matrix multiplication, or as a free function that only takes slices as input.

We need to rewrite the free function to take ndarray vectors as input, and write the operator for sparse matrix/dense vector in term of this function.

[toastronics]

Recap, to make sure I understand what's being proposed here -

The change is related to definitions in src/sparse/prod.rs and src/sparse/csmat.rs.

"Special case of the sparse/dense matrix multiplication":

• impl Mul<&ArrayBase<Ix2>> for &CsMatBase;
• the free function it calls, {csr, csc}_mulacc_dense_rowmaj(lhs, rhs, out), where lhs is an N x M matrix;
• in the special case where rhs and out are M x 1 matrices.

"Free function that only takes slices as input": mul_acc_mat_vec_{csr, csc}

Proposed change:

• modify the interface of mul_acc_mat_vec_{csr,csc}(mat, in_vec, res_vec) such that in_vec and res_vec are ArrayBase<Ix1> instead of slices;
• add an impl Mul<&ArrayBase<Ix1>> for &CsMatBase using the modified mul_acc_mat_vec_{csr,csc}.

Is this understanding correct?

Another possible option, avoiding an API change, would be to keep mul_acc_mat_vec_{csr,csc} as-is and take advantage of the ndarray::ArrayView method into_slice(&self) -> Option<&'a [N]>. However, this method returns None if the ArrayView is not "contiguous and in standard order"; an impl Mul<&ArrayBase<Ix1>> for &CsMatBase which made use of into_slice() would not be as general as the type signature implies. Due to this restriction, the proposed change in the signature of mul_acc_mat_vec_{csr,csc} seems like the better option if this minor break in the API is acceptable; it is simple enough to construct an ArrayBase<Ix1> from a slice if required for any existing external callers of mul_acc_mat_vec_{csr,csc}.

A more general thought, related to this change: what about the ndarray Dot trait? I appreciate the extra convenience of being able to perform matrix multiplication as (A * B) * C, but it seems like the distinction between sprs "Mul = matrix multiplication" and ndarray "Mul = elementwise multiplication, Dot = matrix multiplication" could get confusing when mixing sprs and ndarray matrices. Obviously adding Dot and changing Mul in sprs to behave like ndarray's Mul would be a severe and often silently manifesting breaking change, but changing the sprs Mul to Dot without adding an element-wise Mul would not be a silent break. If this difference with ndarray behavior remains in sprs, it could use some explicit calling out.

[vbarrielle]

Proposed change:

* modify the interface of `mul_acc_mat_vec_{csr,csc}(mat, in_vec, res_vec)` such that `in_vec` and `res_vec` are `ArrayBase<Ix1>` instead of slices;

* add an `impl Mul<&ArrayBase<Ix1>> for &CsMatBase` using the modified `mul_acc_mat_vec_{csr,csc}`.

Is this understanding correct?

Yes that was roughly the idea, I thought it could also be possible to avoid having a true breaking change by having mul_acc_mat_vec_{csr,csc} be able to take either a slice or an ndarray::ArrayView of dimension 1. This could be done using a trait bound on the type of the arguments. Since this method only requires to be able to randomly access its input, I was originally intending a bound on std::ops::Index{Mut}, but that seems hard because I don't think ndarray supports indexing a one dimensional array using just integers.

Another possible option, avoiding an API change, would be to keep mul_acc_mat_vec_{csr,csc} as-is and take advantage of the ndarray::ArrayView method into_slice(&self) -> Option<&'a [N]>. However, this method returns None if the ArrayView is not "contiguous and in standard order"; an impl Mul<&ArrayBase<Ix1>> for &CsMatBase which made use of into_slice() would not be as general as the type signature implies. Due to this restriction, the proposed change in the signature of mul_acc_mat_vec_{csr,csc} seems like the better option if this minor break in the API is acceptable; it is simple enough to construct an ArrayBase<Ix1> from a slice if required for any existing external callers of mul_acc_mat_vec_{csr,csc}.

I think I'd rather have a minor break than have a risk of panic due to non-contiguity while calling a multiplication operator.

A more general thought, related to this change: what about the ndarray Dot trait? I appreciate the extra convenience of being able to perform matrix multiplication as (A * B) * C, but it seems like the distinction between sprs "Mul = matrix multiplication" and ndarray "Mul = elementwise multiplication, Dot = matrix multiplication" could get confusing when mixing sprs and ndarray matrices. Obviously adding Dot and changing Mul in sprs to behave like ndarray's Mul would be a severe and often silently manifesting breaking change, but changing the sprs Mul to Dot without adding an element-wise Mul would not be a silent break. If this difference with ndarray behavior remains in sprs, it could use some explicit calling out.

To me the difference is somehow natural, ndarray features arrays, not matrices, so it's unsurprising that element-wise multiplication is their default. But sprs is a matrix library, so it seemed natural to default to the matrix multiplication. You're right that it should be pointed out explicitly. An ergonomic way to leverage element-wise multiplication would be to be able to call e.g. my_mat.array_view() * other_mat.array_view(), similar to what's done in C++ in the eigen library.

[toastronics]

It turns out that ArrayBase<Ix1> does support indexing with just an integer, although it takes some digging to find it:

https://docs.rs/ndarray/0.12.0/ndarray/trait.NdIndex.html

This has impl NdIndex<Ix1> for Ix, where Ix is defined as type Ix = usize.

Unfortunately to due the requiment to use the len method in mul_acc_mat_vec_{csc,csr}, just adding Index and IndexMut bounds won't work, but a bound on something like the following trait would:

pub trait SliceLike<N>:
    Index<usize, Output = N> + IndexMut<usize, Output = N>
{
    fn len(&self) -> usize;
}

impl<N> SliceLike<N> for [N] {
    fn len(&self) -> usize {
        self.len()
    }
}

impl<N> SliceLike<N> for Vec<N> {
    fn len(&self) -> usize {
        self.len()
    }
}

impl<N, DS> SliceLike<N> for ArrayBase<DS, Ix1>
where
    DS: ndarray::Data<Elem = N> + ndarray::DataMut<Elem = N>,
{
    fn len(&self) -> usize {
        self.len()
    }
}

To me the difference is somehow natural, ndarray features arrays, not matrices, so it's unsurprising that element-wise multiplication is their default. But sprs is a matrix library, so it seemed natural to default to the matrix multiplication. You're right that it should be pointed out explicitly. An ergonomic way to leverage element-wise multiplication would be to be able to call e.g. my_mat.array_view() * other_mat.array_view(), similar to what's done in C++ in the eigen library.

Sounds good to me. On balance, having the more natural product behavior for 'sparse * sparse' and 'sparse * dense' probably outweighs concerns about mixing 'sparse * dense' and 'dense * dense'. Mostly I just wanted to make sure that this difference with ndarray was intentionally retained, since I couldn't find any discussion about it.

[vbarrielle]

Now that you mentionned it, I recall I've written a trait to express that a datatype is a vector of a given dimension: https://docs.rs/sprs/0.6.2/sprs/vec/trait.VecDim.html

It's probably a good idea to see if it is possible to use that and an IndexMut bound, as in https://github.com/vbarrielle/sprs/blob/master/src/sparse/linalg/trisolve.rs#L41

The trait would need to be extended to ndarray types though.

[toastronics]

I looked into implementing VecDim for ArrayBase<Ix1>. Due to the blanket implementation of VecDim for T: AsRef<[N]>, an impl VecDim<N> for ArrayBase<S, Ix1> mirroring impl DenseVector<N> for ArrayBase<S, Ix1> is not allowed (the compiler notes that ArrayBase may later impl AsRef<[N]>, causing the correct implementation of the trait to not be unique). This may be something that would be allowed in the future for with specialization, but for now an impl VecDim for ArrayBase<Ix1> would likely need to occur in the same way as DenseVector (with the AsRef impl replaced by explicit &[N] and Vec impls).

@vbarrielle, how would you like to proceed? Would changing the impl's of VecDim to be like DenseVector, with explicit impl's for

• &[N],
• Vec,
• and [N] (required to compile sparse/linalg/trisolve.rs L336)

and dropping the AsRef<[N]> impl, be an acceptable breaking change?

Alternatively, mul_acc_mat_vec_{csc,csr} could use a trait bound on DenseVector and IndexMut instead of VecDim and IndexMut.

[vbarrielle]

Hi @toastronics, sorry for answering so late.

Would changing the impl's of VecDim to be like DenseVector, with explicit impl's for

* `&[N]`,
* `Vec`,
* and `[N]` (required to compile `sparse/linalg/trisolve.rs` L336)

and dropping the AsRef<[N]> impl, be an acceptable breaking change?

That looks like an acceptable breaking change, since any code that relied on the AsRef blanket impl can simply call .as_ref() explicitly. Thanks for exploring this, I had forgotten about the lack of specialization.