In this article we consider the smoothing problem for hidden Markov models (HMM).Given a hidden Markov chain {Xn} n≥0 and observations {Yn} n≥0 , our objective is to compute E[ϕ(X0, . . . , X k )|y0, . . . , yn] for some real-valued, integrable functional ϕ and k fixed, k n and for some realisation (y0, . . . , yn) of (Y0, . . . , Yn). We introduce a novel application of the multilevel Monte Carlo (MLMC) method with a coupling based on the Knothe-Rosenblatt rearrangement. We prove that this method can approximate the afore-mentioned quantity with a mean square error (MSE) of O( 2 ), for arbitrary > 0 with a cost of O( −2 ). This is in contrast to the same direct Monte Carlo method, which requires a cost of O(n −2 ) for the same MSE. The approach we suggest is, in general, not possible to implement, so the optimal transport methodology of[17] is used, which directly approximates our strategy. We show that our theoretical improvements are achieved, even under approximation, in several numerical examples.