Creating a sparse matrix from lists of sub matrices (Python)

Question

This is my first SO question ever. Let me know if I could have asked it better :)

I am trying to find a way to splice together lists of sparse matrices into a larger block matrix.

I have python code that generates lists of square sparse matrices, matrix by matrix. In pseudocode:

Lx = [Lx1, Lx1, ... Lxn]
Ly = [Ly1, Ly2, ... Lyn]
Lz = [Lz1, Lz2, ... Lzn]

Since each individual Lx1, Lx2 etc. matrix is computed sequentially, they are appended to a list--I could not find a way to populate an array-like object "on the fly".

I am optimizing for speed, and the bottleneck features a computation of Cartesian products item-by-item, similar to the pseudocode:

M += J[i,j] * [ Lxi *Lxj + Lyi*Lyj + Lzi*Lzj ]

for all combinations of 0 <= i, j <= n. (J is an n-dimensional square matrix of numbers).

It seems that vectorizing this by computing all the Cartesian products in one step via (pseudocode):

L = [ [Lx1, Lx2, ...Lxn],
      [Ly1, Ly2, ...Lyn],
      [Lz1, Lz2, ...Lzn] ]
product = L.T * L

would be faster. However, options such as np.bmat, np.vstack, np.hstack seem to require arrays as inputs, and I have lists instead.

Is there a way to efficiently splice the three lists of matrices together into a block? Or, is there a way to generate an array of sparse matrices one element at a time and then np.vstack them together?

Reference: Similar MATLAB code, used to compute the Hamiltonian matrix for n-spin NMR simulation, can be found here:

http://spindynamics.org/Spin-Dynamics---Part-II---Lecture-06.php

Where are the little matrices coming from, and in what sparse matrix format? — user2357112, Apr 08 '16 at 15:56
Each Lxn matrix is computed by repetitive krons of these base matrices: `code` sigma_x = csc_matrix(np.matrix([[0, 1/2], [1/2, 0]])) sigma_y = csc_matrix(np.matrix([[0, -1j/2], [1j/2, 0]])) sigma_z = csc_matrix(np.matrix([[1/2, 0], [0, -1/2]])) unit = csc_matrix(np.matrix([[1, 0], [0, 1]])) `code` — Geoffrey Sametz, Apr 08 '16 at 16:03
Each Lxn matrix is computed by repetitive krons of these base matrices: `sigma_x = csc_matrix(np.matrix([[0, 1/2], [1/2, 0]]))`, `sigma_y = csc_matrix(np.matrix([[0, -1j/2], [1j/2, 0]]))`, `sigma_z = csc_matrix(np.matrix([[1/2, 0], [0, -1/2]]))`, `unit = csc_matrix(np.matrix([[1, 0], [0, 1]]))` — Geoffrey Sametz, Apr 08 '16 at 16:12

user2357112 · Answer 1 · 2016-04-09T02:00:44.843

1

This is scipy.sparse.bmat:

L = scipy.sparse.bmat([Lx, Ly, Lz], format='csc')
LT = scipy.sparse.bmat(zip(Lx, Ly, Lz), format='csr') # Not equivalent to L.T
product = LT * L

edited Apr 09 '16 at 02:00

answered Apr 08 '16 at 16:09

user2357112

260,549
28
431
505

That blends everything into one matrix. I actually need to keep the boundaries between individual sub matrices Lx[n]/Ly[n]/Lz[n] . For example, L[1,0] will contain [Lx1, Ly1, Lz1] x [Lx0, Ly0, Lz0]. Then, the number in J[1,0] will be multiplied with this matrix stored in L[1,0]. I suppose you could multiply the J values with slices of the mega-matrix (thus re-creating the submatrices), but that seems clunky. – Geoffrey Sametz Apr 08 '16 at 17:18
@GeoffreySametz: You can't "keep the boundaries" like you're thinking of, especially with sparse matrices, which are inherently two-dimensional. Even with dense arrays, "keeping the boundaries" would make things really clunky to work with. – user2357112 Apr 08 '16 at 17:27
Thanks. I've started implementing a slice approach and will see how that goes. – Geoffrey Sametz Apr 08 '16 at 17:44
It turns out this won't work, because the transpose of the sum matrix isn't the same as the transpose of the block elements: http://math.stackexchange.com/questions/246289/transpose-of-block-matrix – Geoffrey Sametz Apr 09 '16 at 00:46
1

@GeoffreySametz: Oh, indeed. That's not really `L.T` you want there. I interpreted your pseudocode differently from how you did. The new version of the answer should produce the correct result. – user2357112 Apr 09 '16 at 02:02
That is much more elegant than my solution, although I had to covert the tuples to lists for bmat to understand them: `L_row = bmat([list(l) for l in zip(Lx, Ly, Lz)], format='csr')` Unfortunately, the bottleneck of the final dot product remains slower than the original unvectorized code. – Geoffrey Sametz Apr 09 '16 at 03:01

score 0 · Answer 2 · answered Apr 09 '16 at 02:16

I have a "vectorized" solution, but it's almost twice as slow as the original code. Both the bottleneck shown above, and the final dot product shown in the last line below, take about 95% of the calculation time according to kernprof tests.

    # Create the matrix of column vectors from these lists
L_column = bmat([Lx, Ly, Lz], format='csc')
# Create the matrix of row vectors (via a transpose of matrix with
# transposed blocks)
Lx_trans = [x.T for x in Lx]
Ly_trans = [y.T for y in Ly]
Lz_trans = [z.T for z in Lz]
L_row = bmat([Lx_trans, Ly_trans, Lz_trans], format='csr').T
product = L_row * L_column

Geoffrey Sametz · Accepted Answer · 2016-04-10T02:38:18.290

I was able to get a tenfold speed increase by not using sparse matrices and using an array of arrays.

Lx = np.empty((1, nspins), dtype='object') Ly = np.empty((1, nspins), dtype='object') Lz = np.empty((1, nspins), dtype='object')

These are populated with the individual Lx arrays (formerly sparse matrices) as they are generated. Using the array structure allows the transpose and Cartesian product to perform as desired:

Lcol = np.vstack((Lx, Ly, Lz)).real Lrow = Lcol.T # As opposed to sparse version of code, this works! Lproduct = np.dot(Lrow, Lcol)

The individual Lx[n] matrices are still "bundled", so Product is an n x n matrix. This means in-place multiplication of the n x n J array with Lproduct works:

scalars = np.multiply(J, Lproduct)

Each matrix element is then added on to the final hamiltonian matrix:

for n in range(nspins): for m in range(nspins): M += scalars[n, k].real

Creating a sparse matrix from lists of sub matrices (Python)

3 Answers3