  We introduce a new variational method for the numerical homogenization of
divergence form elliptic, parabolic and hyperbolic equations with arbitrary
rough ($L^\infty$) coefficients. Our method does not rely on concepts of
ergodicity or scale-separation but on compactness properties of the solution
space and a new variational approach to homogenization. The approximation space
is generated by an interpolation basis (over scattered points forming a mesh of
resolution $H$) minimizing the $L^2$ norm of the source terms; its
(pre-)computation involves minimizing $\mathcal{O}(H^{-d})$ quadratic (cell)
problems on (super-)localized sub-domains of size $\mathcal{O}(H \ln (1/ H))$.
The resulting localized linear systems remain sparse and banded. The resulting
interpolation basis functions are biharmonic for $d\leq 3$, and polyharmonic
for $d\geq 4$, for the operator $-\diiv(a\nabla \cdot)$ and can be seen as a
generalization of polyharmonic splines to differential operators with arbitrary
rough coefficients. The accuracy of the method ($\mathcal{O}(H)$ in energy norm
and independent from aspect ratios of the mesh formed by the scattered points)
is established via the introduction of a new class of higher-order Poincar\'{e}
inequalities. The method bypasses (pre-)computations on the full domain and
naturally generalizes to time dependent problems, it also provides a natural
solution to the inverse problem of recovering the solution of a divergence form
elliptic equation from a finite number of point measurements.
