Hmm... this is a bit tricky. I'm not sure if this way will be any faster or memory efficient, but it might be worth comparing to the current approach since (I think) there's less back-and-forth between numpy and Python.

Let's assume we are pulling out a subset of rows.

1. Construct an empty `lil_matrix`, since we know the subset's final shape. `mat = lil_matrix(...)`
2. For each row we want to keep, get the starting and ending indices of its nonzero values and the columns that they're in.
3. Fill in the corresponding row in the `lil_matrix`. For example, if we were filling in the first row:

    ```python
    mat[0, h5_indices[start:end]] = h5_data[start:end]
    ```
4. Once we've processed all rows we want to keep, convert the lil_matrix to csr.

I think this approach will only work when subsetting either observations or samples, but that should be okay since we don't have a need to subset both at the same time.