Extracting the multiphase structure of the neutral interstellar medium (ISM) is key to understand the star formation in galaxies. The radiative condensation of the diffuse warm neutral medium producing a thermally unstable lukewarm medium and a dense cold medium is closely related to the initial step which leads the HI-to-H2 transition and the formation of molecular clouds. Up to now the mapping of these phases out of 21 cm emission hyper-spectral cubes has remained elusive mostly due to the velocity blending of individual cold structures present on a given line of sight. To address this problem, we developed a new Gaussian decomposition algorithm named ROHSA (Regularized Optimization for Hyper-Spectral Analysis) based on a multi-resolution process from coarse to fine grid. ROHSA uses a regularized non-linear least-square criterion to take into account simultaneously the spatial coherence of the emission and the multiphase nature of the gas. In order to obtain a solution with spatially smooth parameters, the optimization is performed on the whole data cube at once. The performances of ROHSA were tested on a synthetic observation computed from numerical simulations of thermally bi-stable turbulence. An application on a 21 cm observation of a high Galactic latitude region from the GHIGLS survey is presented. The evaluation of ROHSA on synthetic 21 cm observations shows that it is able to recover the multiphase nature of the HI and physically meaningful information about the underlying three-dimensional fields (density, velocity and temperature). The application on a real 21 cm observation of a high Galactic latitude field produces a picture of the multiphase HI, with isolated, filamentary and narrow structures and wider, diffuse and space filling components. The test-case field used here contains significant intermediate-velocity clouds that were well mapped out by the algorithm.