This paper is devoted to analysis of block multi-indexed higher-order covariance matrices, which can be used for the least-squares estimation problem. The formulation of linear and nonlinear least squares estimation problems is proposed, showing that their statements and solutions lead to generalized `normal equations', employing covariance matrices of the underlying processes. Then, we provide a class of efficient algorithms to estimate higher-order statistics (generalized multi-indexed covariance matrices), which are necessary taking in mind practical aspects of the nonlinear treatment of the least-squares estimation problem. The algorithms are examined for different higher-order and non-Gaussian processes (time-series) and an impact of signal properties on covariance matrices is analysed.