We introduce a systematic approach to construct localized finite element basis functions for second and higher order elliptic problems with heterogeneous multiscale coefficients. The basis functions are energy minimizing functions on local patches. On a regular mesh with mesh size h, the basis functions have supports of diameter h(log(1/h)) and have the optimal convergence rate O(h^k) in the energy norm for elliptic operators of order 2k. From the perspective of operator compression, these localized basis functions provide an accurate approximation of the principal eigen-space of the corresponding elliptic operators, and thus provide an efficient and optimal way to compress any order elliptic operators with localized basis. From the perspective of the sparse PCA, our results show that a large set of covariance functions can be approximated by a rank-n operator with optimal accuracy, i.e., on the order of its n-th eigenvalue arranged in the descending order. The range space of the rank-n operator can be spanned by this localized basis. In contrast to the traditional error analysis of Finite Element Methods that utilizes the interpolation type polynomial approximation property in the Sobolev space, we use the projection type polynomial approximation property in the Sobolev space in establishing the error analysis for our new method. Our method provides a solid mathematical foundation to build new multigrid methods for solving linear systems and eigenvalue decomposition for elliptic operators. This is a joint work with Pengchuan Zhang.