title: Parallel Construction of Compressed Sparse Column Matrix
body: <p>Compressed Sparse Column Matrix (csc) with m columns, as defined <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html" rel="nofollow noreferrer">here</a>.</p>
  
  <p>It has three arrays: Data, Indices, and Indptr.</p>
  
  <ul>
  <li>indptr is a 1D array with a size of m+1</li>
  <li>indices is row indices of all non-zero entries of the matrix.</li>
  <li>indices[indptr[i]:indptr[i + 1]] is row indices of all non-zero entries of the i th column. The values of these non-zero entries are defined in data[indptr[i]:indptr[i + 1]]</li>
  </ul>
  
  <p>Compressed sparse row (csr) matrix with n rows has a similar <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html" rel="nofollow noreferrer">definition</a>.</p>
  
  <p>Problem:</p>
  
  <p>Cpu works on n jobs. The p th job produces data for the q th row of a sparse matrix (csc or csr) with n rows and m columns. p and q don't have to be equal, as long as which job corresponds to which row is recorded.</p>
  
  <p>How to efficiently collect the data from the jobs and produce a csc matrix?</p>
  
  <p>Current solutions:</p>
  
  <p>Approach A.</p>
  
  <ol>
  <li>Empty arrays called shared_data, shared_indices, shared_indptr (they defines a csr matrix), shared_ordering are created. Sizes of shared_data and shared_indices depends on an estimation of number of nonzero entries. Set shared_indptr[0] = 0. These arrays are shared by jobs. </li>
  <li>Counters of number of non-zero entries (count_nz) and number of rows entered to the csr matrix (count_row) are also shared by jobs.</li>
  <li>Accessing resources described in 1 and 2 requires acquiring a lock.</li>
  <li>Each job knows the column indices of the nonzero entries (job_indices) that it tries to put into the csr matrix. It also knows the values of the nonzero entries (job_data), and the number of the nonzero entries (job_nz).</li>
  <li>Each job tries to acquire lock.</li>
  <li>Upon acquiring the lock, a job does the following. Then release lock </li>
  </ol>
  
  <p>.</p>
  
  <pre><code>after acquiring lock
  shared_data[counter_nz:(counter_nz + job_nz)] = job_data
  shared_indices[counter_nz:(counter_nz + job_nz)] = job_indices
  shared_indptr[count_row + 1] = shared_indptr[count_row] + job_nz
  shared_ordering[counter_row] = job_id
  counter_nz  += job_nz
  counter_row += 1
  release lock.
  </code></pre>
  
  <ol start="7">
  <li><p>create a csr matrix object by wrapping the shared arrays (like scipy.sparse.csr_matrix((shared_data, shared_indices, shared_indptr)).</p></li>
  <li><p>Convert the csr to a csc matrix (the goal is to make a csc matrix  instead of a csr matrix ...)</p></li>
  </ol>
  
  <p>The problem with this approach is the lock really kills performance. It is like a dozen people tries to grab the same tool to work on different parts of a project.</p>
  
  <p>Approach B</p>
  
  <p>Each cpu works on m jobs and build a csc matrix with the method described in A. This way uses much more memory, but the matrices are generated more quickly.</p>
  
  <p>Approach C</p>
  
  <p>One cpu works on m job and build a csc matrix. This takes roughly twice more time than approach A on a 12 cpu machine.</p>
  
  <p>Approach D (Update 06/05/2016)</p>
  
  <p>There are s cpu cores. The shared_data, shared_indptr, and shared_indices arrays of the final csc matrix and the shared_ordering array are shared by all cpu cores. No lock needed. Each cpu core runs a single subprocess (instead of starting and stopping subprocess all the time) Each cpu core process 1/s of the m jobs and each cpu construct a small csc matrix with the method described in A. After a cpu finished all jobs, it broadcast the three sizes of the data that it put into each of the three arrays of its small csc matrix to other cpus. Once a cpu receives this information from all other cpu (there are s - 1 info to receive from other cpus), it calculates the position of its three arrays of its small csc matrix in the three arrays of the final csc matrix. Then it copies its small csc matrix to the shared csc matrix. </p>
  
  <p>After that, the cpu sends a signal to other cpu. Once a cpu receives s - 1 signals, it starts the next cycle. (probably don't need this 2nd synchronization.)</p>
  
  <p>Hopefully the synchronizing would waste much less time than the lock in A. And that this approach scales linearly from 1 cpu core to 12 cpu core. If anyone has experience with this, feel free to give suggestions.</p>
  
  <p>A draw back of approach D is that it uses twice the amount of memory than the memory require to hold the final csc matrix. Half of the memory is for the shared arrays. The other half is used by all the arrays of the small csc matrices... This seems ok. In approach B, each cpu only has about 1/s of whole memory to use. Approach D is better.</p>
  
  <p>Can keep the arrays of the small matrics as buffer so no time waste for recreating numpy arrays.</p>
  
  <p>I am trying to work with B and optimize each job's memory usage per cpu core. But if approach A can avoid killing performance by a smarter design of lock, I would use A... I guess I will go with D.
  Related question:
  <a href="https://stackoverflow.com/questions/22796118/efficiently-populate-scipy-sparse-matrix-from-subset-of-dictionary?rq=1">Efficiently populate SciPy sparse matrix from subset of dictionary</a></p>