Hierarchical cubature is a new method for achieving linear scaling computation of the exchange-correlation matrix central to Density Functional Theory. Hierarchical cubature combines a k-dimensional generalization of the binary search tree with adaptive numerical integration involving an entirely Cartesian grid. Hierarchical cubature overcomes strong variations in the electron density associated with nuclear cusps through multiresolution rather than spherical-polar coordinate transformations. This unique Cartesian representation allows use of the exact integration error during grid construction, supporting O(log N) range-queries that exploit locality of the Cartesian Gaussian based electron density. Convergence is controlled by tau (r), which bounds the local integration error of the electron density. An early onset of linear scaling is observed for RB3LYP/6-31G** calculations on water clusters, commencing at (H2O)(30) and persisting with decreasing values of tau (r). Comparison with nuclear weight schemes suggests that the new method is competitive on the basis of grid points per atom. Systematic convergence of the RPBE0/6-31G** Ne-2 binding curve is demonstrated with respect to tau (r).