Matrix Storage and Design

16.1. Matrix Storage and Design#

16.1.1. Sparse, Structured, and Dense Matrices#

In numerical computation, we often encounter matrices that have special properties. We can broadly classify them into three types:

  • A sparse matrix is a matrix with enough zero-valued entries that it is computationally efficient to store and manipulate the matrix by taking advantage of these zeros. [Wilkinson]

\[\begin{split} \begin{bmatrix} \bullet & & \bullet & \bullet & \bullet \\ & \bullet & & & \\ \bullet & & \bullet & & \\ \bullet & & & \bullet & \\ \bullet & & & & \bullet \end{bmatrix} \end{split}\]

  • A structured matrix has entries that follow a predictable mathematical pattern, which can be exploited for storage and computation (e.g., Toeplitz, Circulant, or Diagonal matrices).

\[\begin{split} \begin{bmatrix} a & b & c & d & e \\ b & a & b & c & d \\ c & b & a & b & c \\ d & c & b & a & b \\ e & d & c & b & a \end{bmatrix} \end{split}\]

  • A dense matrix is neither sparse nor structured, meaning most of its entries are non-zero and do not follow a simple pattern. These are stored as standard 2D arrays.

\[\begin{split} \begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \\ \end{bmatrix} \end{split}\]

16.1.2. Design Goals for Sparse Matrices#

When designing sparse matrix formats and libraries, we have several (sometimes conflicting) goals:

  • Usability: Most operations (*, \, +, etc.) should give the same results for sparse and dense matrices. The type should propagate automatically (e.g., sparse + dense = dense).

  • Simplicity & Robustness: The interface should be easy to use and handle edge cases (like all-zero matrices) gracefully.

  • Storage Cost: The memory required should be proportional to the number of non-zeros, \(\text{nnz}\). We write this as \(O(\text{nnz})\).

  • Performance: The time for sparse operations (like matrix-vector multiply) should be proportional to the number of non-zeros, \(O(\text{nnz})\), not the dense size, \(O(n^2)\).

16.1.3. Common Data Structures for Sparse Matrices#

When we move from a dense Array to a sparse format, we must decide how to store only the non-zero values. As a running example, let’s consider how to store the following 5x5 matrix, which has 10 non-zero entries (nnz=10).

\[\begin{split} \begin{bmatrix} 5 & 0 & -3 & -2 & 7 \\ 0 & 5 & 0 & 0 & 0 \\ -2 & 0 & -1 & 0 & 0 \\ -4 & 0 & 0 & -10 & 0 \\ 0 & 0 & 0 & 0 & 9 \end{bmatrix} \end{split}\]

16.1.3.1. 1. Full (Dense) Storage#

  • Used in the standard Array type.

  • Stores all matrix entries, including the zeros.

  • Memory Usage: \(O(m \cdot n)\) for an \(m \times n\) matrix.

16.1.3.2. 2. List of Lists (LIL)#

  • Store only the non-zero values, organized by row or by column.

  • This example shows a column-oriented LIL:

  Column 1 -> Rows [1,3,4], Values [5,-2,-4]
  Column 2 -> Rows [2], Values [5]
  Column 3 -> Rows [1,3], Values [-3,-1]
  Column 4 -> Rows [1,4], Values [-2,-10]
  Column 5 -> Rows [1,5], Values [7,9]

  • Pro: Very easy to insert new elements (good for incremental matrix construction).

  • Pro: This is essentially the adjacency list format for the graph represented by the matrix.

  • Con: Can be inefficient for matrix-vector products compared to other formats, as it involves many small, separate lists.

  • Memory Usage: \(O(\text{nnz} + n)\) for a column-oriented format with \(n\) columns (one list per column, plus storage for the non-zeros).

16.1.3.3. 3. Coordinate List (COO)#

  • Store three lists of length \(\text{nnz}\): one for row indices (I), one for column indices (J), and one for the values (V).

  (Row, Col, Val)
  [1, 1, 5]
  [3, 1, -2]
  [4, 1, -4]
  [2, 2, 5]
  [1, 3, -3]
  [3, 3, -1]
  [1, 4, -2]
  [4, 4, -10]
  [1, 5, 7]
  [5, 5, 9]

  • Pro: Also very simple and easy to construct incrementally (just append to the three lists).

  • Con: Must be sorted or converted to another format (like CSC) for efficient arithmetic or look-ups. Operations on unsorted COO are very slow.

  • Memory Usage: \(O(\text{nnz})\) (specifically, \(3 \times \text{nnz}\) storage).

16.1.3.4. 4. Compressed Sparse Column (CSC)#

This is the most common format for high-performance computing. It is essentially the COO format, but with the column indices compressed.

Instead of storing the column index for every non-zero, we store three arrays:

  • nzval: A list of all non-zero values, ordered column by column.

  • rowval: A list of the row indices for each corresponding value in nzval.

  • colptr: An array of \(n+1\) pointers (indices) that tell us where each column starts in the nzval and rowval arrays.

  nzval  = [5,-2,-4, 5, -3,-1, -2,-10, 7, 9]   # Non-zero values
  rowval = [1, 3, 4, 2,  1, 3,  1,  4, 1, 5]   # Row indices for each value
  colptr = [1, 4, 5, 7, 9, 11]                # Pointers to the start of each column

  • The colptr array contains \(n+1\) entries. The non-zeros for column k are stored in nzval and rowval from index colptr[k] up to colptr[k+1]-1.

  • This design means colptr[1] is always 1, and colptr[n+1] is always nnz+1. The number of elements in column k is simply colptr[k+1] - colptr[k].

  • Pro: Highly efficient for arithmetic operations, column slicing, and matrix-vector products.

  • Con: Very inefficient for incremental construction (inserting a new non-zero requires rebuilding large parts of the arrays).

  • Con: Element look-up \(A[i,j]\) is relatively slow, \(O(\log(\text{nnz}_j))\) (where \(\text{nnz}_j\) is the number of non-zeros in column \(j\)), as you must search the column.

  • Memory Usage: \(O(\text{nnz} + n)\).

  • Note: This is the format used by Julia’s SparseMatrixCSC type and many other libraries (e.g., in Python’s scipy.sparse). A similar Compressed Sparse Row (CSR) format is also common.

Contents