Thermal rate constants k(T) and cumulative reaction probabilities N(E) can be computed as a sum of correlation functions C-nm = <(n)/f((H) over tilde/phi (m)). In this paper we discuss the use of two different Krylov subspace methods to compute these correlation functions for large systems. The first approach is based on the Lanczos algorithm to transform the Hamiltonian to tridiagonal form. As shown by Mandelshtam (J. Chem, Phys. 1998, 108, 9999) and Chen and Guo (J, Chem. Phys. 1999, 111, 9944), ail correlation functions can be computed from a single recursion. The second approach treats a number of linear systems of equations using a Krylov subspace solver. Here the quasiminimal residual (QMR) method was used. For the first approach, we found that we needed the same number of Lanczos recursions as the size of the matrix. If Ilo re-orthogonalization is used, the number of recursions grows further. The linear solver approach, on the other hand, converges fast for each linear system, but many systems must be solved.