@@ -79,24 +79,27 @@ double L_BFGS::ChooseScalingFactor(const size_t iterationNum,
7979{
8080 typedef typename CubeType::elem_type CubeElemType;
8181
82+ constexpr const CubeElemType tol =
83+ 100 * std::numeric_limits<CubeElemType>::epsilon ();
84+
8285 double scalingFactor;
8386 if (iterationNum > 0 )
8487 {
8588 int previousPos = (iterationNum - 1 ) % numBasis;
8689 // Get s and y matrices once instead of multiple times.
8790 const arma::Mat<CubeElemType>& sMat = s.slice (previousPos);
8891 const arma::Mat<CubeElemType>& yMat = y.slice (previousPos);
89-
92+
9093 const CubeElemType tmp = arma::dot (yMat, yMat);
91- const CubeElemType denom = (tmp != CubeElemType ( 0 ) ) ? tmp : CubeElemType (1 );
92-
94+ const CubeElemType denom = (tmp >= tol ) ? tmp : CubeElemType (1 );
95+
9396 scalingFactor = arma::dot (sMat , yMat) / denom;
9497 }
9598 else
9699 {
97100 const CubeElemType tmp = arma::norm (gradient, " fro" );
98-
99- scalingFactor = (tmp != CubeElemType ( 0 ) ) ? (1.0 / tmp) : 1.0 ;
101+
102+ scalingFactor = (tmp >= tol ) ? (1.0 / tmp) : 1.0 ;
100103 }
101104
102105 return scalingFactor;
@@ -135,16 +138,18 @@ void L_BFGS::SearchDirection(const MatType& gradient,
135138 for (size_t i = iterationNum; i != limit; i--)
136139 {
137140 int translatedPosition = (i + (numBasis - 1 )) % numBasis;
138-
141+
139142 const arma::Mat<CubeElemType>& sMat = s.slice (translatedPosition);
140143 const arma::Mat<CubeElemType>& yMat = y.slice (translatedPosition);
141-
144+
142145 const CubeElemType tmp = arma::dot (yMat, sMat );
143-
144- rho[iterationNum - i] = (tmp != CubeElemType (0 )) ? (1.0 / tmp) : CubeElemType (1 );
145-
146- alpha[iterationNum - i] = rho[iterationNum - i] * arma::dot (sMat , searchDirection);
147-
146+
147+ rho[iterationNum - i] = (tmp != CubeElemType (0 )) ? (1.0 / tmp) :
148+ CubeElemType (1 );
149+
150+ alpha[iterationNum - i] = rho[iterationNum - i] *
151+ arma::dot (sMat , searchDirection);
152+
148153 searchDirection -= alpha[iterationNum - i] * yMat;
149154 }
150155
@@ -410,7 +415,8 @@ L_BFGS::Optimize(FunctionType& function,
410415 //
411416 // But don't do this on the first iteration to ensure we always take at
412417 // least one descent step.
413- // TODO: to speed this up, investigate use of arma::norm2est() in Armadillo 12.4
418+ // TODO: to speed this up, investigate use of arma::norm2est() in Armadillo
419+ // 12.4
414420 if (arma::norm (gradient, 2 ) < minGradientNorm)
415421 {
416422 Info << " L-BFGS gradient norm too small (terminating successfully)."
@@ -442,7 +448,7 @@ L_BFGS::Optimize(FunctionType& function,
442448 << std::endl;
443449 break ;
444450 }
445-
451+
446452 // Build an approximation to the Hessian and choose the search
447453 // direction for the current iteration.
448454 SearchDirection (gradient, itNum, scalingFactor, s, y, searchDirection);
0 commit comments