Given a set of Kohn-Sham orbitals from an insulating system, we present a simple, robust, efficient and highly parallelizable method to construct a set of, optionally orthogonal, localized basis functions for the associated subspace. Our method explicitly uses the fact that density matrices associated with insulating systems decay exponentially along the off-diagonal direction in the real space representation. Our method avoids the usage of an optimization procedure, and the localized basis functions are constructed directly from a set of selected columns of the density matrix (SCDM). Consequently, the only adjustable parameter in our method is the truncation threshold of the localized basis functions. Our method can be used in any electronic structure software package with an arbitrary basis set. We demonstrate the numerical accuracy and parallel scalability of the SCDM procedure using orbitals generated by the Quantum ESPRESSO software package. We also demonstrate a procedure for combining SCDM with Hockney's algorithm to efficiently perform Hartree-Fock exchange energy calculations with near linear scaling.