Multiplication of sparse 3D matrix with 2D matrix

Question

I have sparse COO 3D A matrix NxNxM (The third dimension is M) and dense 2D B matrix (NxN), which will be an optimization variable for my solver. I want to multiply A and B to obtain (NxNXM) matrix. Then I need to sum the elements of the resulting 3D matrix to create 2D matrix (NxN).

For clarification, If I wouldn’t use a sparse matrix, I could do it with the below code. It is a multiplication of 3D and 2D matrices.

np.sum(np.einsum('ijk,jk->ijk', A, B))

where A is 3D matrix and B is 2D matrix.

My objective function is to minimize the sum of the elements of the multiplication of sparse A with decision variable matrix B.

How can I do this?

I am using pydata sparse library to create sparse matrix.

What does pydata provide? With the scipy.sparse package you have to use its own multiplication methods. numpy functions do it wrong. — hpaulj
You can only multiply 2D matrices with scipy.sparse but my A matrix is 3D and it is sparse. — Emin

CJR CJR · Accepted Answer · 2020-10-26T17:48:42

https://sparse.pydata.org/en/stable/

Fortunately in many cases we can leverage the existing scipy.sparse algorithms if we can intelligently transpose and reshape our multi-dimensional array into an appropriate 2-d sparse matrix, perform a modified sparse matrix operation, and then reshape and transpose back. These reshape and transpose operations can all be done at numpy speeds by modifying the arrays of coordinates. After scipy.sparse runs its operations (often written in C) then we can convert back to using the same path of reshapings and transpositions in reverse.

All it's doing is this; you should just cut out the middleman. You really aren't gaining anything by doing this as a 3d array, it's just harder to work with.

from scipy import sparse
import numpy as np    

M, N = 1000, 100

A = sparse.random(M, N ** 2, density=0.01, format='csr')
B = np.random.rand(N ** 2)

C = A.multiply(B).sum(axis=0).A

>>> A
<1000x10000 sparse matrix of type '<class 'numpy.float64'>'
    with 100000 stored elements in Compressed Sparse Row format>
>>> B.shape
(10000,)
>>> C.shape
(1, 10000)

Multiplication of sparse 3D matrix with 2D matrix

1 Answers