# Cropping in Sparse Matrix

Let $A$ be an $M \times N$ sparse matrix stored in compressed column format, in a C-like programming environment. I am interested in the best solution to get a sub-matrix of $A$. In MATLAB notation, this would read:

A_sub = A(i:i+m-1, j:j+n-1)

where $A_{sub}$ is again a sparse matrix, but this time being $m \times n$.

• I guess solutions should make use of additional data structures such as sparse accumulators to achieve linear time. – Tolga Birdal Jan 3 '14 at 10:17
• Linear time in what? Naïve approaches will achieve linear time in the number of nonzero entries. (Not being familiar with sparse accumulators, I'm interested in hearing what advantages they provide.) – Geoff Oxberry Jan 3 '14 at 19:02
• We call this MatGetSubMatrix() in PETSc. The serial algorithm is simple; the parallel implementation is more complicated. – Jed Brown Jan 3 '14 at 19:43

While you mentioned using C, Yousef Saad's package SPARSKIT has a lot of sparse matrix algorithms if you're comfortable with Fortran 77.

The compressed column format consists of three arrays, ja, ia (integers) and a (floats/doubles). For each column index j, ja(j) stores the location in the array ia where the list of all non-zero rows in column j starts.

The first step in extracting a sub-matrix A_sub from A is to count how many non-zero entries there are in each of the n columns; for that purpose, you can make an array count of all zeros and then do something to the effect of:

int k, l;
for (k=0; k<n; k++) {
for (l=A.ja[k+j]; l<A.ja[k+j+1]; l++) {
if (A.ia[l]>=i && A.ia[l]<i+m) {
count[k]++;
}
}
}

Next, you'd fill out A_sub.ja as follows:

A_sub.ja[0] = 0;
for (k=0; k<n; k++) {
A_sub.ja[k+1] = A_sub.ja[k]+count[k];
}

The next step is to fill out the arrays ia and a, which requires a second pass through A.

If you're being space-conscious, instead of making a new array count, you can store the same values starting at the address A_sub.ja[1] with no effect on the output. SPARSKIT does this if you want to see an example.

With the Eigen matrix class library (Eigen), you can efficiently extract a submatrix with a single line:

Eigen::SparseMatrix A(m,n);
Eigen::SparseMatrix B = A.block(startRow, startCol, numRows, numCols);

If A already exists in compressed sparse column (or row) format, you can extract a submatrix with code similar to this:

Eigen::MappedSparseMatrix<double> Amap(A.m, A.n, A.nnz, A.colPtrs, A.rowPtrs, A.vals);
Eigen::SparseMatrix B = Amap.block(startRow, startCol, numRows, numCols);