4

I want to create a toeplitz matrix of toeplitz matrix. H1, H2 and H3 are toeplitz matrices already. My result should look like that: H1 0 0 H2 H1 0 H3 H2 H1 0 H3 H2 0 0 H3

The existing toeplitz-function only accepts vector, so I can't use it for matrix. Currently I'm using vstack to create the first column, then second column etc. and then I use hstackto merge all columns. This takes a lot of effort, since I have to specifically add np.zeros matrices at certain places. I can't think of a better way to concatenate numpy arrays, since there are only a few functions for that and none of them really fits my problem.

pythonFriend
  • 85
  • 7
  • Could you show your existing code? It would help because it would show what result you are trying to achieve and what you don't like about the result you are getting. – roadrunner66 Apr 07 '16 at 00:03
  • There is nothing special about my code. I just have a loop which uses vstack to stack my list of np.arrays (H1,H2,H3). Then I have to manually vstack my np.zero-array with same shape as one of the H to create my first column. Then I have to manually vstack EMPTY H1 H2 H3 EMPTY for the second column, and same for the third. In the end I use a loop to hstack all of my created columns. The problem is the whole stacking part, which I can't do dynamically. My final toeplitz of toeplitz matrix constists of 25 diffents Hs and the number of EMPTY arrays will be 200+ for each column. – pythonFriend Apr 07 '16 at 00:42

3 Answers3

5

Instead of nested calls to vstack and hstack, it will be more efficient to preallocate the final array, and then use a nested loop to fill in the array. You can initially use a higher dimensional array to keep the code clean.

For example, this script

import numpy as np

H1 = np.array([[11, 11], [11, 11]])
H2 = np.array([[22, 22], [22, 22]])
H3 = np.array([[33, 33], [33, 33]])

inputs = (H1, H2, H3)

# This assumes all the arrays in `inputs` have the same shape,
# and that the data type of all the arrays is the same as H1.dtype.
nh = len(inputs)
nrows = 2*nh - 1
m, n = H1.shape
# T is a 4D array.  For a given i and j, T[i, :, j, :] is a 2D array
# with shape (m, n).  T can be intepreted as a 2D array of 2D arrays. 
T = np.zeros((nrows, m, nh, n), dtype=H1.dtype)
for i, H in enumerate(inputs):
    for j in range(nh):
        T[i + j, :, j, :] = H

# Partially flatten the 4D array to a 2D array that has the desired
# block structure.
T.shape = (nrows*m, nh*n)

print(T)

prints

[[11 11  0  0  0  0]
 [11 11  0  0  0  0]
 [22 22 11 11  0  0]
 [22 22 11 11  0  0]
 [33 33 22 22 11 11]
 [33 33 22 22 11 11]
 [ 0  0 33 33 22 22]
 [ 0  0 33 33 22 22]
 [ 0  0  0  0 33 33]
 [ 0  0  0  0 33 33]]

(Note that the result is not a Toeplitz matrix; it is a block Toeplitz matrix.)

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • That's actually what i wanted, but I struggle to understand to parameter which decides the final block toeplitz dimension. In the end it should have 256*256 dimensions. But if i change nh or nrows, I get out of bounds errors. – pythonFriend Apr 07 '16 at 10:21
  • *"...if i change nh or nrows, I get out of bounds errors."* `nh` and `nrows` are determined by the number of input arrays. As you can see in the code, `nh` is simply `len(inputs)` (i.e. the number of input arrays), and `nrows` is `2*nh - 1`. In the code that I wrote, you can't change these independently. The implementation is based on my interpretation of the question. If you need something more flexible, edit the question and explain. – Warren Weckesser Apr 07 '16 at 19:01
0

Here's an alternative approach, for anyone interested in this problem

from pylab import *
import scipy.linalg

H1 = array([[11, 11], [11, 11]])
H2 = array([[22, 22], [22, 22]])
H3 = array([[33, 33], [33, 33]])

# Setup blocks 
t = array([zeros_like(H1), H1, H2, H3])
# Create index to each block, using toeplitz
idx = scipy.linalg.toeplitz(r_[1, 2, 3, zeros(2)], r_[1, zeros(2)]).astype(int)
# Index into blocks, transpose and reshape to get re-ordered array
#  copy is used to ensure memory is nicely ordered
T = t[idx].transpose(0, 2, 1, 3).reshape(10, 6).copy()

Most of the time is spent in scipy.linalg.toeplitz, making it slower than filling in memory in an array for the small matrices used here, so I'd recommend profiling before using this approach.

user2699
  • 2,927
  • 14
  • 31
0

Here's an example using the Kronecker product that affords some flexibility for the number of rows and columns (nr, nc) desired:

import numpy as np

H1 = np.array([[11, 11], [11, 11]])
H2 = np.array([[22, 22], [22, 22]])
H3 = np.array([[33, 33], [33, 33]])

inputs = (H1, H2, H3)
nr = 5
nc = 3
m, n = H1.shape

T = np.zeros([m*nr,n*nc])
for ii, H in enumerate(inputs):
    T = T + np.kron(np.eye(nr,nc,-ii),H)

print(T)

prints

[[11. 11.  0.  0.  0.  0.]
 [11. 11.  0.  0.  0.  0.]
 [22. 22. 11. 11.  0.  0.]
 [22. 22. 11. 11.  0.  0.]
 [33. 33. 22. 22. 11. 11.]
 [33. 33. 22. 22. 11. 11.]
 [ 0.  0. 33. 33. 22. 22.]
 [ 0.  0. 33. 33. 22. 22.]
 [ 0.  0.  0.  0. 33. 33.]
 [ 0.  0.  0.  0. 33. 33.]]
gKhagb
  • 1
  • 1