We consider the natural generalization of the Schr{o}dinger equation to Markovian open system dynamics: the so-called the Lindblad equation. We give a quantum algorithm for simulating the evolution of an $n$-qubit system for time $t$ within precision $epsilon$. If the Lindbladian consists of $mathrm{poly}(n)$ operators that can each be expressed as a linear combination of $mathrm{poly}(n)$ tensor products of Pauli operators then the gate cost of our algorithm is $O(t, mathrm{polylog}(t/epsilon)mathrm{poly}(n))$. We also obtain similar bounds for the cases where the Lindbladian consists of local operators, and where the Lindbladian consists of sparse operators. This is remarkable in light of evidence that we provide indicating that the above efficiency is impossible to attain by first expressing Lindblad evolution as Schr{o}dinger evolution on a larger system and tracing out the ancillary system: the cost of such a textit{reduction} incurs an efficiency overhead of $O(t^2/epsilon)$ even before the Hamiltonian evolution simulation begins. Instead, the approach of our algorithm is to use a novel variation of the linear combinations of unitaries construction that pertains to channels.