New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

COO.__matmul__ gives the wrong shape of resultant array #202

Closed
gongliyu opened this Issue Oct 22, 2018 · 6 comments

Comments

Projects
None yet
2 participants
@gongliyu
Copy link
Contributor

gongliyu commented Oct 22, 2018

Please see the following, COO.matmul is not consistent with np.ndarray.matmul

In [1]: import numpy as np

In [2]: import sparse

In [3]: npx = np.random.random((1,2,3,4))

In [4]: npy = np.random.random((1,2,4,3))

In [5]: npz = npx @ npy

In [6]: npz.shape
Out[6]: (1, 2, 3, 3)

In [7]: spx = sparse.random((1,2,3,4), 0.8)

In [8]: spy = sparse.random((1,2,4,3), 0.8)

In [9]: spz = spx @ spy

In [10]: type(spz)
Out[10]: numpy.ndarray

In [11]: spz.shape
Out[11]: (1, 2, 3, 1, 2, 3)
@hameerabbasi

This comment has been minimized.

Copy link
Collaborator

hameerabbasi commented Oct 22, 2018

Unfortunately, the fix to this will have to wait until a CSR replacement is in, as I can see no easy way out.

@gongliyu

This comment has been minimized.

Copy link
Contributor

gongliyu commented Oct 22, 2018

Can we utilize the scipy.sparse matmul functionality to implement it?
Please see the following workaround,

def _matmul_2_2(a, b):
    """
    """
    assert a.ndim == 2 and b.ndim == 2
    if has_pydata_sparse():
        if is_pydata_sparse(a):
            a = as_scipy_coo(a).tocsr()
        if is_pydata_sparse(b):
            b = as_scipy_coo(b).tocsc()
    res = a @ b
    if has_pydata_sparse() and (is_pydata_sparse(a) or is_pydata_sparse(b)):
        res = as_pydata_coo(res)
    return res
    
    
def matmul(a, b, *, use_recursion=True):
    """
    """
    if a.ndim < 2:
        raise ValueError('a.ndim < 2')
    if b.ndim < 2:
        raise ValueError('b.ndim < 2')

    if a.shape[-1] != b.shape[-2]:
        raise ValueError('a.shape[-1] != b.shape[-2]')

    if a.ndim == 2 and b.ndim == 2:
        return _matmul_2_2(a, b)

    if not is_scipy_sparse(a):
        while a.ndim < b.ndim:
            a = a[np.newaxis]

    if not is_scipy_sparse(b):
        while a.ndim > b.ndim:
            b = b[np.newaxis]

    for i, j in zip(a.shape[:-2], b.shape[:-2]):
        if i != 1 and j != 1 and i != j:
            raise ValueError('shapes of a and b are not broadcastable')

    if use_recursion:
        def _matmul(a, b):
            if a.ndim == 2:
                assert b.ndim == 2
                return _matmul_2_2(a, b)
            res = []
            length = 0
            if not is_scipy_sparse(a):
                length = max(length, a.shape[0])
            if not is_scipy_sparse(b):
                length = max(length, b.shape[0])
            for i in range(length):
                a_next = (a if is_scipy_sparse(a) else a[0]
                          if a.shape[0] == 1 else a[i])
                b_next = (b if is_scipy_sparse(b) else b[0]
                          if b.shape[0] == 1 else b[i])
                res.append(_matmul(a_next, b_next))
            return base_ops.stack(res)
        return _matmul(a,b)
    else:
        if not (is_scipy_sparse(a) or is_scipy_sparse(b)):
            for i in range(a.ndim-2):
                if a.shape[i] == 1 and b.shape[i] != 1:
                    a = concatenate([a]*b.shape[i], axis=i)
                if a.shape[i] != 1 and b.shape[i] == 1:
                    b = concatenate([b]*a.shape[i], axis=i)

        shape1 = a.shape[:-2]
        if not is_scipy_sparse(a):
            a = a.reshape([np.prod(shape1)]+list(a.shape[-2:]))
        if not is_scipy_sparse(b):
            b = b.reshape([np.prod(shape1)]+list(b.shape[-2:]))
        res = []
        for i in range(a.shape[0]):
            a_i = a if is_scipy_sparse(a) else a[i]
            b_i = b if is_scipy_sparse(b) else b[i]
            tmp = _matmul_2_2(a_i, b_i)
            tmp = as_pydata_coo(tmp)
            res.append(tmp)
        return base_ops.stack(res).reshape(shape1+(a.shape[-2],b.shape[-1]))
@hameerabbasi

This comment has been minimized.

Copy link
Collaborator

hameerabbasi commented Oct 22, 2018

Well, we're looping in Python at the end of this... It'll be slow but it works. If you're willing to make that into a PR we'll accept it. The definition is in sparse/coo/common.py:dot. You can use scipy.sparse.csr_matrix.dot along with the sparse.COO.tocsr to implement the 2-D version.

hameerabbasi added a commit that referenced this issue Oct 23, 2018

Fix sparse.dot returning incorrect shape for N-D arrays. (#204)
* fix #202

* fix 202 COO.__matmul__ returns wrong shaped array

* try fix some syntax error on python 2.7

* fix some travis-ci check error

* fix a CI check error

* add functionality to squeeze a and b before matmul

* improve coverage

* remove some extra lines

* add docs for matmul
@gongliyu

This comment has been minimized.

Copy link
Contributor

gongliyu commented Dec 11, 2018

There is still an bug remaining for this issue:

In [1]: import numpy as np

In [2]: import sparse

In [3]: npx = np.random.random((3,4))

In [4]: npy = np.random.random((1,2,4,3))

In [5]: np.matmul(npx, npy).shape
Out[5]: (1, 2, 3, 3)

In [6]: spx = sparse.random((3,4), 0.8)

In [7]: spy = sparse.random((1,2,4,3), 0.8)

In [8]: sparse.matmul(spx, spy).shape
Out[8]: (3, 1, 2, 3)

It is caused by 187-165 in common.py

 # When one of the input is less than 2-d, it is equivalent to dot
    if a.ndim <= 2 or b.ndim <= 2:
        return dot(a, b)

The piece of code was intend to speed up using dot. However, in this situation, matmul is not equivalent to dot. So we should reshape the result after dot.

@hameerabbasi Would you please reopen this issue?

Thanks!

@hameerabbasi hameerabbasi reopened this Dec 11, 2018

@hameerabbasi

This comment has been minimized.

Copy link
Collaborator

hameerabbasi commented Dec 11, 2018

Would you like to make a PR or should I?

@gongliyu

This comment has been minimized.

Copy link
Contributor

gongliyu commented Dec 11, 2018

I will make a PR. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment