(1) wavelet based methods [@Beylkin:1995; @Beylkin:1998; @Beylkin:1998b; @DorobantuEngquist1988; @EngquistRunborg02] (2) the Fast Multipole Method [@GreengardRokhlin:1987] and (3) Hierarchical matrices [@HackbuschGrasedyck:2002; @Bebendorf:2008]. Although methods based on (classical) wavelets achieve a multiresolution compression of the solution space of in $L^2$ and although approximate wavelets and approximate $L^2$ projections can stabilize hierarchical basis methods [@VassilevskiWang1997a; @VassilevskiWang1998], their applications to are limited by the facts that (a) the underlying wavelets can perform arbitrarily badly [@BaOs:2000] in their $H^1_0(\Omega)$ approximation of the solution space and (b) the operator does not preserve the orthogonality between subscales/subbands with classical wavelets. The Fast Multipole Method and hierarchical matrices exploit the property that sub-matrices of the inverse discrete operator are low rank away from the diagonal. This low rank property can be rigorously proven for (based on the approximation of its Green’s function by sums of products of harmonic functions [@Bebendorf2005]) and leads to provable convergence (with rough coefficients), up to the pre-specified level of accuracy $\epsilon$ in $L^2$-norm, in $\mathcal{O}(N \ln^{6} N \ln^{2d+2} \frac{1}{\epsilon})$ operations ([@Bebendorf2005] and [@Bebendorf:2008 Thm. 2.33 and Thm. 4.28]). Can the problem of finding a fast solver for be, to some degree, reformulated as an Uncertainty Quantification/Decision Theory problem that could, to some degree, 