[FIXED] allocate memory in python for large scipy.sparse matrix operations

Issue

Is there a way I can allocate memory for scipy sparse matrix functions to process large datasets?

Specifically, I'm attempting to use Asymmetric Least Squares Smoothing (translated into python here and the original here) to perform a baseline correction on a large mass spec dataset (length of ~60,000).

The function (see below) uses the scipy.sparse matrix operations.

def baseline_als(y, lam, p, niter):
  L = len(y)
  D = sparse.csc_matrix(np.diff(np.eye(L), 2))
  w = np.ones(L)
  for i in xrange(niter):
    W = sparse.spdiags(w, 0, L, L)
    Z = W + lam * D.dot(D.transpose())
    z = spsolve(Z, w*y)
    w = p * (y > z) + (1-p) * (y < z)
  return z

I have no problem when I pass data sets that are 10,000 or less in length:

baseline_als(np.ones(10000),100,0.1,10)

But when passing larger data sets, e.g.

baseline_als(np.ones(50000), 100, 0.1, 10)

I get a MemoryError, for the line

  D = sparse.csc_matrix(np.diff(np.eye(L), 2))

Solution

Try changing

D = sparse.csc_matrix(np.diff(np.eye(L), 2))

to

diag = np.ones(L - 2)
D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2)

D will be a sparse matrix in DIAgonal format. If it turns out that being in CSC format is important, convert it using the tocsc() method:

D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2).tocsc()

The following example shows that the old and new versions generate the same matrix:

In [67]: from scipy import sparse

In [68]: L = 8

Original:

In [69]: D = sparse.csc_matrix(np.diff(np.eye(L), 2))

In [70]: D.A
Out[70]: 
array([[ 1.,  0.,  0.,  0.,  0.,  0.],
       [-2.,  1.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  1., -2.],
       [ 0.,  0.,  0.,  0.,  0.,  1.]])

New version:

In [71]: diag = np.ones(L - 2)

In [72]: D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2)

In [73]: D.A
Out[73]: 
array([[ 1.,  0.,  0.,  0.,  0.,  0.],
       [-2.,  1.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  1., -2.],
       [ 0.,  0.,  0.,  0.,  0.,  1.]])

Answered By - Warren Weckesser