[WIP] Outer product implementation #690

krtab · 2019-08-20T14:59:30Z

Hi !
This is my first contribution to ndarray ! :)
This deals with issue #652.

There are two proposed function here, more or less general. I think that a lot of choices can be made (should the f function be able to work on chunks or not), and efficiency wise, I am not too sure about what I have. Anyhow, I think it is a nice starting point, so let me know what you think about it, whether interface or implementation wise.

Also I didn't really know where to put it, so if you want to move the code somewhere else, please do.

bluss · 2019-08-20T18:02:36Z

src/linalg/impl_linalg.rs

+    Sb: Data<Elem = T>,
+    I: Dimension,
+    T: crate::ScalarOperand,
+    for<'a> &'a ArrayBase<Sb, I>: std::ops::Mul<T, Output = Array<T, I>>,


This seems inefficient since it's using a multiplication that produces an Array result, can we use a reusable scratch array for this or some other way to make this mostly in-place?

Maybe using Zip to make the assignment directly without using assign

Hi !

Yes indeed, this is not the most efficient. Yet, I haven't really found a function allowing me to directly perform a multiplication into a specified array/buffer. A possibility would be to use MulAssign but then I would have to add a pass to copy the first matrix and then mulassign each chunk, which I'm not sure would be more efficient. An other possibility is to just go and write ad-hoc code for the multiplication but it creates code duplication and if the std::ops implementation ever benefit from performance improvements (for example by using blas), these function won't.

Maybe using Zip to make the assignment directly without using assign

Not sure of what you mean here.

Let me know what you think, and thanks for your review.

Hey @bluss, I'm allowing myself to bump the discussion, do you have any opinion, let me know if something isn't clear.

krtab · 2020-03-03T19:01:28Z

This new API is both easier to use and also takes into account the comments of @bluss. Apart from function names, I'm satisfied with it. Let me know if you have any comments, and I'll start writing documentation!

krtab · 2020-04-19T11:21:09Z

@LukeMathWalker @jturner314 I'm not sure @bluss is around anymore, do you have any input ?

bluss · 2020-04-19T11:32:20Z

I won't be around much, but I'm trying to get 0.13.1 done when I'm here, and maybe get started on the next release.

Why do we have both generic dim and dynamic dim functions? If you replace the dimension type parameter with IxDyn, both functions have the same type signature. That should mean we only need one function of each kind.

Your methods are not callable by crate users, because they are not public. I suppose that tells us this is a draft PR, and that's ok.

bluss · 2020-04-19T11:33:49Z

In general, should this be in ndarray-linalg or here? I think as much as possible can go in ndarray-linalg, but have we defined the boundaries? Since this is more about the dimensionality of the operation and less about the linalg, maybe the ndarray crate is a good fit(?). Maybe @LukeMathWalker knows

bluss · 2020-04-19T11:36:32Z

src/linalg/impl_linalg.rs

+        .collect();
+
+    unsafe {
+        let mut res: ArrayD<T> = ArrayBase::uninitialized(res_dim);


If you have time, this can be ported to Array::maybe_uninit now. Look into using Zip::apply_assign_into - see example in tests? I'd still use T: Copy as the restriction, as long as we don't handle drop of partially filled arrays on panics.

src/linalg/impl_linalg.rs

krtab · 2020-04-20T16:48:36Z

Hi @bluss,

Thanks for your comments.

Why do we have both generic dim and dynamic dim functions? If you replace the dimension type parameter with IxDyn, both functions have the same type signature. That should mean we only need one function of each kind.

I made it so because in general, if you take the kronecker product of two matrices, one with n dimensions and the other with m dimensions, the dimensions of the result will be max(n,m). As the "max" operation is not implemented on types having the Dimension Trait, we have to loose the dimensionality information at the type-level. But the special case were both matrices have the same dimensionality is quite common, and allows us to preserve the dimensionality info at the type-level, hence my choice of two implementations.

Your methods are not callable by crate users, because they are not public. I suppose that tells us this is a draft PR, and that's ok.

Yes, in general I think it should be up to you the crate maintainers to decide the API.

I have implemented the modifications for both MaybeUninit, FnMut and apply_assign_into.

I'd still use T: Copy as the restriction, as long as we don't handle drop of partially filled arrays on panics.

I'm not sure I fully understand what you mean.

In general, I think you should double check my unsafe code, I'm pretty new to this so I'm not so confident it is ok.

bluss · 2020-04-20T17:13:22Z

In one case the dimensionality inputs are D and D (where D is a type parameter for the dimension) and the output dimensionality is D.

For the dynamic case, the dim of the inputs are IxDyn, IxDyn and the output has type IxDyn, so it looks like it can be handled by the first case (substitute D = IxDyn and they are the same). There might be something I'm missing.

Edit: (Reading more) seems like it will work, both operations are really the same, just adapting to the dimensionality. Making one interface with type parameter D will work, but it might be more tricky to code - and be like if <dynamic dimensional> { ... } else { .. }

bluss · 2020-04-20T17:20:16Z

src/linalg/impl_linalg.rs

+{
+    let mut res_dim = a.raw_dim();
+    let mut res_dim_view = res_dim.as_array_view_mut();
+    res_dim_view *= &b.raw_dim().as_array_view();


Obviously I kind of love this, that we can use array methods on dimensions 🙂.

The unsafe code looks good, we only need to prove that all elements are assigned to, and they will be if the exact chunks don't leave any uneven remainder. And that looks good to me, B's shape evenly divides the result's shape, right?

Ah no, we also have to check the product. We need to use a saturating multiplication, otherwise we can overflow, and then the above doesn't hold?

krtab · 2020-04-21T12:43:03Z

In one case the dimensionality inputs are D and D (where D is a type parameter for the dimension) and the output dimensionality is D.

For the dynamic case, the dim of the inputs are IxDyn, IxDyn and the output has type IxDyn, so it looks like it can be handled by the first case (substitute D = IxDyn and they are the same). There might be something I'm missing.

Edit: (Reading more) seems like it will work, both operations are really the same, just adapting to the dimensionality. Making one interface with type parameter D will work, but it might be more tricky to code - and be like if <dynamic dimensional> { ... } else { .. }

It is indeed quite tricky to code : I think the problem can be seen as: what should the type of res_dim be:

If we are working with dynamic dimensions we can collect into Vec<Ix> which is IntoDimension<Dim=IxDyn>, this allows us to leverage all the std iterator functions and even perform the checked multiplication that is needed as you pointed out.
If we are working with static Dimension, then res_dim should have the same type as the input ones, which is [Ix,n] (or a n-tuple of Ix which has the right IntoDimension). But we can't collect() into an [Ix,n] so that's why I had used the array multiplication trick. Unfortunately this trick won't allow the overflow check so we will have to do without it.

I ended up with the current solution which could replace the zeros with some uninitialized mem but that would require unsafe code I think.

I can also try to optimize a bit for the case where the number of dimensions are the same.

I also fixed a bug in the chunks/inputs shape compatibility for different ndims, and added two new tests for that.

bluss · 2020-04-21T14:23:26Z

I think both problems can be solved 😸 . Since I've written most of ndarray I have done a fair bit of weird dimension hacks. Not sure if it needs to be solved right now, I can't commit that much time right now, otherwise I'd do it.

Creating a new dimension value: D::zeros(n). Iterating it and filling it with its values by using .slice_mut() (shared ref iteration, .slice()).

bluss · 2020-04-21T14:36:12Z

src/linalg/impl_linalg.rs

+
+    // Reshape input arrays to compatible shapes
+    let a_reshape = a.view().into_shape(a_shape_completed).unwrap();
+    let b_reshape = b.view().into_shape(b_shape_completed.clone()).unwrap();


into_shape is unfortunately not general enough to always succeed here :( I'm sorry if I have pushed you over to this solution with into_shape. We'll need to think about what we can do, because we know the completed shape is compatible.

What's the issue with it? The layout constraint?

It only supports c/f-contiguous arrays, unfortunately.

The shape completed code can be much simpler I guess. Make views av and bv and use .insert_axis() until they both have the same number of axes.

hm no, insert_axis changes the type, that's no good :(

One can check if it's a dyn dimension - then convert to explicitly typed ArrayD and convert back to Array<_, D> later with .into_dimensionality()? Is that truly the best way? Not sure.

bluss · 2020-04-21T14:36:57Z

src/linalg/impl_linalg.rs

+            Zip::from(&b_reshape)
+                .apply_assign_into(res_chunk, |&b_elem| MaybeUninit::new(f(a_elem, b_elem)))
+        });
+    // This is safe because the exact chunks covered exactly the res


krtab · 2023-01-03T10:24:48Z

I'm closing this as too old. If anyone reworks on this feel free to let me know.

bluss reviewed Aug 20, 2019

View reviewed changes

LukeMathWalker mentioned this pull request Sep 8, 2019

Are we ndarray yet? #597

Open

31 tasks

krtab requested a review from bluss March 30, 2020 09:03

bluss reviewed Apr 19, 2020

View reviewed changes

src/linalg/impl_linalg.rs Outdated Show resolved Hide resolved

krtab added 6 commits April 20, 2020 17:20

First try for various outter product implementations

8c5738d

Added allow dead code while WIP to pass CI.

71129cb

Allow beta specific clippy:type_repetition_in_bounds

26b1087

One allocation less.

3d9933e

Tydied-up code a bit

566909e

New API, more clear and more efficient.

6dc1a67

krtab force-pushed the master branch from 0ed6fd2 to 6dc1a67 Compare April 20, 2020 15:21

krtab added 2 commits April 20, 2020 18:39

Changed to FnMut + MaybeUninit

551deeb

Use apply_assign_into

cc24d1e

bluss reviewed Apr 20, 2020

View reviewed changes

Change interface + fix bugs in different ndim case

d0098f4

bluss reviewed Apr 21, 2020

View reviewed changes

emmatyping mentioned this pull request Nov 6, 2021

Implement Kronecker product #1105

Merged

zgbkdlm mentioned this pull request Jan 20, 2022

Outer product #1148

Open

krtab closed this Jan 3, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Outer product implementation #690

[WIP] Outer product implementation #690

krtab commented Aug 20, 2019

bluss Aug 20, 2019

krtab Aug 20, 2019 •

edited

Loading

krtab Nov 13, 2019

krtab commented Mar 3, 2020

krtab commented Apr 19, 2020

bluss commented Apr 19, 2020

bluss commented Apr 19, 2020

bluss Apr 19, 2020

bluss Apr 21, 2020

krtab commented Apr 20, 2020 •

edited

Loading

bluss commented Apr 20, 2020 •

edited

Loading

bluss Apr 20, 2020

bluss Apr 20, 2020

krtab commented Apr 21, 2020 •

edited

Loading

bluss commented Apr 21, 2020

bluss Apr 21, 2020

krtab Apr 21, 2020

bluss Apr 21, 2020

bluss Apr 21, 2020

bluss Apr 21, 2020

bluss Apr 21, 2020

bluss Apr 21, 2020

krtab commented Jan 3, 2023

[WIP] Outer product implementation #690

[WIP] Outer product implementation #690

Conversation

krtab commented Aug 20, 2019

Choose a reason for hiding this comment

krtab Aug 20, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

krtab commented Mar 3, 2020

krtab commented Apr 19, 2020

bluss commented Apr 19, 2020

bluss commented Apr 19, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

krtab commented Apr 20, 2020 • edited Loading

bluss commented Apr 20, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

krtab commented Apr 21, 2020 • edited Loading

bluss commented Apr 21, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

krtab commented Jan 3, 2023

krtab Aug 20, 2019 •

edited

Loading

krtab commented Apr 20, 2020 •

edited

Loading

bluss commented Apr 20, 2020 •

edited

Loading

krtab commented Apr 21, 2020 •

edited

Loading