  In this work, we show that solvers of elliptic boundary value problems in $d$
dimensions can be approximated to accuracy $\epsilon$ from only
$\mathcal{O}\left(\log(N)\log^{d}(N / \epsilon)\right)$ matrix-vector products
with carefully chosen vectors (right-hand sides). The solver is only accessed
as a black box, and the underlying operator may be unknown and of an
arbitrarily high order. Our algorithm (1) has complexity
$\mathcal{O}\left(N\log^2(N)\log^{2d}(N / \epsilon)\right)$ and represents the
solution operator as a sparse Cholesky factorization with
$\mathcal{O}\left(N\log(N)\log^{d}(N / \epsilon)\right)$ nonzero entries, (2)
allows for embarrassingly parallel evaluation of the solution operator and the
computation of its log-determinant, (3) allows for
$\mathcal{O}\left(\log(N)\log^{d}(N / \epsilon)\right)$ complexity computation
of individual entries of the matrix representation of the solver that, in turn,
enables its recompression to an $\mathcal{O}\left(N\log^{d}(N /
\epsilon)\right)$ complexity representation. As a byproduct, our compression
scheme produces a homogenized solution operator with near-optimal approximation
accuracy. By polynomial approximation, we can also approximate the continuous
Green's function (in operator and Hilbert-Schmidt norm) to accuracy $\epsilon$
from $\mathcal{O}\left(\log^{1 + d}\left(\epsilon^{-1}\right)\right)$ solutions
of the PDE. We include rigorous proofs of these results. To the best of our
knowledge, our algorithm achieves the best known trade-off between accuracy
$\epsilon$ and the number of required matrix-vector products.
