Kronecker product

[emmatyping]

Myself (and another user on the rust discord) have run into the need for the kronecker product of two matrices.

I currently have code like:

fn kron(a: &Array2<f64>, b: &Array2<f64>) -> Array2<f64> {
    let dima = a.shape()[0];
    let dimb = b.shape()[0];
    let dimout = dima * dimb;
    let mut out = Array2::zeros((dimout, dimout));
    for (mut chunk, elem) in out.exact_chunks_mut((dimb, dimb)).into_iter().zip(a.iter()) {
        chunk.assign(&(*elem * b));
    }
    out
}

(note: in my code I only work on unitary matrices, so this is not general in a few ways)

It would be nice to have a general implementation of this in ndarray (numpy as np.kron).

[Jvanrhijn]

I myself have run into the need for an outer product of two vectors, which is of course a special case of the kronecker product. Looks like this used to exist in older versions of ndarray_linalg, but was removed at some point.

On a related note, numpy also has support for some other outer operations, i.e. if B is a binary operation between two elements of a field, and a and b are vectors over this field, then B(a, b)[i, j] == B(a[i], b[j]). This would be quite nice to have in ndarray as well. something like

fn outer_operation<T>(op: impl Fn(T, T) -> T, a: &Array1<T>, b: &Array1<T>) -> Array2<T>

[LukeMathWalker]

@termoshtt can you provide the rationale behind the removal of outer from ndarray-linalg?

Kronecker product looks like a good addition - I have added the good first issue label to the issue and linked it in #597. If you fancy submitting a PR we can work on it from there 👍

[notmgsk]

After some failed attempts I was able to come up with a working generic implementation of kron, starting from the snippet above (@ethanhs)

use ndarray::{arr2, Array2, LinalgScalar};
use num_complex::Complex64;

pub fn kron<T>(a: &Array2<T>, b: &Array2<T>) -> Array2<T>
where
    T: LinalgScalar,
{
    let dima = a.shape()[0];
    let dimb = b.shape()[0];
    let dimout = dima * dimb;
    let mut out = Array2::zeros((dimout, dimout));
    for (mut chunk, elem) in out.exact_chunks_mut((dimb, dimb)).into_iter().zip(a.iter()) {
        let v: Array2<T> = Array2::from_elem((dimb, dimb), *(elem)) * b;
        chunk.assign(&v);
    }
    out
}

fn main() {
    let a = arr2(&[[1.0, 0.0], [0.0, 1.0]]);
    let b = arr2(&[[0.0, 1.0], [1.0, 0.0]]);
    let c = kron(&a, &b);
    println!("{:#}", c);

    let d = arr2(&[
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
    ]);
    let e = arr2(&[
        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
    ]);
    let f = kron(&d, &e);
    println!("{:#}", f);
}

I'm not even a week in to learning Rust, so take the above with a grain of salt. And any critique is welcome.

[emmatyping]

Note that a Kronecker between non-square matrices will be a little different, but yes that looks correct for square matrices, which is what I think we both care about since we are using it for quantum computing purposes. For the more general kronecker it is kron(A, B) = C where A = x*y, B=m*n C= xm*yn of course.

The more general implementation can be made in a very similar fashion, you just need to get the second dimension of b. Do you want to make a PR? If not I am happy to.

[notmgsk]

Note that a Kronecker between non-square matrices will be a little different, but yes that looks correct for square matrices, which is what I think we both care about since we are using it for quantum computing purposes. For the more general kronecker it is kron(A, B) = C where A = xy, B=mn C= xm*yn of course.

The more general implementation can be made in a very similar fashion, you just need to get the second dimension of b. Do you want to make a PR? If not I am happy to.

Sometimes I forget non-square matrices exist :)

Sure, I'll give it a try.

[bluss]

Superseded by #1105