CODEWITHSUNDEEP
pythonpython-3.xnumpy
Aenaon
Aenaon

Reputation: 3573

2D sparse matrix multiplied by 3D array

Is there any way to multiply a 2D sparse matrix by a 3D numpy array please? For example I have this function 

def myFun(x, p):
    r = 2
    out = x * np.log(p) + r * np.log(1-p)
    return out

where x is an array of dimension 3500, 90 and p another array with dimensions 3500, 90, 70. At the moment both x and p are dense arrays and I am just broadcasting when I call the function:

out = myFun(x[..., None], p)

However array x is quite sparse, only 7% of its elements are non-zero. On the other side, p doesnt have any zeros, only floats between zero and one. I am hoping though that with a sparse matrix (from scipy.sparse probably) I will see a speed improvement. However, I do not know how to do this operation or if this more efficient please.

I am using python 3.

Many thanks

Upvotes: 1

Views: 184

Answers (2)

max9111
max9111

Reputation: 6482

You can try the following implementation. For this simple function this looks like a bit exaggerated, but I also had troubles to get numexpr to work with Intel SVML (otherwise I would prefer numexpr). This solution should give 0.07s per call on a Quadcore i7 and should scale quite well on more cores. Please also note that the first call has a compilation overhead of about 0.5s.

Installing Intel SVML

import numpy as np
import numba as nb

x = np.random.uniform(-13, 1, (3500, 90, 1)).clip(0, None)
p = np.random.random((3500, 90, 70))

@nb.njit(parallel=True,fastmath=True)
def nb_myFun_sp(x, p):
    out=np.empty(p.shape,p.dtype)
    r = 2.

    for i in nb.prange(p.shape[0]):
      for j in range(p.shape[1]):
        if x[i,j,0]!=0.:
          x_=x[i,j,0]
          for k in range(p.shape[2]):
            out[i,j,k] = x_ * np.log(p[i,j,k]) + r * np.log(1.-p[i,j,k])
        else:
          for k in range(p.shape[2]):
            out[i,j,k] = r * np.log(1.-p[i,j,k])

    return out

@nb.njit(parallel=True,fastmath=True)
def nb_myFun(x, p):
    out=np.empty(p.shape,p.dtype)
    r = 2.

    for i in nb.prange(p.shape[0]):
      for j in range(p.shape[1]):
        x_=x[i,j,0]
        for k in range(p.shape[2]):
          out[i,j,k] = x_ * np.log(p[i,j,k]) + r * np.log(1.-p[i,j,k])
    return out

Upvotes: 1

Paul Panzer
Paul Panzer

Reputation: 53029

You can exploit the sparseness of x using the where keyword.

def sprse(x, p):
    r = 2
    out = x * np.log(p, where=x.astype(bool)) + r * np.log(1-p)
    return out

from timeit import timeit
x = np.random.uniform(-13, 1, (3500, 90, 1)).clip(0, None)
p = np.random.random((3500, 90, 70))
assert np.all(sprse(x, p)==myFun(x, p))
def f():
    return myFun(x, p)
print(timeit(f, number=3))
def f():
    return sprse(x, p)
print(timeit(f, number=3))

Sample run:

5.171174691990018
3.2122434769989923

Upvotes: 1

Related Questions