Outer product

[zgbkdlm]

I found that ArrayBase has no equivalent method to np.outer of numpy. This is quite surprising because outer product is a standard feature and is trivial to implement.

For example, by refering to the implementation of np.outer, we can implement ArrayBase's outer as:

pub fn outer(x: &Array<f64, Ix1>, y: &Array<f64, Ix1>) -> Array<f64, Ix2> {
    let (size_x, size_y) = (x.shape()[0], y.shape()[0]);
    let x_reshaped = x.view().into_shape((size_x, 1)).unwrap();
    let y_reshaped = y.view().into_shape((1, size_y)).unwrap();
    x_reshaped.dot(&y_reshaped)
}

then

let a = array![1., 2., 3.];
let b = array![1., 2., 3.];

let outer = a.outer(&b);

There is a similar issue/pull request #690, but it implements Kronecker product not outer.

[jturner314]

A useful generalization is numpy.ufunc.outer. That would be more work to implement efficiently, though, because we couldn't reuse the matrix multiplication implementation.

[vladfaust]

For 2 dimensional arrays:

pub fn outer_product2(a: &ArrayView2<f32>, b: &ArrayView2<f32>) -> Array4<f32> {
    let (m, n) = a.dim();
    let (p, q) = b.dim();

    let a_reshaped = a.to_shape((m, n, 1, 1)).unwrap();
    let b_reshaped = b.to_shape((1, 1, p, q)).unwrap();

    let prod = &a_reshaped * &b_reshaped;

    prod.to_owned()
}

#[test]
fn test_outer_product2() {
    // -6	10
    // -8	5
    let a = array![[-6., 10.], [-8., 5.]];

    // 7	-6
    // 0	-2
    let b = array![[7., -6.], [0., -2.]];

    // |-42 36|	|70 -60|
    // |0   12|	|0  -20|
    // -------- --------
    // |-56 48| |35 -30|
    // |0   16| |0  -10|
    let expected = array![
        [[[-42., 36.], [0., 12.]], [[70., -60.], [0., -20.]]],
        [[[-56., 48.], [0., 16.]], [[35., -30.], [0., -10.]]]
    ];

    let result = outer_product2(&a.view(), &b.view());
    assert_eq!(result, expected);

    // |-42 70| |36 -60|
    // |-56 35| |48 -30|
    // -------- --------
    // |0    0| |12 -20|
    // |0    0| |16 -10|
    let expected = array![
        [[[-42., 70.], [-56., 35.]], [[36., -60.], [48., -30.]]],
        [[[0., 0.], [0., 0.]], [[12., -20.], [16., -10.]]]
    ];

    let result = outer_product2(&b.view(), &a.view());
    assert_eq!(result, expected);
}