  We describe a randomized algorithm for producing a near-optimal hierarchical
off-diagonal low-rank (HODLR) approximation to an $n\times n$ matrix
$\mathbf{A}$, accessible only though matrix-vector products with $\mathbf{A}$
and $\mathbf{A}^{\mathsf{T}}$. We prove that, for the rank-$k$ HODLR
approximation problem, our method achieves a $(1+\beta)^{\log(n)}$-optimal
approximation in expected Frobenius norm using $O(k\log(n)/\beta^3)$
matrix-vector products. In particular, the algorithm obtains a
$(1+\varepsilon)$-optimal approximation with $O(k\log^4(n)/\varepsilon^3)$
matrix-vector products, and for any constant $c$, an $n^c$-optimal
approximation with $O(k \log(n))$ matrix-vector products. Apart from
matrix-vector products, the additional computational cost of our method is just
$O(n \operatorname{poly}(\log(n), k, \beta))$. We complement the upper bound
with a lower bound, which shows that any matrix-vector query algorithm requires
at least $\Omega(k\log(n) + k/\varepsilon)$ queries to obtain a
$(1+\varepsilon)$-optimal approximation.
  Our algorithm can be viewed as a robust version of widely used "peeling"
methods for recovering HODLR matrices and is, to the best of our knowledge, the
first matrix-vector query algorithm to enjoy theoretical worst-case guarantees
for approximation by any hierarchical matrix class. To control the propagation
of error between levels of hierarchical approximation, we introduce a new
perturbation bound for low-rank approximation, which shows that the widely used
Generalized Nystr\"om method enjoys inherent stability when implemented with
noisy matrix-vector products. We also introduced a novel randomly perforated
matrix sketching method to further control the error in the peeling algorithm.
