In this article we consider the smoothing problem for hidden Markov models. Given a hidden Markov chain {X-n}(n >= 0) and observations {Y-n}(n >= 0), our objective is to compute E[phi(X-0, ..., X-k)vertical bar y(0), ..., y(n)] for some real-valued, integrable functional phi and k fixed, k << n and for some realization (y(0), ..., y(n)) of (Y-0, ..., Y-n). We introduce a novel application of the multilevel Monte Carlo method with a coupling based on the Knothe-Rosenblatt rearrangement. We prove that this method can approximate the aforementioned quantity with a mean square error (MSE) of O(epsilon(2)) for arbitrary epsilon > 0 with a cost of O(epsilon(-2)). This is in contrast to the same direct Monte Carlo method, which requires a cost of O(n epsilon(-2)) for the same MSE. The approach we suggest is, in general, not possible to implement, so the optimal transport methodology of [A. Spantini, D. Bigoni, and Y. Marzouk, T. Mach. Learn. Res., 19 (2018), pp. 2639-2709; M. Parno, T. Moselhy, and Y. Marzouk, SIAM/ASA J. Uncertain. Quantif., 4 (2016), pp. 1160-1190] is used, which directly approximates our strategy. We show that our theoretical improvements are achieved, even under approximation, in several numerical examples.