An accurate characterization of tissue residue function R(t) in bolus-tracking magnetic resonance imaging is of crucial importance to quantify cerebral hemodynamics. R(t) estimation requires to solve a deconvolution problem. The most popular deconvolution method is singular value decomposition (SVD). However, SVD is known to bear some limitations, e.g., R(t) profiles exhibit nonphysiological oscillations and take on negative values. In addition, SVD estimates are biased in presence of bolus delay and dispersion. Recently, other deconvolution methods have been proposed, in particular block-circulant SVD (cSVD) and Tikhonov regularization (TIKH). Here we propose a new method based on nonlinear stochastic regularization (NSR). NSR is tested on simulated data and compared with SVD, cSVD, and TIKH in presence and absence of bolus dispersion. A clinical case in one patient has also been considered. NSR is shown to perform better than SVD, cSVD, and TIKH in reconstructing both the peak and the residue function, in particular when bolus dispersion is considered. In addition, differently from SVD, cSVD, and TIKH, NSR always provides positive and smooth R(t).