Numerical simulation of dynamical systems have been a successful means for studying complex physical phenomena. However, in large-scale settings, the system dimension makes the computations infeasible due to memory and time limitations, and ill-conditioning. The remedy is model reduction. This dissertation focuses on projection methods to efficiently construct reduced order models for large linear dynamical systems. A modified cyclic low-rank Smith method is introduced to compute low-rank approximations to solutions of large-scale Lyapunov equations. Unlike the original cyclic low-rank Smith method of Penzl, the number of columns in the modified approximant does not necessarily increase at each step and is much lower. Fundamental convergence results are established for the errors in the approximate solutions and also in the approximate Hankel singular values. For positive real balancing, this work derives a multiplicative error bound and develops a modified scheme with an absolute error bound for a certain subclass of positive real systems. Moreover, a frequency weighted balancing method with guaranteed stability and a simple Hinfinity error bound is introduced. Unlike the existing approaches, the method avoids the explicit computation of the input and output weightings. This dissertation derives an exact expression for the H2 norm of the error system of the Lanczos procedure, the first such result for Krylov based methods. The resulting expression shows that the H2 error is due to the mismatch at the mirror images of the poles of the original and reduced systems, and hence suggests choosing the mirror images as the interpolation points for the rational Krylov method. In addition two algorithms are proposed to overcome the rank deficiencies occurring in the MIMO version of the rational Krylov method. Finally, a novel model reduction algorithm by least-squares is developed, one of the cornerstones of this dissertation. The method is a projection and combines Krylov and singular value decomposition methods. The reduced model is asymptotically stable, matches a certain number of moments; and minimizes a weighted H2 error in the discrete time case. The effectiveness of the proposed approaches is tested by means of various numerical experiments.