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.