Generating banded matrices using numpy

Question

I'm using the following piece of code to create a banded matrix from a generator g:

def banded(g, N):
    """Creates a `g` generated banded matrix with 'N' rows"""
    n=len(g)
    T = np.zeros((N,N+n-1))
    for x in range(N):
        T[x][x:x+n]=g
    return T

The usage is simple as:

banded([1,2,3], 3)

and it returns

[1, 2, 3, 0, 0]
[0, 1, 2, 3, 0]
[0, 0, 1, 2, 3]

It will used mostly to solve a finite difference model with a given stencil for example (-1, 1)

Is there a better way to generate that stencil? I could not find any good NumPy Function for that.

By better I mean, faster, using less memory, removing the loop from python and sending to Numpy stack. Any (or all) of those are improvements.

Please define "better." Shorter code? More general (in what way)? Faster? Consuming less memory? Able to work in higher dimensions? — John Zwinck, Sep 23 '18 at 07:38
Faster, consuming less memory, removing the loop from python and sending to numpy stack. Any (or all) of those are better. — Lin, Sep 23 '18 at 07:38
OK. What's a realistic input in the real program? How long is `g` typically and how large is `N`? — John Zwinck, Sep 23 '18 at 07:40
`g` is not very large, varying from 2 to 32 elements (is a regular finite difference stencil). But `N` can be slightly large. Like `100000`, anythow `N` that really dictates the size of matrix. And the size grows `N**2`. — Lin, Sep 23 '18 at 07:57

Divakar · Accepted Answer · 2018-09-23T08:24:55.120

Here's one with np.lib.stride_tricks.as_strided to give us a 2D view into a zeros padded 1D version of the input and as such pretty memory efficient and hence performant too. This trick had been explored numerous times - 1,2.

Thus, the implementation would be -

def sliding_windows(a, W):
    a = np.asarray(a)
    p = np.zeros(W-1,dtype=a.dtype)
    b = np.concatenate((p,a,p))
    s = b.strides[0]
    strided = np.lib.stride_tricks.as_strided
    return strided(b[W-1:], shape=(W,len(a)+W-1), strides=(-s,s))

Sample runs -

In [99]: a = [1,2,3]

In [100]: sliding_windows(a, W=3)
Out[100]: 
array([[1, 2, 3, 0, 0],
       [0, 1, 2, 3, 0],
       [0, 0, 1, 2, 3]])

In [101]: a = [1,2,3,4,5]

In [102]: sliding_windows(a, W=3)
Out[102]: 
array([[1, 2, 3, 4, 5, 0, 0],
       [0, 1, 2, 3, 4, 5, 0],
       [0, 0, 1, 2, 3, 4, 5]])

With the same philosophy, but less messier version, we can also leverage np.lib.stride_tricks.as_strided based scikit-image's view_as_windows to get sliding windows. More info on use of as_strided based view_as_windows.

from skimage.util.shape import view_as_windows

def sliding_windows_vw(a, W):
    a = np.asarray(a)
    p = np.zeros(W-1,dtype=a.dtype)
    b = np.concatenate((p,a,p))
    return view_as_windows(b,len(a)+W-1)[::-1]

WOW! This method is the fastest by 4 orders of magnitude. – Lin Sep 23 '18 at 08:03 — Lin, Sep 23 '18 at 08:03

score 3 · Answer 2 · answered Sep 23 '18 at 07:40

3

You could use scipy.sparse.diags:

Input:

diags([1, 2, 3], [0, 1, 2], shape=(3,5)).toarray()

Output:

array([[ 1.,  2.,  3.,  0.,  0.],
      [ 0.,  1.,  2.,  3.,  0.],
      [ 0.,  0.,  1.,  2.,  3.]])

The second list, [0, 1, 2], is an offset list. It tells how offset from the diagonal you want a certain element to be.

answered Sep 23 '18 at 07:40

Joe Patten

1,664
1
9
15

Joe, yes, scipy sparse is one idea. A good one indeed. In fact the second parameter can be a `range`. – Lin Sep 23 '18 at 07:42
For sake of comparison I tried to generate the same matrix using the sparse method and the method that I used on example, the sparse matrix was just slightly faster in matrix allocation (for a 100,000x100,000) matrix. Of course I believe that sparseness will payoff when trying to solve the inverse. – Lin Sep 23 '18 at 07:59
Just for sake of information, I have used `scipy.sparse.diags(g,range(len(g), shape=(N, N+len(g)-1)).toarray()`, where `N` and `g` are the same from the question. – Lin Sep 23 '18 at 08:00

score 3 · Answer 3 · answered Apr 21 '20 at 17:48

You can also use scipy's toeplitz function, which is very similar to its matlab counterpart. It is also smart around the shape, do do not worry about the it.

import scipy.linalg as scl

# define first column and first line
column1 = [1,0,0]
line1 = [1,2,3,0,0]

scl.toeplitz(column1, line1)

If you want variable sizes, use your size parameter (N) to add zeros to the columns and lines dynamically. The following is my implementation for finite differences with a minimum size of 3.

col = [1,0,0] + [0] * (size -3)
lin = [1,2,1] + [0] * (size -3)
m = scl.toeplitz(col, lin)

Output:

array([[1., 2., 1., 0., 0.],
        [0., 1., 2., 1., 0.],
        [0., 0., 1., 2., 1.],
        [0., 0., 0., 1., 2.],
        [0., 0., 0., 0., 1.]])

Generating banded matrices using numpy

3 Answers3

Linked

Related