In recent years, fractional partial differential equation (FPDE) has been widely applied in options pricing problems, which better explains many important empirical facts of financial markets. However, the vast majority of the literatures focus on pricing the single asset option under the FPDE framework. In this paper, a two-dimensional FPDE governing the valuation of rainbow options is established, where two underlying assets are assumed to follow independent exponential Lévy processes, and its boundary conditions are determined by solving one-dimensional FPDEs. A second-order accurate finite difference scheme is proposed to discretize the two-dimensional FPDE. Given the block Toeplitz with Toeplitz block structure of the coefficient matrix, a fast preconditioned Krylov subspace method is developed for solving the resulting linear system with O(NlogN) computational complexity per iteration, where N is the matrix size. The proposed preconditioner accelerates the convergence of the iterative method with theoretical analysis. Numerical examples are given to demonstrate the accuracy and efficiency of our proposed numerical methods.