Kronecker product needs a SparseMatrix implementation

[stared]

The current implementation of the Kronecker product #704 silently changes sparse matrices into dense matrices.

Example:

v1 = math.matrix([0, 2, 3, 0, 0], 'sparse')
v2 = math.matrix([0, 1, 0 -1], 'sparse')
v = math.kron(v1, v2)

Returns v, which is a DenseMatrix even though v1 and v2 areSparseMatrix. I would expect SparseMatrix instead.

Is it intended behavior?

[josdejong]

Good question. Not all functions do have support for SparseMatrices. If they don't have support, they convert the SparseMatrix input into a DenseMatrix, do their calculation, and return a DenseMatrix.

What you can do is wrap the function call so you get a SparseMatrix as result, like:

v = math.sparse(math.kron(v1, v2))

[mike239x]

As of right now you can just use something like this:

#!/usr/bin/env node
const { sparse } = require('mathjs');

function kron(m1, m2) {
    let re     = sparse([[0]]);
    let height = m1.size()[0] * m2.size()[0];
    let width  = m1.size()[1] * m2.size()[1];
    re.resize([height, width]);
    m1.forEach( (val1, pos1) => {
        m2.forEach( (val2, pos2) => {
            let i = pos1[0] * m2.size()[0] + pos2[0];
            let j = pos1[1] * m2.size()[1] + pos2[1];
            re.set([i, j], val1 * val2);
        }, skipZeros = true);
    }, skipZeros = true);
    return re;
}

I believe this is not the fastest way to do it, but it gets the job done.

Alternatively, I also wrote how to do that manually, but for this code to work you would have to fix #1166 first:

#!/usr/bin/env node
const { sparse } = require('mathjs');

function kron(m1, m2) {
    let _index    = [];
    let _values   = [];
    let _ptr      = [0];
    let nnz_num   = 0;
    let m1_height = m1.size()[0];
    let m1_width  = m1.size()[1];
    let m2_height = m2.size()[0];
    let m2_width  = m2.size()[1];
    let _size     = [m1_height * m2_height, m1_width * m2_width];
    // looping over columns of m1
    for (let m1_j = 0; m1_j < m1._ptr.length - 1; m1_j++) {
        let m1_from = m1._ptr[m1_j];
        let m1_to   = m1._ptr[m1_j+1];
        // looping over columns of m2
        for (let m2_j = 0; m2_j < m2._ptr.length - 1; m2_j++) {
            let m2_from = m2._ptr[m2_j];
            let m2_to   = m2._ptr[m2_j+1];
            // looping over elements of the current column of m1
            for (let m1_index = m1_from; m1_index < m1_to; m1_index++) {
                let m1_i   = m1._index[m1_index];
                let m1_val = m1._values[m1_index];
                // looping over elements of the current column of m2
                for (let m2_index = m2_from; m2_index < m2_to; m2_index++) {
                    let m2_i   = m2._index[m2_index];
                    let m2_val = m2._values[m2_index];
                    _index.push(m1_i * m2_height + m2_i);
                    _values.push(m1_val * m2_val);
                }
            }
            nnz_num += (m1_to - m1_from) * (m2_to - m2_from);
            _ptr.push(nnz_num);
        }
    }
    return /*sparse_matrix*/ {_values, _index, _ptr, _size};
}

// a small test:
let m = sparse([
    [0,2],
    [1,3],
]);
console.log(m);
let n = sparse([
    [1,3],
    [2,1],
]);
console.log(n);
let mn = sparse([
    [0,0,2,6],
    [0,0,4,2],
    [1,3,3,9],
    [2,1,6,3],
])
console.log(mn);
let re = kron(m,n);
console.log(re);

[josdejong]

Thanks for sharing Mike! Would you be interested in working out a SparseMatrix implementation of kron? (implementation, tests, docs)

[mike239x]

@josdejong if I manage my free time properly, then yes, but so far I am not good at this...

[josdejong]

:+1: no worries, only if you like it.