We analyze statistical features of the 'optimization landscape' in a random version of one of the simplest constrained optimization problems of the least-square type: finding the best approximation for the solution of a system of M linear equations in N unknowns: (a(k), x)= b(k), k = 1, . . . , M on the N-sphere x(2) = N. We treat both the N-component vectors ak and parameters bk as independent mean zero real Gaussian random variables. First, we derive the exact expressions for the mean number of stationary points of the least-square loss function in the overcomplete case M > N in the framework of the Kac-Rice approach combined with the random matrix theory for Wishart ensemble. Then we perform its asymptotic analysis as N -> infinity at a fixed alpha = M/N > 1 in various regimes. In particular, this analysis allows to extract the large deviation function for the density of the smallest Lagrange multiplier lambda(min) associated with the problem, and in this way to find its most probable value. This can be further used to predict the asymptotic mean minimal value epsilon(min) of the loss function as N -> infinity. Finally, we develop an alternative approach based on the replica trick to conjecture the form of the large deviation function for the density of epsilon(min) at N >> 1 and any fixed ratio alpha = M/N > 0. As a by-product, we find the compatibility threshold alpha(c) < 1 which is the value of a beyond which a large random linear system on the N-sphere becomes typically incompatible.