The relaxation spectra, from which other material functions used to describe mechanical properties of materials can be uniquely determined, are important for modeling the rheological properties of polymers used in chemistry, food technology, medicine, cosmetics, and many other industries. The spectrum, being not directly accessible by measurement, is recovered from relaxation stress or oscillatory shear data. Only a few models and identification methods take into account the non-negativity of the real spectra. In this paper, the problem of recovery of non-negative definite relaxation spectra from discrete-time noise-corrupted measurements of relaxation modulus obtained in the stress relaxation test is considered. A new hierarchical identification scheme is developed, being applicable both for relaxation time and frequency spectra. Finite-dimensional parametric classes of models are assumed for the relaxation spectra, described by a finite series of power-exponential and square-exponential basis functions. The related models of relaxation modulus are given by compact analytical formula, described by the products of power of time and the modified Bessel functions of the second kind for the time spectrum, and by recurrence formulas based on products of power of time and complementary error functions for frequency spectrum. The basis functions are non-negative. In result, the identification task was reduced to a finite-dimensional linear-quadratic problem with non-negative unknown model parameters. To stabilize the solution, an additional smoothing constraint is introduced. Dual approach was used to solve the stated optimal identification task resulting in the hierarchical two-stage identification scheme. In the first stage, dual problem is solved in two levels and the vector of non-negative model parameters is computed to provide the best fit of the relaxation modulus to experiment data. Next, in second stage, the optimal non-negative spectrum model is determined. A complete scheme of the hierarchical computations is outlined; it can be easily implemented in available computing environments. The model smoothness is analytically studied, and the applicability ranges are numerically examined. The numerical studies have proved that using developed models and algorithm, it is possible to determine non-negative definite unimodal and bimodal relaxation spectra for a wide class of polymers. However, the examples also demonstrated that if the basis functions are non-negative and the model is properly selected for a given type of the real spectrum (unimodal, multimodal), the optimal model determined without non-negativity constraint can be non-negative in the dominant range of its arguments, especially in the wide neighborhood of the spectrum peaks.