Developing a computationally efficient but accurate building energy simulation (BES) model is important for many purposes. Model order reduction (MOR) methods are attractive and much more reliable than identification approaches, since it directly extract a lower-dimensional model from a detailed physics-based model without any pre-simulations. However, because of computational and data storage requirements, there are challenges of applying these methods to a large-scale building. To overcome the problem, this paper introduces the Krylov subspace method to the building science field. Technical issues of applying the method to building applications are addressed and a suitable algorithm that overcomes those challenges is presented. To demonstrate the reliability of the algorithm, comparisons between the resulted reduced-order model (ROM) and a high-fidelity model from a commercial BES software for a 60-zone case study building are provided. The ROM was a factor of 100 faster than the high fidelity model but with high accuracy.