  There has been a rise in the popularity of algebraic methods for graph
algorithms given the development of the GraphBLAS library and other sparse
matrix methods. An exemplar for these approaches is Breadth-First Search (BFS).
The algebraic BFS algorithm is simply a recurrence of matrix-vector
multiplications with the $n \times n$ adjacency matrix, but the many redundant
operations over nonzeros ultimately lead to suboptimal performance. Therefore
an optimal algebraic BFS should be of keen interest especially if it is easily
integrated with existing matrix methods.
  Current methods, notably in the GraphBLAS, use a Sparse Matrix masked-Sparse
Vector (SpMmSpV) multiplication in which the input vector is kept in a sparse
representation in each step of the BFS, and nonzeros in the vector are masked
in subsequent steps. This has been an area of recent research in GraphBLAS and
other libraries. While in theory these masking methods are asymptotically
optimal on sparse graphs, many add work that leads to suboptimal runtime. We
give a new optimal, algebraic BFS for sparse graphs, thus closing a gap in the
literature.
  Our method multiplies progressively smaller submatrices of the adjacency
matrix at each step. Let $n$ and $m$ refer to the number of vertices and edges,
respectively. On a sparse graph, our method takes $O(n)$ algebraic operations
as opposed to $O(m)$ operations needed by theoretically optimal sparse matrix
approaches. Thus for sparse graphs it matches the bounds of the best-known
sequential algorithm and on a Parallel Random Access Machine (PRAM) it is
work-optimal. Our result holds for both directed and undirected graphs.
Compared to a leading GraphBLAS library our method achieves up to 24x faster
sequential time and for parallel computation it can be 17x faster on large
graphs and 12x faster on large-diameter graphs.
