Working with banded matrices in Python

AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.

These notebooks are all available on Github.

Illustration of a few Python tools for banded matrices.

For more information, see the documentation:

In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
from scipy.sparse import diags

# import a module with a new local name for brevity:
import scipy.sparse.linalg as sp_linalg

You can create a simple diagonal matrix using the numpy diag command, or with eye for the identity matrix, e.g.

In [4]:
A = diag([1,2,3,4.])
A
Out[4]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  2.,  0.,  0.],
       [ 0.,  0.,  3.,  0.],
       [ 0.,  0.,  0.,  4.]])
In [5]:
inv(A)
Out[5]:
array([[ 1.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.5       ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.33333333,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.25      ]])
In [6]:
norm(A, inf)
Out[6]:
4.0

The numpy inv and norm functions work on this as a full matrix (with many elements equal to 0).

sparse diagonals

We can form a banded matrix using sparse storage as follows. See the documentation for more info.

In [7]:
diagonals = [[1, 2, 3, 4], [10, 20, 30], [100, 200, 300]]
A = diags(diagonals, offsets=[0, -1, 2], format='csc')

If you try to print A this way, it tells you what sort of object it is:

In [8]:
A
Out[8]:
<4x4 sparse matrix of type '<class 'numpy.float64'>'
	with 9 stored elements in Compressed Sparse Column format>

Inside a print function, it converts to a telling you the nonzeros:

In [9]:
print(A)

  (0, 0)	1.0
  (1, 0)	10.0
  (1, 1)	2.0
  (2, 1)	20.0
  (2, 2)	3.0
  (3, 2)	30.0
  (0, 2)	100.0
  (3, 3)	4.0
  (1, 3)	200.0

The shape attribute is also defined:

In [10]:
A.shape
Out[10]:
(4, 4)

To see it as the full matrix, you can convert it into an ordinary numpy array with the toarray method:

In [11]:
B = A.toarray()
B
Out[11]:
array([[   1.,    0.,  100.,    0.],
       [  10.,    2.,    0.,  200.],
       [   0.,   20.,    3.,    0.],
       [   0.,    0.,   30.,    4.]])

You could compute norm(B) with the numpy norm command, but trying to compute norm(A) gives an error. Instead you have to use the version of norm from scipy.sparse.linalg:

In [12]:
sp_linalg.norm(A,inf)
Out[12]:
212.0

Similarly to compute the inverse:

In [13]:
Ainv = sp_linalg.inv(A)
Ainv
Out[13]:
<4x4 sparse matrix of type '<class 'numpy.float64'>'
	with 16 stored elements in Compressed Sparse Column format>

Note that this returns a "sparse matrix" even though the inverse is in fact dense:

In [14]:
Ainv.toarray()
Out[14]:
array([[  6.00047994e-01,   3.99952006e-02,  -3.99952006e-03,
         -1.99976003e+00],
       [ -5.99928009e-04,   5.99928009e-05,   4.99940007e-02,
         -2.99964004e-03],
       [  3.99952006e-03,  -3.99952006e-04,   3.99952006e-05,
          1.99976003e-02],
       [ -2.99964004e-02,   2.99964004e-03,  -2.99964004e-04,
          1.00017998e-01]])

We get the same thing by computing the ordinary numpy inv of the dense array B:

In [15]:
inv(B)
Out[15]:
array([[  6.00047994e-01,   3.99952006e-02,  -3.99952006e-03,
         -1.99976003e+00],
       [ -5.99928009e-04,   5.99928009e-05,   4.99940007e-02,
         -2.99964004e-03],
       [  3.99952006e-03,  -3.99952006e-04,   3.99952006e-05,
          1.99976003e-02],
       [ -2.99964004e-02,   2.99964004e-03,  -2.99964004e-04,
          1.00017998e-01]])

Modifying elements

You can modify an element that is already nonzero in the sparse A:

In [16]:
A[0,0] = 5.
A.toarray()
Out[16]:
array([[   5.,    0.,  100.,    0.],
       [  10.,    2.,    0.,  200.],
       [   0.,   20.,    3.,    0.],
       [   0.,    0.,   30.,    4.]])

You can also modify an element that was 0., but this will require the internal representation to change so a warning message is printed. This is no big deal if we're just modifying a few elements, e.g. to impose boundary conditions, but if you are constructing a very large sparse matrix for a finite element problem one element at a time, for example, then this can be slow.

In [17]:
A[3,0] = -50.
A.toarray()

/Users/rjl/miniconda/envs/rbook2/lib/python3.6/site-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)
Out[17]:
array([[   5.,    0.,  100.,    0.],
       [  10.,    2.,    0.,  200.],
       [   0.,   20.,    3.,    0.],
       [ -50.,    0.,   30.,    4.]])

You can suppress printing all warnings with the command below. See the logging documentation for more info.

In [18]:
import logging
logging.captureWarnings(True)

matrix formats

Above we specified format='csc' in creating A. If you don't specify a format it uses dia.

In [19]:
diagonals = [[1, 2, 3, 4], [10, 20, 30], [100, 200, 300]]
A2 = diags(diagonals, offsets=[0, -1, 2])
In [20]:
type(A2), type(A)
Out[20]:
(scipy.sparse.dia.dia_matrix, scipy.sparse.csc.csc_matrix)

You can still apply many operations to it, e.g. compute the norm:

In [21]:
sp_linalg.norm(A2, inf)
Out[21]:
212.0

But this format does not allow you to add new nonzeros, this would give an error:

In [22]:
#A2[3,0] = -5.