In this article, we consider high-dimensional data which contains a low-dimensional non-Gaussian structure contaminated with Gaussian noise. Motivated by the joint diagonalization algorithms, we propose a linear dimension reduction procedure called joint low-dimensional approximation (JLA) to identify the non-Gaussian subspace. The method uses matrices whose non-zero eigen spaces coincide with the non-Gaussian subspace. We also prove its global consistency, that is the true mapping to the non-Gaussian subspace is achieved by maximizing the contrast function defined by such matrices. As examples, we will present two implementations of JLA, one with the fourth-order cumulant tensors and the other with Hessian of the characteristic functions. A numerical study demonstrates validity of our method. In particular, the second algorithm works more robustly and efficiently in most cases.