We propose a new method to analyze and represent data recorded on a domain of gen-eral shape in R d by computing the eigenfunctions of Laplacian defined over there and expanding the data into these eigenfunctions. Instead of directly solving the eigenvalue problem on such a domain via the Helmholtz equation (which can be quite complicated and costly), we find the integral operator commuting with the Laplacian and diagonal-ize that operator. Although our eigenfunctions satisfy neither the Dirichlet nor the Neu-mann boundary condition, computing our eigenfunctions via the integral operator is simple and has a potential to utilize modern fast algorithms to accelerate the compu-tation. We also show that our method is better suited for small sample data than the Karhunen-Loève Transform/Principal Component Analysis. In fact, our eigenfunctions depend only on the shape of the domain, not the statistics of the data. As a further ap-plication, we demonstrate the use of our Laplacian eigenfunctions for solving the heat equation on a complicated domain. Key words: Laplacian eigenfunctions, boundary conditions, Green’s function, spectral decomposition, Karhunen-Loève transform, principal component analysis, heat equation