implement dot for higher-dimensional arrays with Ix2

[EbrahimEldesoky]

closes #1587

right now calling .dot() on anything above Ix2 doesn't work, so if you have a 3-D or 5-D array and want to multiply it by a matrix you have to manually reshape, dot, then reshape back.

this adds Dot<ArrayRef<A, Ix2>> for Ix3, Ix4, Ix5, Ix6, and IxDyn. the implementation flattens all axes except the last into a single "rows" dimension, delegates to the existing 2-D dot (which already handles BLAS), then reshapes the result back to the correct output shape.

for example:

use ndarray::prelude::*;
use ndarray::linalg::Dot;

let x = Array5::<f64>::zeros((3, 2, 5, 9, 12));
let y = Array2::<f64>::zeros((12, 13));
let z = x.dot(&y); // shape: (3, 2, 5, 9, 13)

[EbrahimEldesoky]

hi @akern40 please when you have time check this pr

[nilgoyette]

One thing you could do before someone checks your PR is remove the formatting changes. It's hard[er] to know what you actually modified. We should probably format the whole project with a specific rule but it hasn't been done yet.

[EbrahimEldesoky]

One thing you could do before someone checks your PR is remove the formatting changes. It's hard[er] to know what you actually modified. We should probably format the whole project with a specific rule but it hasn't been done yet.

One thing you could do before someone checks your PR is remove the formatting changes. It's hard[er] to know what you actually modified. We should probably format the whole project with a specific rule but it hasn't been done yet.

Thank you @nilgoyette, Please take a look.

[nilgoyette]

I don't know how common this operation is used (I never had to use it). If it is common, then yes we should add this feature to ndarray. Is it the same behavior as numpy dot?

[EbrahimEldesoky]

I don't know how common this operation is used (I never had to use it). If it is common, then yes we should add this feature to ndarray. Is it the same behavior as numpy dot?

Yep, this lines up exactly with how numpy's dot works — for an N-D array against a 2-D matrix, it takes the sum-product over the last axis of the first array and the first axis of the second.
It's a pattern that shows up constantly in ML/data science too, like passing batched inputs shaped (Batch, Seq, Features) through a (Features, Outputs) weight matrix. Right now you'd have to flatten, multiply, then reshape everything back manually, which gets tedious fast and is an easy place to introduce bugs

[nilgoyette]

I'm not a fan of the reference_nd_dot function, but I understand that not using it would force you to hardcode long and unreadable arrays.

Good job. LGTM.

[EbrahimEldesoky]

I'm not a fan of the reference_nd_dot function, but I understand that not using it would force you to hardcode long and unreadable arrays.

Good job. LGTM.

thank you

[bluss]

The relevant numpy documentation for dot is here:

https://numpy.org/doc/stable/reference/generated/numpy.dot.html

---

I'm basically not the maintainer here anymore - I think this is a good idea if we think enough about the interface and how to get the best out of performance. We should avoid whole-array-copying traps if possible.

The existing use of to_shape is pretty good for this, but will unfortunately cause copying in many cases.

Maybe this is a good enough compromise for now (with some improvements)?

---

In theory we could do it without any extra copying using matrixmultiply or blas like this:

Numpy had their example like this:

dot(a, b)[i,j,m] = sum(a[i,j,:] * b[:,m])

If axes i, j are contiguous together we can fold them together and compute as:

dot(a, b)[ij,m] = sum(a[ij,:] * b[:,m])

where ij is the combined axis. That's what's happening in this PR, more or less.

If axes i, j are not contiguous together (we can't reshape away that axis for free), then you could repair that by breaking up that axis in to multiple reductions, like this:

dot(a, b)[0,j,m] = sum(a[0,j,:] * b[:,m])  // sgemm i=0
dot(a, b)[1,j,m] = sum(a[1,j,:] * b[:,m])  // sgemm i=1
...         // and so on for the lenght of the axis

and then the operation still doesn't need extra copying of the ndarray. (of course both blas and matrixmultiply will copy the elements, efficiently, to their packing buffers anyway, but that's a necessary part of the matrix multiplication implementation, and usually quite optimized.)

---

Also for interest, general tensordot is already implemented for ndarray by the crate https://docs.rs/ndarray-einsum/ but I don't know its implementation. I think it actually might be a bit similar to the solution in this PR.

[EbrahimEldesoky]

@bluss, thanks for the guidance Please take a look at the updates Is there anything else?