Fastest numerical method to solve Lindblad Master Equation?

The Lindblad Master Equation is a generalization of the Schrodinger Equation for open quantum systems, given by

$$
frac{mathrm{d} rho}{mathrm{d}t} = -i left[ H, rhoright] + sum_k gamma_k left( L_k rho L_k^dagger – frac{1}{2} left{ L_k^dagger L_k, rho right}right) = mathcal{L}(rho)
$$
where $rho$ is the density matrix describing the state of the system, $H$ is some Hamiltonian, and $L_i$ are the jump operators. All of these objects are of size $N times N$ where $N$ is the Hilbert size of the system.

Question: Assuming that everything is time-independent (constant Hamiltonian and jump operators), what would be the fastest numerical method to solve this equation? More specifically, starting from a known initial density matrix $rho_0$, how can I get the density matrix at time $t_f$, $rho(t_f)$?

Starting point: I already know of two possible methods.

(1) Computing the complete Lindbladian propagator, i.e. $Lambda(t_f) = exp(t_f mathcal{L})$ and then doing the simple superoperator-matrix product $rho_f = Lambda(t_f) rho_0$. However, $mathcal{L}$ is a superoperator of size $N times N times Ntimes N$, so computing its exponential does not seem numerically very fast.

(2) Making some Kraus operator propagation, i.e.

$$
rho(t+mathrm{d}t) = sum_nu M_nu rho(t) M_nu^dagger
$$
where $M_nu$ are the Kraus operators, given by
$$
M_0 = I – mathrm{d}t left(i H + frac{1}{2} sum_k L_k^dagger L_kright) quad text{and} quad M_nu = sqrt{mathrm{d}t} L_nu
$$
which is a trace preserving method up to second order in $mathrm{d}t$. The good thing about this is that it only requires matrix-matrix products. However, we do need to take sufficiently many time steps, and thus perform quite a lot of matrix-matrix products.

Is there any known fast method for such problems?