be solved as such in an automated fashion? Can discovery be computed? Although these questions may seem unorthodox their answer appears to be positive: this paper shows that this reformulation is possible and leads to a multigrid/multiresolution method/algorithm solving , up to the pre-specified level of accuracy $\epsilon$ in $H^1$-norm (i.e. finding $u^{\app}$ such that $\|u-u^{\app}\|_{H^1_0(\Omega)}\leq \epsilon \|g\|_{H^{-1}(\Omega)}$ for an arbitrary $g$ decomposed over $N$ degrees of freedom), in $\mathcal{O}\big(N \ln^{3d}\big( \max(\frac{1}{\epsilon},N^{1/d})\big)\big)$ operations (for $\epsilon\sim N^{-1/d}$, the hierarchical matrix method achieves $\epsilon$-accuracy in $L^2$ norm in $\mathcal{O}(N \ln^{2d+8} N)$ operations and the proposed multiresolution method achieves $\epsilon$-accuracy in $H^1$ norm in $\mathcal{O}(N \ln^{3d} N)$ operations). For subsequent solves (i.e. if needs to be solved for more than one $g$) then the proposed multiresolution method achieves accuracy $\epsilon \approx N^{-\frac{1}{d}}$ in $H^1$-norm in $\mathcal{O}(N \ln^{d+1} N)$ operations (we refer to Subsection \[subseccomplexity\] and in particular to Table \[tabcomplexity\] for a detailed complexity analysis of the proposed method, which can also achieve sublinear complexity if one only requires $L^2$-approximations).

The core mechanism supporting the complexity of the method presented here is the fast decomposition of $H^1_0(\Omega)$ into a direct sum of linear subspaces that are orthogonal (or near-orthogonal) with respect to the 