We introduce a new variational method for the numerical homogenization of divergence formelliptic, parabolic and hyperbolic equations with arbitrary rough (L∞)coefficients. Our method does not rely on concepts of ergodicity or scale-separation buton compactness properties of the solution space and a new variational approach tohomogenization. The approximation space is generated by an interpolation basis (overscattered points forming a mesh of resolution H) minimizing the L2 norm of thesource terms; its (pre-)computation involves minimizing 𝒪(H−d)quadratic (cell) problems on (super-)localized sub-domains of size 𝒪(H ln(1/H)).The resulting localized linear systems remain sparse and banded. The resultinginterpolation basis functions are biharmonic for d ≤ 3, and polyharmonic ford ≥ 4, forthe operator −div(a∇·) and can be seen as a generalization ofpolyharmonic splines to differential operators with arbitrary rough coefficients. Theaccuracy of the method (𝒪(H)in energy norm and independent from aspect ratios of the mesh formed by the scatteredpoints) is established via the introduction of a new class ofhigher-order Poincaré inequalities. The method bypasses (pre-)computations on the fulldomain and naturally generalizes to time dependent problems, it also provides a naturalsolution to the inverse problem of recovering the solution of a divergence form ellipticequation from a finite number of point measurements.