First of all, I don't know anything about Stein-Haff estimator other than what I saw from a few seconds of Googling in https://stat.duke.edu/~berger/papers/yang.pdf , which contains the quote "This estimator has two problems. First, the intuitively compatible ordering $\phi_1 \geq \phi_2 \geq \dots \geq \phi_p$ is frequently violated. Second, and more serious, some of the $\phi_i$ may even be negative. Stein suggests an isotonizing algorithm to avoid these problems. ... The details of this isotonizing algorithm can be found in Lin and Perlman (1985)".
That reference is: LIN, S. P. and PERLMAN,M. D. (1985). A Monte Carlo comparison of four estimators for a covariance matrix. In Multivariate Analysis 6 (P. R. Krishnaiah, ed.) 411-429. North-Holland, Amsterdam.
However, I do know about optimization. Isotonizing constraints can be placed on a least squares problem, making it into a (linearly constrained) convex Quadratic Programming (QP) problem, which is easy to formulate and numerically solve using off the shelf software. If an L^1 norm is used for the regression, or even an L^1 penalty being added to an L^2 objective, that is still a convex QP. In the case in which the objective is solely L^1, it would actually be a Linear Programming (LP) problem, which is a special case of a convex QP.
As for the negative eigenvalues, presuming those are still possible after adding the isotonizing constraints, that can be dealt with by imposing a semidefinite constraint on the covariance matrix. I.e., imposing the constraint that the minimum eigenvalue $\geq 0$. You could actually set a minimum eigenvalue value other than 0 if you so desire, and you would need to do this if you want to ensure the covariance matrix is nonsingular as you seem to suggest is desired or required. Addition of this semidefinite constraint turns the whole optimization problem into a convex semidefinite program (SDP), or technically, something convertible thereinto.
Formulating and numerically solving such a convex SDP, i.e., objective in your choice of norm (L^p for any $p \geq 1$), plus any objective penalty in your choice of norm $(p \geq 1)$ which need not be the same as the other norm, plus isotonizing (linear) constraints, plus semi-definite constraint, is VERY easy and straightforward using a tool such as CVX http://cvxr.com/cvx/ . This should be very fast executing unless dimension of the covariance matrix (what you called p, not what I called p) is in the thousands or greater. YALMIP http://users.isy.liu.se/johanl/yalmip/ could be used instead of CVX (which only allows formulation and solution of convex optimization problems, except for the optional specification of integer constraints). YALMIP allows for a greater choice of optimization solvers and a greater range of problems (non-convex) which can be formulated and solved than CVX, but has a steeper learning curve.