diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py
index 219cea714c1..df98806720c 100644
--- a/odl/operator/pspace_ops.py
+++ b/odl/operator/pspace_ops.py
@@ -15,6 +15,7 @@
 from odl.operator.operator import Operator
 from odl.operator.default_ops import ZeroOperator
 from odl.space import ProductSpace
+from odl.util import COOMatrix
 
 
 __all__ = ('ProductSpaceOperator',
@@ -171,7 +172,16 @@ def __init__(self, operators, domain=None, range=None):
             if not all(isinstance(op, Operator) for op in operators.data):
                 raise ValueError('sparse matrix `operator` contains non-'
                                  '`Operator` entries')
+            # scipy sparse matrix not supported (deprecated due to API changes)
+            # keep now for backward compatibility
+            print('Warning: scipy.sparse.spmatrix is deprecated.')
+            self.__ops = COOMatrix(operators.data,
+                                   (operators.row, operators.col),
+                                   operators.shape)
+
+        elif isinstance(operators, COOMatrix):
             self.__ops = operators
+
         else:
             self.__ops = self._convert_to_spmatrix(operators)
 
@@ -229,8 +239,9 @@ def __init__(self, operators, domain=None, range=None):
     @staticmethod
     def _convert_to_spmatrix(operators):
         """Convert an array-like object of operators to a sparse matrix."""
+
         # Lazy import to improve `import odl` time
-        import scipy.sparse
+        # import scipy.sparse
 
         # Convert ops to sparse representation. This is not trivial because
         # operators can be indexable themselves and give the wrong impression
@@ -279,8 +290,8 @@ def _convert_to_spmatrix(operators):
         # in `coo_matrix.__init__`
         data_arr = np.empty(len(data), dtype=object)
         data_arr[:] = data
-        return scipy.sparse.coo_matrix((data_arr, (irow, icol)),
-                                       shape=(nrows, ncols))
+
+        return COOMatrix(data_arr, (irow, icol), (nrows, ncols))
 
     @property
     def ops(self):
@@ -382,7 +393,7 @@ def derivative(self, x):
         data[:] = deriv_ops
         indices = [self.ops.row, self.ops.col]
         shape = self.ops.shape
-        deriv_matrix = scipy.sparse.coo_matrix((data, indices), shape)
+        deriv_matrix = COOMatrix(data, indices, shape)
         return ProductSpaceOperator(deriv_matrix, self.domain, self.range)
 
     @property
@@ -431,7 +442,7 @@ def adjoint(self):
         data[:] = adjoint_ops
         indices = [self.ops.col, self.ops.row]  # Swap col/row -> transpose
         shape = (self.ops.shape[1], self.ops.shape[0])
-        adj_matrix = scipy.sparse.coo_matrix((data, indices), shape)
+        adj_matrix = COOMatrix(data, indices, shape)
         return ProductSpaceOperator(adj_matrix, self.range, self.domain)
 
     def __getitem__(self, index):
@@ -731,6 +742,7 @@ class BroadcastOperator(Operator):
     ReductionOperator : Calculates sum of operator results.
     DiagonalOperator : Case where each operator should have its own argument.
     """
+
     def __init__(self, *operators):
         """Initialize a new instance
 
@@ -901,6 +913,7 @@ class ReductionOperator(Operator):
     DiagonalOperator : Case where each operator should have its own argument.
     SeparableSum : Corresponding construction for functionals.
     """
+
     def __init__(self, *operators):
         """Initialize a new instance.
 
@@ -1145,7 +1158,7 @@ def __init__(self, *operators, **kwargs):
         data = np.empty(n_ops, dtype=object)
         data[:] = operators
         shape = (n_ops, n_ops)
-        op_matrix = scipy.sparse.coo_matrix((data, (irow, icol)), shape)
+        op_matrix = COOMatrix(data, (irow, icol), shape)
         super(DiagonalOperator, self).__init__(op_matrix, **kwargs)
 
         self.__operators = tuple(operators)
diff --git a/odl/util/__init__.py b/odl/util/__init__.py
index c43c403e4b7..2f261bfab28 100644
--- a/odl/util/__init__.py
+++ b/odl/util/__init__.py
@@ -17,6 +17,7 @@
 from .testutils import *
 from .utility import *
 from .vectorization import *
+from .sparse import *
 
 __all__ = ()
 __all__ += graphics.__all__
@@ -26,3 +27,4 @@
 __all__ += testutils.__all__
 __all__ += utility.__all__
 __all__ += vectorization.__all__
+__all__ += sparse.__all__
diff --git a/odl/util/sparse.py b/odl/util/sparse.py
new file mode 100644
index 00000000000..85cb7dfeb1c
--- /dev/null
+++ b/odl/util/sparse.py
@@ -0,0 +1,54 @@
+
+"""
+Sparse matrix representation for creating product space operators.
+"""
+
+__all__ = ('COOMatrix',)
+
+
+class COOMatrix():
+    """
+    Custom coo matrix representation for creating product space operators.
+
+    The columns, rows and data are stored in separate lists such that
+    A[i[k], j[k]] = data[k].
+
+    Note that, the class is only used as a container and does not provide
+    any matrix operations. Further, no checks are performed on the data thus
+    duplicate and out-of-order indices are allowed and the user is responsible
+    for ensuring the correct shape of the matrix.
+
+    """
+
+    def __init__(self, data, index, shape):
+
+        # type check
+        if len(data) != len(index[0]) or len(data) != len(index[1]):
+            raise ValueError('data and index must have the same length')
+
+        self.__data = data
+        self.__index = index
+        self.__shape = shape
+
+    @property
+    def row(self):
+        """Return the row indices of the matrix."""
+        return self.__index[0]
+
+    @property
+    def col(self):
+        """Return the column indices of the matrix."""
+        return self.__index[1]
+
+    @property
+    def shape(self):
+        """Return the shape of the matrix."""
+        return self.__shape
+
+    @property
+    def data(self):
+        """Return the data of the matrix."""
+        return self.__data
+
+    def __repr__(self):
+        return f"COO matrix({self.data}, {self.index})"