[Scipy-tickets] [SciPy] #1859: Sparse matrix subtraction with another matrix seems to regard the other matrix as a scalar

SciPy Trac scipy-tickets@scipy....
Fri Mar 1 19:13:11 CST 2013


#1859: Sparse matrix subtraction with another matrix seems to regard the other
matrix as a scalar
------------------------------------+---------------------------------------
 Reporter:  photon.engine           |       Owner:  jakevdp    
     Type:  defect                  |      Status:  new        
 Priority:  normal                  |   Milestone:  Unscheduled
Component:  scipy.sparse            |     Version:  0.10.1     
 Keywords:  sparse scalar subtract  |  
------------------------------------+---------------------------------------
 The bulk of the problematic code below initializes a banded array for
 numerically solving the oscillator equation; it can be ignored if desired.
 The checks toward the end don't show any problems, but I couldn't subtract
 two matrices of the same shape on the last line; it seems SciPy treats
 "rhs" as a scalar.

 {{{
 """
 Solve the simple harmonic oscillator equation.
 """

 # Python 3.2.3, SciPy 0.10.1+dfsg1-4, NumPy 1:1.6.2-1ubuntu1, Ubuntu 12.10

 import numpy
 import scipy.linalg
 import scipy.sparse

 boundary = 1
 points = numpy.linspace(-boundary, boundary, 2001)
 step = points[1] - points[0]

 # left-hand side: second derivative
 side_diagonal = numpy.ones(len(points)) * (step ** -2)
 main_diagonal = numpy.ones(len(points)) * (-2 * step ** -2)
 # left-hand side: other terms
 main_diagonal += numpy.ones(len(points)) * 9
 # right-hand side
 rhs = numpy.zeros(len(points))

 # diagonals in the banded matrix
 diagonals = numpy.vstack([side_diagonal, main_diagonal, side_diagonal])
 # boundary conditions
 diagonals[0, 1] = 1 / step
 diagonals[1, 0] = -1 / step
 diagonals[-1, -2] = -1 / step
 diagonals[-2, -1] = 1 / step
 rhs[0] = 1
 rhs[-1] = 0

 # solve the equation
 solution = scipy.linalg.solve_banded((1, 1), diagonals, rhs)
 # check the solution
 matrix = scipy.sparse.spdiags(diagonals, (1, 0, -1), len(points),
 len(points))

 # checks
 print(matrix.shape, solution.shape, numpy.dot(matrix, solution).shape,
   rhs.shape)
 # (2001, 2001) (2001,) (2001,) (2001,)
 print(type(matrix), type(solution), type(numpy.dot(matrix, solution)),
   type(rhs))
 # <class 'scipy.sparse.dia.dia_matrix'> <class 'numpy.ndarray'> <class
 'numpy.ndarray'> <class 'numpy.ndarray'>
 import scipy.sparse.sputils
 print(scipy.sparse.sputils.isscalarlike(matrix),
   scipy.sparse.sputils.isscalarlike(solution),
   scipy.sparse.sputils.isscalarlike(numpy.dot(matrix, solution)),
   scipy.sparse.sputils.isscalarlike(rhs))
 # False False False False

 # File "/usr/lib/python3/dist-packages/scipy/sparse/compressed.py", line
 193
 # NotImplementedError: adding a scalar to a sparse matrix is not supported
 print(numpy.array(numpy.dot(matrix, solution)) - rhs)
 }}}

-- 
Ticket URL: <http://projects.scipy.org/scipy/ticket/1859>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list