4

I have a numpy array C of size 1000 x 100. I have figured out how to compute the matrix of outer products of its rows transposed

Create array of outer products in numpy

This yields a 1000x100x100 numpy array D. I'd like to store each of the D[i,:,:] as one of the block diagonals of a sparse matrix of size 1000*100^2 x 1000*100^2. I am able to accomplish this using scipy.sparse.block_diag

E = spspar.block_diag(D,"lil")

I am wondering if there is a more efficient way of doing this. Some of the entries of E need to be manipulated at first (hence "lil"). The resulting matrix will be used to solve a sequence of sparse linear systems.

Thanks in advance.

Community
  • 1
  • 1
bigcat
  • 103
  • 1
  • 7
  • 1
    Your code takes less than five seconds to run on my machine, and seems to work fine. What are your requirements/goals for efficiency? – cge Apr 30 '15 at 23:39
  • I need to solve y = (E + Q(k))*x for x iteratively. 'k' is an iteration index, and Q(k) another very sparse matrix. The code also takes less than 5 seconds on my machine. I was hoping it would be faster. – bigcat May 01 '15 at 01:29
  • How is the time for constructing `E` compared to the solution time? `sparse` is pretty good at converting formats to suit the operation. For example the math is probably being done with `csr`. You might also time how long it takes to convert among the the different formats. The default for `block_diag` is `coo`. `coo` to `csr` is fast. to/from `lil` is slower. – hpaulj May 01 '15 at 01:36
  • `block_diag` code is Python, so you can study it. It looks like it turns each submatrix (block) into a `coo`, and then collects all their `row`,`col`,`data` values into 3 large arrays. It builds the final `coo` from that. You might gain some speed by constructing those `coo` inputs yourself, but it might not be worth the effort. – hpaulj May 01 '15 at 02:32
  • 1
    I tried streamlining `block_diag` and only got a 30-40% speedup. – hpaulj May 01 '15 at 03:29
  • @hpaul, thanks. I am wondering if the only reasonable speed up would come from a C/C++ implementation. I will implement the algorithm in Python first though to see how slow it is. – bigcat May 01 '15 at 03:42
  • The step of converting each sub matrix to sparse takes about 50% of the time, e.g. `L=[sparse.coo_matrix(d) for d in D]`. Converting `D` directly to a sparse matrix takes the same time, `sparse.coo_matrix(Dbig.reshape(1000,-1))` – hpaulj May 01 '15 at 04:04
  • 1
    I wonder if there would be a way to calculate E + Q(k) from D without actually constructing E at all, perhaps through clever indexing tricks. If so, that could be considerably faster. – cge May 01 '15 at 08:42

0 Answers0