9

It seems numpy's einsum function does not work with scipy.sparse matrices. Are there alternatives to do the sorts of things einsum can do with sparse matrices?

In response to @eickenberg's answer: The particular einsum I'm wanting to is numpy.einsum("ki,kj->ij",A,A) - the sum of the outer products of the rows.

drevicko
  • 14,382
  • 15
  • 75
  • 97
  • 1
    Yes, there are alternatives, but there's no general "sparse einsum". It depends on what you need to do. – Fred Foo Apr 27 '14 at 13:54

2 Answers2

5

A restriction of scipy.sparse matrices is that they represent linear operators and are thus kept two dimensional, which leads to the question: Which operation are you seeking to do?

All einsum operations on a pair of 2D matrices are very easy to write without einsum using dot, transpose and pointwise operations, provided that the result does not exceed two dimensions.

So if you need a specific operation on a number of sparse matrices, it is probable that you can write it without einsum.

UPDATE: A specific way to implement np.einsum("ki, kj -> ij", A, A) is A.T.dot(A). In order to convince yourself, please try the following example:

import numpy as np
rng = np.random.RandomState(42)
a = rng.randn(3, 3)
b = rng.randn(3, 3)
the_einsum_ab = np.einsum("ki, kj -> ij", a, b)
the_a_transpose_times_b = a.T.dot(b)
# We write a test in order to assert equality
from numpy.testing import assert_array_equal
assert_array_equal(the_einsum_ab, the_a_transpose_times_b)  # This passes, so equality

This result is slightly more general. Now if you use b = a you obtain your specific result.

eickenberg
  • 14,152
  • 1
  • 48
  • 52
4

einsum translates the index string into a calculation using the C version of np.nditer. http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html is a nice introduction to nditer. Note especially the Cython example at the end.

https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py is a Python simulation of the einsum.

scipy.sparse has its own code (ultimately in C) to perform the basic operations, summation, matrix multiplication, etc. Sparse matricies have their own data structures. They can be lists, dictionaries, or a set of numpy arrays. Numpy notation can be used because sparse has the appropriate __xxx__ methods.

A sparse matrix is a matrix, a 2d array object. A sparse einsum could be written, but it would end up using the sparse matrix multiplication, not nditer. So at best it would be a notational convenience.

Sparse csr_matrix.dot is:

def dot(self, other):
    """Ordinary dot product
    ...
    """
    return self * other

A=sparse.csr_matrix([[1,2],[3,4]])
A.dot(A.T).A
(A*A.T).A
A.__rmul__(A.T).A
A.__mul__(A.T).A
np.einsum('ij,kj',A.A,A.A)
# array([[ 5, 11],
#        [11, 25]])
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • Can the "optimize=True" parameter in einsum somehow optimize sparse operations? I work a lot with np.einsum, in many dimensions where many of them are mostly sparse. – Bernardo Costa Apr 08 '21 at 12:47
  • 2
    @BernardoCosta, `einsum` does not deal with sparsity. It doesn't work with `scipy.sparse` matrices, and doesn't examine the dense arrays for sparsity either. Its `optimize` mode looks at the order of calculation when given more than 2 arrays. As the docs say, look at `einsum_path` for more details on that. `einsum` has changed since I first wrote this answer; it now has more in common with the `np.matmul/@`. function – hpaulj Apr 08 '21 at 15:35