I can straight forwardly create the matrix C as (where the elements 0, A, B are 2X2 matrices)
[ 0 A 0 0 0 0 ... 0 0
B 0 A 0 0 0 ... 0 0
0 B 0 A 0 0 ... 0 0
0 0 B 0 A 0 .... 0 0
..........
0 0 0 0 0 .... B 0]
where 0 is the 2X2 matrix consisting of zeros, A, B are given 2X2 matrices (that are to be placed along the upper and lower diagonals respectively) and the matrix C is of size 2JX2J
However, it occurred to me that it is probably more efficient to directly make use of sparse matrix notation and then set the non-zero elements using something like (with J = 5)
def TheMatrix(J,alpha):
A = [[0.5, alpha],[ alpha,0.5]]
B =[[0.5,-alpha],[-alpha,0.5]]
upperDiag = A*np.ones(J-1)
lowerDiag = B*np.ones(J-1)
diagonals = [upperDiag, lowerDiag] # non-zero diagonals
ASparse = sps.diags(diagonals,offsets=[1, -1])
ASparse = sps.csc_matrix(ASparse)
However, this does not seem to work. Does anyone have any suggestions on how to proceed? In detail, the error reported (for the case J =5) is
File "myprg.py", line 16, in TheMatrix(J,alpha)
upperDiag = A*np.ones(J-1)
~^~~~~~~~~~~~~
ValueError: operands could not be broadcast together with shapes (2,2) (5,)
So, I understand that the need is to somehow replace np.ones(J-1) to something that makes it compatible with A, B being 2X2 matrices and not scalars.
4
I include the following in case it helps others: To create the sparse matrix I took up the comments from delfstack
Because a sparse matrix is a list of the positions and values of the non-zero elements, it is straight forward to accept that the sparse matrix can be formed from this information and delfstack shows how this works with the input matrix “Ain”
Ain = np.array([[16, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 5], [0, 0, 0, 0]])
delfstack suggests
Asparse = [] # initialise Asparse to the empty element
rows, cols = Ain.shape # the input matrix
for i in range(rows):
for j in range(cols):
if Ain[i][j] != 0: # identify the non-zero elements
triplet = [i, j, Ain[i][j]]
Asparse.append(triplet) # append to update Asparse
(In my example I don’t actually create the matrix A).
This seemed to mostly work, however I did hit an unexpected issue when I started to check my result by reconstructing the original matrix. This threw me until I started printing out data. When I printed out a particular triplet, it gave the result
triplet = [0, 0, np.int64(16)]
and there were problems adding np.int64(16) to 0. To overcome this I had to make the minor change
triplet = [i, j, Ain[i][j]]
triplet = list(map(int,triplet))
Asparse.append(triplet)
leading to
triplet = [0, 0, 16]
and this seemed to work.
Note added: Not sure why, but trying this today and the code now seems to work without the change.
For completeness I ran the programme using the following
sys version 3.11.9 (tags/v3.11.9:de54cf5, Apr 2 2024, 10:12:12) [MSC v.1938 64 bit (AMD64)]
python version 3.11.9
numpy version 2.0.0
scipy version 1.14.0