Creating a special matrix in numpy

Question

[a b c       ] 
[  a b c     ]
[    a b c   ]
[      a b c ]

Hello

For my economics course we are suppose to create an array that looks like this. The problem is I am an economist not a programmer. We are using numpy in python. Our professor says college is not preparing us for the real world and wants us to learn programming (which is a good thing). We are not allowed to use any packages and must come up with an original code. Does anybody out there have any idea how to make this matrix. I have spent hours trying codes and browsing the internet looking for help and have been unsuccessful.

Of what type are the blank elements? And what does this have to do with economics? — xavier, Sep 15 '14 at 02:12
the blanks are 0. We are writing our own hodrick prescott filter — thegodfather, Sep 15 '14 at 02:17

chthonicdaemon · Answer 1 · 2014-09-15T15:38:37.777

7

This kind of matrix is called a Toeplitz matrix or constant diagonal matrix. Knowing this leads you to scipy.linalg.toeplitz:

import scipy.linalg
scipy.linalg.toeplitz([1, 0, 0, 0], [1, 2, 3, 0, 0, 0])

=>

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

edited Sep 15 '14 at 15:38

answered Sep 15 '14 at 15:28

chthonicdaemon

19,180
2
52
66

tom10 · Answer 2 · 2014-09-15T15:59:58.780

The method below fills one diagonal at a time:

import numpy as np
x = np.zeros((4, 6), dtype=np.int)
for i, v in enumerate((6,7,8)):
    np.fill_diagonal(x[:,i:], v)

array([[6, 7, 8, 0, 0, 0],
       [0, 6, 7, 8, 0, 0],
       [0, 0, 6, 7, 8, 0],
       [0, 0, 0, 6, 7, 8]])

or you could do the one liner:

x = [6,7,8,0,0,0]
y = np.vstack([np.roll(x,i) for i in range(4)])

Personally, I prefer the first since it's easier to understand and probably faster since it doesn't build all the temporary 1D arrays.

Edit:
Since a discussion of efficiency has come up, it might be worthwhile to run a test. I also included time to the toeplitz method suggested by chthonicdaemon (although personally I interpreted the question to exclude this approach since it uses a package rather than using original code -- also though speed isn't the point of the original question either).

import numpy as np
import timeit
import scipy.linalg as sl

def a(m, n):    
    x = np.zeros((m, m), dtype=np.int)
    for i, v in enumerate((6,7,8)):
        np.fill_diagonal(x[:,i:], v)

def b(m, n):
    x = np.zeros((n,))
    x[:3] = vals
    y = np.vstack([np.roll(x,i) for i in range(m)])

def c(m, n):
    x = np.zeros((n,))
    x[:3] = vals
    y = np.zeros((m,))
    y[0] = vals[0]
    r = sl.toeplitz(y, x)
    return r

m, n = 4, 6
print timeit.timeit("a(m,n)", "from __main__ import np, a, b, m, n", number=1000)
print timeit.timeit("b(m,n)", "from __main__ import np, a, b, m, n", number=1000)
print timeit.timeit("c(m,n)", "from __main__ import np, c, sl, m, n", number=1000)

m, n = 1000, 1006
print timeit.timeit("a(m,n)", "from __main__ import np, a, b, m, n", number=1000)
print timeit.timeit("b(m,n)", "from __main__ import np, a, b, m, n", number=1000)
print timeit.timeit("c(m,n)", "from __main__ import np, c, sl, m, n", number=100)

# which gives:
0.03525209  # fill_diagonal
0.07554483  # vstack
0.07058787  # toeplitz

0.18803215  # fill_diagonal
2.58780789  # vstack
1.57608604  # toeplitz

So the first method is about a 2-3x faster for small arrays and 10-20x faster for larger arrays.

I like the one liner. I don't think it would be any slower since the array slicing is expensive too. — xavier, Sep 15 '14 at 02:43
@xavier: slicing is cheap because it's just a different view on the data and not a copy, though it's hard to say which will win for such small arrays (here, for example, the header for the view could be larger than the data). http://docs.scipy.org/doc/numpy/glossary.html#term-view — tom10, Sep 15 '14 at 13:38
I would be hesitant to even add the one liner as it reinforces poor numpy code. As shown `np.roll` has to allocate two new 1d arrays for each row, allocate a 2d array for the vstack, and finally copy all of the 1d arrays over the 2d array. — Daniel, Sep 15 '14 at 13:58

score 2 · Answer 3 · edited May 23 '17 at 12:00

2

This is a simplified tridiagonal matrix. So it is essentially a this question

def tridiag(a, b, c, k1=-1, k2=0, k3=1):
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

a = [1, 1]; b = [2, 2, 2]; c = [3, 3]
A = tridiag(a, b, c)
print(A)

Result:

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

edited May 23 '17 at 12:00

Community

1
1

answered Sep 15 '14 at 15:48

BKay

1,397
1
15
26

score 1 · Answer 4 · answered Sep 15 '14 at 02:31

Something along the lines of

import numpy as np
def createArray(theinput,rotations) :
    l = [theinput]
    for i in range(1,rotations) :
        l.append(l[i-1][:])
        l[i].insert(0,l[i].pop())
    return np.array(l)

print(createArray([1,2,3,0,0,0],4))
"""
[[1 2 3 0 0 0]
 [0 1 2 3 0 0]
 [0 0 1 2 3 0]
 [0 0 0 1 2 3]]
"""

Eelco Hoogendoorn · Answer 5 · 2014-09-15T22:57:51.373

If you care about efficiency, it is hard to beat this:

import numpy as np

def create_matrix(diags, n):
    diags = np.asarray(diags)
    m = np.zeros((n,n+len(diags)-1), diags.dtype)
    s = m.strides
    v = np.lib.index_tricks.as_strided(
        m,
        (len(diags),n),
        (s[1],sum(s)))
    v[:] = diags[:,None]
    return m

print create_matrix(['a','b','c'], 8)

Might be a little over your head, but then again that's good inspiration ;)

Or even better: a solution which has both O(n) storage and runtime requirements, rather than all the other solutions posted thus far, which are O(n^2)

import numpy as np

def create_matrix(diags, n):
    diags = np.asarray(diags)
    b = np.zeros(len(diags)+n*2, diags.dtype)
    b[n:][:len(diags)] = diags
    s = b.strides[0]
    v = np.lib.index_tricks.as_strided(
        b[n:],
        (n,n+len(diags)-1),
        (-s,s))
    return v

print create_matrix(np.arange(1,4), 8)

Sigve Karolius · Answer 6 · 2015-11-17T19:57:44.990

This is an old question, however some new input can always be useful.

I create tridiagonal matrices in python using list comprehension.

Say a matrix that is symmetric around "-2" and has a "1" on either side:

           -2   1   0
Tsym(3) =>  1  -2   1 
            0   1  -2

This can be created using the following "one liner":

Tsym = lambda n: [ [ 1 if (i+1==j or i-1==j) else -2 if j==i else 0 for i in xrange(n) ] for j in xrange(n)] # Symmetric tridiagonal matrix (1,-2,1)

A different case (that several of the other people answering has solved perfectly fine) is:

             1   2   3   0   0   0 
Tgen(4,6) => 0   1   2   3   0   0
             0   0   1   2   3   0
             0   0   0   1   2   3

Can be made using the one liner shown below.

Tgen = lambda n,m: [ [ 1 if i==j else 2 if i==j+1 else 3 if i==j+2 else 0 for i in xrange(m) ] for j in xrange(n)] # General tridiagonal matrix  (1,2,3)

Feel free to modify to suit your specific needs. These matrices are very common when modelling physical systems and I hope this is useful to someone (other than me).

score 0 · Answer 7 · answered Feb 04 '18 at 13:16

Hello since your professor asked you not to import any external package, while most answers use numpy or scipy. You better use only python List to create 2D array (compound list), then populate its diagonals with the items you wish, Find the code below

def create_matrix(rows = 4, cols = 6):
    mat = [[0 for col in range(cols)] for row in range(rows)] # create a mtrix filled with zeros of size(4,6) 
    for row in range(len(mat)): # gives number of lists in the main list, 
        for col in range(len(mat[0])): # gives number of items in sub-list 0, but all sublists have the same length
            if row == col:
                mat[row][col] = "a"
            if col == row+1:
                mat[row][col] = "b"
            if col == row+2:
                mat[row][col] = "c"
    return mat

create_matrix(4, 6)

[['a', 'b', 'c', 0, 0, 0],

[0, 'a', 'b', 'c', 0, 0],

[0, 0, 'a', 'b', 'c', 0],

[0, 0, 0, 'a', 'b', 'c']]

score 0 · Answer 8 · answered Aug 04 '21 at 00:08

Creating Band Matrix

Check out the definition for it in wiki : https://en.wikipedia.org/wiki/Band_matrix

You can use this function to create band matrices like diagonal matrix with offset=1 or tridiagonal matrix (The one you are asking about) with offset=1 or Pentadiagonal Matrix with offset=2

def band(size=10, ones=False, low=0, high=100, offset=2):
    shape = (size, size)
    n_matrix = np.random.randint(low, high, shape) if not ones else np.ones(shape,dtype=int)
    n_matrix = np.triu(n_matrix, -1*offset)
    n_matrix = np.tril(n_matrix, offset)
    return n_matrix

In your case you should use this

rand_tridiagonal = band(size=6,offset=1)
print(rand_tridiagonal)

Creating a special matrix in numpy

8 Answers8

Linked