Seismic datasets contain valuable information that originate from areas of interest in the subsurface; such seismic reflections are however inevitably contaminated by other events created by waves reverberating in the overburden. Multi-Dimensional Deconvolution (MDD) is a powerful technique used at various stages of the seismic processing sequence to create ideal datasets deprived of such overburden effects. Whilst the underlying forward problem is well defined for a single source, a successful inversion of the MDD equations requires availability of a large number of sources alongside prior information introduced in the form of physical preconditioners (e.g., reciprocity). In this work, we reinterpret the cost function of time-domain MDD as a finite-sum functional, and solve the associated inverse problem by means of stochastic gradient descent algorithms, where gradients are computed using a small subset of randomly selected sources. Through synthetic and field data examples, the proposed method is shown to converge more stably than the conventional approach based on full gradients. Stochastic MDD represents a novel, efficient, and robust strategy to deconvolve seismic wavefields in a multi-dimensional fashion.