more intuitive way to fill up three elements right of diagonal vectors

Question

I need to fill up the matrix in the figure. My way is to use the loop which is not east to read. When the dimension is large, it not efficient too. Is there any way to make it easier to read and probably quicker (say using vectorization or some function I don't know?)

# The code is in python, but I am open to R as well.    
    S=6
    p=[0.1,0.3,0.6] # the pi in the figure
    # The reason I use loop for p is that it handles a flexible dimension of p
    mat = np.zeros((S, S))
    p = np.array(p) 
    for i in range(S):
        for j, x in enumerate(p): 
                if i + j < S-1:
                    mat[i+j][i] = x

                elif i + j == S-1:
                    mat[S-1][i] = p[j:].sum()
                else:
                    pass

    mat.T

[This](http://stackoverflow.com/questions/5852495/how-do-i-find-the-scalar-product-of-a-numpy-matrix) will be helpful. Also, if you use `for`, change the values on yellow parts after `for`. No need to use `if` inside `for`, which is cumbersome for only two change. — Sangbok Lee, Apr 01 '17 at 15:36

score 1 · Answer 1 · answered Apr 01 '17 at 15:36

Here is an example - fill in the values you need (and adjust the final column)...

mat <- matrix(0,nrow=10,ncol=10)
diag(mat) <- 1
diag(mat[1:(ncol(mat)-1),-1]) <- 2
diag(mat[1:(ncol(mat)-2),-(1:2)]) <- 3

mat
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    2    3    0    0    0    0    0    0     0
 [2,]    0    1    2    3    0    0    0    0    0     0
 [3,]    0    0    1    2    3    0    0    0    0     0
 [4,]    0    0    0    1    2    3    0    0    0     0
 [5,]    0    0    0    0    1    2    3    0    0     0
 [6,]    0    0    0    0    0    1    2    3    0     0
 [7,]    0    0    0    0    0    0    1    2    3     0
 [8,]    0    0    0    0    0    0    0    1    2     3
 [9,]    0    0    0    0    0    0    0    0    1     2
[10,]    0    0    0    0    0    0    0    0    0     1

Thanks, decomposing it to three diagonals is smart. – MLE Apr 01 '17 at 15:41 — MLE, Apr 01 '17 at 15:41

DSM · Accepted Answer · 2017-04-01T16:32:09.187

As always, there are lots of equivalent ways to build a matrix. When I first saw it, I thought "oh, it's like the top part of a Toeplitz", but summing the diagonals in p should also work (note I've included the bottom-right change). Some googling revealed yet another way, using scipy.sparse.diags:

import numpy as np
import scipy.sparse
from scipy.linalg import triu, toeplitz

def build_toep(S, p):
    out = triu(toeplitz(p + [0]*(S-len(p))))
    out[-2:,-1] = [1 - p[1], 1]
    return out

def build_diag(S, p):
    out = sum(np.diag([v]*(S-i), i) for i,v in enumerate(p))
    out[-2:,-1] = [1 - p[1], 1]
    return out

def build_sparse(S, p):
    out = scipy.sparse.diags(p, range(len(p)), shape=(S, S)).toarray()
    out[-2:,-1] = [1 - p[1], 1]
    return out

which gives

In [150]: S, p = 6, [0.1, 0.3, 0.6]

In [151]: build_toep(S, p)
Out[151]: 
array([[ 0.1,  0.3,  0.6,  0. ,  0. ,  0. ],
       [ 0. ,  0.1,  0.3,  0.6,  0. ,  0. ],
       [ 0. ,  0. ,  0.1,  0.3,  0.6,  0. ],
       [ 0. ,  0. ,  0. ,  0.1,  0.3,  0.6],
       [ 0. ,  0. ,  0. ,  0. ,  0.1,  0.7],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  1. ]])

In [152]: np.allclose(build_toep(S, p), build_diag(S, p))
Out[152]: True

In [153]: np.allclose(build_toep(S, p), build_sparse(S, p))
Out[153]: True

It didn't occurred to me that this is part of toeplitz(never heard it before), thanks. — MLE, Apr 01 '17 at 16:08

more intuitive way to fill up three elements right of diagonal vectors

2 Answers2