Newton optimizer: handle flat directions without producing NaN

[SteveBronder]

Submission Checklist

• Run unit tests: ./runTests.py src/test/unit
• Run cpplint: make cpplint
• Declare copyright holder and open-source license: see below

Summary

Fixes #3425.

make_negative_definite_and_solve divided the gradient projection by fabs(eigenvalue) with no guard against zero or near zero values. This could lead to some directions of the gradient and hessian being flat and causing NaN values to return. Now we check that the absolute of the eigen value is greater than a tolerance defined by the an epsilon scaled by the overall maximum eigenvalue. We reject non-finite step directions, candidate points, and objective values in newton_step instead of accepting them. And the newton service layer not returns error_codes::SOFTWARE with TERM_LSFAIL when the final log density or parameters are not finite.

Adds a flat_target test model plus unit tests at the solve, step, and service layers that reproduce the reported NaN.

Documentation

Updated docs for make_negative_definite_and_solve to reflect the change.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company): Steve Bronder

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[WardBrian]

Seems similar to #3309 which @nhuurre had some thoughts on, might be a good reviewer

[nhuurre]

Where does this number come from? Is machine epsilon the appropriate baseline here? Why does it depend on the number of rows? Maybe this is related to the numerical precision of the solver...

Also you should handle the case where all eigenvalues are zero. Set some minimum absolute tolerance.

[SteveBronder]

So I'm going to be updating this to be more inline with what Eigen has.

double max_abs_eigenvalue = eigenvalues.cwiseAbs().maxCoeff();
double tolerance = std::fmax(max_abs_eigenvalue * 4 * H.rows()
                                 * std::numeric_limits<double>::epsilon(),
                             std::numeric_limits<double>::min());

I had claude do an explainer on where it got those values from and it gave me the below. I think it is easier to block quote claude here as I think it actually does a decent job of explaining where it got each piece of this.

---

This is the same cutoff Eigen's rank-revealing decompositions use to decide
whether a pivot or singular value is numerically zero, written out in one
expression. Eigen computes it in two steps: a dimensionless threshold(),
and a per-call multiplication by the largest pivot or singular value.

The bound. Eigen master's threshold() in
Eigen/src/misc/RankRevealingBase.h
reads:

// Higham's backward error bound: ||ΔA||₂ ≤ c·min(m,n)·u·||A||₂.
// The factor of 4 covers the constant c.
return m_usePrescribedThreshold
           ? m_prescribedThreshold
           : NumTraits<Scalar>::epsilon() * RealScalar(4 * (std::min)(self().rows(), self().cols()));

That is the backward error bound for a floating-point factorization from
N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed.,
SIAM, 2002. A perturbation of the input smaller than this bound is
indistinguishable from rounding error in the factorization itself, so a
pivot below it carries no information about rank. Eigen 3.4 uses the same
formula without the factor of 4; the comment there
(Eigen/src/QR/ColPivHouseholderQR.h:382-383) says it "turns out to be
identical to Higham's formula used already in LDLt."

Piece by piece:

Term in our code	Term in the bound	Where Eigen does it
max_abs_eigenvalue	‖A‖₂	Eigen multiplies threshold() by the largest pivot or singular value at the point of use: abs(m_maxpivot) * threshold() at Eigen/src/QR/ColPivHouseholderQR.h:259 and Eigen/src/LU/FullPivLU.h:334; m_singularValues.coeff(0) * threshold() at Eigen/src/SVD/SVDBase.h:153. For a symmetric matrix the 2-norm is the largest absolute eigenvalue, so eigenvalues.cwiseAbs().maxCoeff() is the same quantity.
4	the constant c	Eigen master, RankRevealingBase.h: "The factor of 4 covers the constant c."
H.rows()	min(m, n)	Eigen master uses (std::min)(rows(), cols()); 3.4 uses m_qr.diagonalSize(). The Hessian is square, so both equal H.rows().
std::numeric_limits<double>::epsilon()	u, the unit roundoff	Eigen uses NumTraits<Scalar>::epsilon(), which for double is std::numeric_limits<double>::epsilon().
std::fmax(..., std::numeric_limits<double>::min())	lower clamp	Eigen/src/SVD/SVDBase.h:153: numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)()). This is what makes an all-zero Hessian produce a positive cutoff instead of zero.
abs_eigenvalue <= tolerance drops the direction	comparison	Eigen counts a pivot toward the rank only when abs(pivot) > premultiplied_threshold (ColPivHouseholderQR.h:262), so <= is treated as zero, matching ours.

Precedent outside Eigen. The same relative cutoff, dimension times
epsilon times the largest singular value or eigenvalue, is the default in
SciPy's pinvh
(rtol = max(a.shape) * eps, applied to the largest absolute eigenvalue,
which is the symmetric eigendecomposition case we have here),
MATLAB's pinv
(max(size(A)) * eps(norm(A))), and LAPACK's
DGELSD
(singular values S(i) <= RCOND * S(1) are treated as zero, with machine
precision used when RCOND < 0).

---

I wish there was an easy way to share claude sessions.

[nhuurre]

If target is flat then gradient is zero and what you substitute for the eigenvalue doesn't matter, as long as it's finite. But it's also possible that the target is linear, and if so, the gradient is nonzero while the hessian is still zero. I think you'd want nonzero movement in that case. So instead of zero you should use inverse tolerance. Unlike pseudo-inverse, such "saturating inverse" is continuous.

Suggested change

	eigenprojections[i] = 0;
	eigenprojections[i] = -eigenprojections[i] / tolerance;

[nhuurre]

This looks like a reimplementation of stan::math::is_scal_finite

[nhuurre]

By the way, and this is a pre-existing issue, but even a simple model like

parameters {
  vector[3] x;
}
model {
  x ~ normal(0, [10,1,0.1]');
}

takes 100 iterations, which is completely unreasonable behavior for a Newton solver. Apparently step size is 1.90735e-06 for every iteration so something is very wrong.

And that something is stan::model::grad_hess_log_prob. I didn't look into the details, but I do know that a finite-difference algorithm should divide by epsilon, like, at some point, and this one never does.
If I change these "multiply-by-half_epsilon" to "divide-by-half_epsilon" (which is the smallest change that makes the algorithm look like it could be correct), stepsize recovers to 1 and the solver converges in "only" 37 iterations.

stan::model::grad_hess_log_prob is not used by anything else since Laplace sampler gets the Hessian from stan::math::internal::finite_diff_hessian_auto. (I would have thought that that internal namespace meant it's not used outside of math but whatever)
If I change Newton solver to use stan::math::internal::finite_diff_hessian_auto then it converges in 2 iterations, with stepsize 1 for both. This is how Newton should behave on a multinormal target.

[WardBrian]

stan::model::grad_hess_log_prob is not used by anything else since Laplace sampler gets the Hessian from stan::math::internal::finite_diff_hessian_auto. (I would have thought that that internal namespace meant it's not used outside of math but whatever)

That ship may have sailed... Even within math, it's not used by any other code, so I'm not sure why it (or finite_diff_hessian_times_vector_auto, which was just following the pattern) is in internal

If I change Newton solver to use stan::math::internal::finite_diff_hessian_auto then it converges in 2 iterations, with stepsize 1 for both. This is how Newton should behave on a multinormal target.

This sounds like it would be worthy of it's own PR

[SteveBronder]

stan::model::grad_hess_log_prob is not used by anything else since Laplace sampler gets the Hessian from stan::math::internal::finite_diff_hessian_auto. (I would have thought that that internal namespace meant it's not used outside of math but whatever)

That ship may have sailed... Even within math, it's not used by any other code, so I'm not sure why it (or finite_diff_hessian_times_vector_auto, which was just following the pattern) is in internal

For the record, anything inside of math's internal namespace has zero API guarantee aka the math library can change these on a whim with no version notice. This is what Eigen does as well (and partly why upgrading Eigen is such a hassle for us)

[WardBrian]

Yeah I understand that. These functions were added entirely on their own with no direct usages in math, so I suspect the internal was more due to a lack of confidence or something to discourage use rather than an API stability concern