Many applications in numerical mathematics, scientific computing, and engineering require solving large-scale matrix problems with Toeplitz structure. Iterative solution techniques based on projection processes onto Krylov subspaces are used to solve those problems. Both the efficiency and robustness of iterative techniques can be improved by using preconditioning. The process of preconditioning is essential to most successful applications of iterative methods. In this thesis, our main goal is to develop preconditioning techniques for solving large-scale Toeplitz systems, Toeplitz matrix exponential, and apply them to option pricing problems in financial engineering. First, we study the mathematical properties of the optimal preconditioner and the generalized superoptimal preconditioner. Several existing results are extended and new properties are developed. We also consider using the normalized preconditioned conjugate gradient method with a tri-diagonal preconditioner to solve a nonsymmetric Toeplitz system, which arises from the discretization of a partial integro-differential equation (PIDE) in option pricing problems. By using the definition of family of generating functions introduced in iii [116], we prove that the tri-diagonal preconditioner leads to a superlinear convergence rate under certain conditions. Numerical results exemplify our theoretical analysis. Apart from Toeplitz systems, the Toeplitz matrix exponential (TME) also plays a key role in various application fields. In this thesis, we exploit the Krylov subspace method with the shift-invert preconditioning technique to approximate the TME. By making use of the Toeplitz structure and the famous Toeplitz matrix inversion formula, we reduce the computational cost to O(n log n) in total compared with the O(n 3 ) complexity for traditional methods. Moreover, for the nonsymmetric TME, a sufficient condition for the approximation error bound is established, which guarantees that the error bound is independent of the norm of the matrix. Numerical results are given to demonstrate the efficiency of the method. Finally, we consider approximating a matrix exponential with block Toeplitz matrix, which arises from the integration of a two dimensional PIDE in option pricing under the stochastic volatility jump diffusion model by the exponential time integration scheme. The shift-invert Arnoldi method is employed to fast approximate the matrix exponential. This results in an inner-outer iteration. To reduce the computational cost, we utilize the matrix splitting technique with multigrid method to deal with the shift-invert matrix-vector product in each inner iteration. Numerical results show that the proposed scheme is robust and efficient even compared with the existing high accurate implicit-explicit Euler based extrapolation scheme.