src/NonLinearOptimization/dogleg.h

0001 namespace Eigen {
0002
0003 namespace internal {
0004
0005 template <typename Scalar>
0006 void dogleg(
0007         const Matrix< Scalar, Dynamic, Dynamic >  &qrfac,
0008         const Matrix< Scalar, Dynamic, 1 >  &diag,
0009         const Matrix< Scalar, Dynamic, 1 >  &qtb,
0010         Scalar delta,
0011         Matrix< Scalar, Dynamic, 1 >  &x)
0012 {
0013     using std::abs;
0014     using std::sqrt;
0015
0016     typedef DenseIndex Index;
0017
0018     /* Local variables */
0019     Index i, j;
0020     Scalar sum, temp, alpha, bnorm;
0021     Scalar gnorm, qnorm;
0022     Scalar sgnorm;
0023
0024     /* Function Body */
0025     const Scalar epsmch = NumTraits<Scalar>::epsilon();
0026     const Index n = qrfac.cols();
0027     eigen_assert(n==qtb.size());
0028     eigen_assert(n==x.size());
0029     eigen_assert(n==diag.size());
0030     Matrix< Scalar, Dynamic, 1 >  wa1(n), wa2(n);
0031
0032     /* first, calculate the gauss-newton direction. */
0033     for (j = n-1; j >=0; --j) {
0034         temp = qrfac(j,j);
0035         if (temp == 0.) {
0036             temp = epsmch * qrfac.col(j).head(j+1).maxCoeff();
0037             if (temp == 0.)
0038                 temp = epsmch;
0039         }
0040         if (j==n-1)
0041             x[j] = qtb[j] / temp;
0042         else
0043             x[j] = (qtb[j] - qrfac.row(j).tail(n-j-1).dot(x.tail(n-j-1))) / temp;
0044     }
0045
0046     /* test whether the gauss-newton direction is acceptable. */
0047     qnorm = diag.cwiseProduct(x).stableNorm();
0048     if (qnorm <= delta)
0049         return;
0050
0051     // TODO : this path is not tested by Eigen unit tests
0052
0053     /* the gauss-newton direction is not acceptable. */
0054     /* next, calculate the scaled gradient direction. */
0055
0056     wa1.fill(0.);
0057     for (j = 0; j < n; ++j) {
0058         wa1.tail(n-j) += qrfac.row(j).tail(n-j) * qtb[j];
0059         wa1[j] /= diag[j];
0060     }
0061
0062     /* calculate the norm of the scaled gradient and test for */
0063     /* the special case in which the scaled gradient is zero. */
0064     gnorm = wa1.stableNorm();
0065     sgnorm = 0.;
0066     alpha = delta / qnorm;
0067     if (gnorm == 0.)
0068         goto algo_end;
0069
0070     /* calculate the point along the scaled gradient */
0071     /* at which the quadratic is minimized. */
0072     wa1.array() /= (diag*gnorm).array();
0073     // TODO : once unit tests cover this part,:
0074     // wa2 = qrfac.template triangularView<Upper>() * wa1;
0075     for (j = 0; j < n; ++j) {
0076         sum = 0.;
0077         for (i = j; i < n; ++i) {
0078             sum += qrfac(j,i) * wa1[i];
0079         }
0080         wa2[j] = sum;
0081     }
0082     temp = wa2.stableNorm();
0083     sgnorm = gnorm / temp / temp;
0084
0085     /* test whether the scaled gradient direction is acceptable. */
0086     alpha = 0.;
0087     if (sgnorm >= delta)
0088         goto algo_end;
0089
0090     /* the scaled gradient direction is not acceptable. */
0091     /* finally, calculate the point along the dogleg */
0092     /* at which the quadratic is minimized. */
0093     bnorm = qtb.stableNorm();
0094     temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
0095     temp = temp - delta / qnorm * numext::abs2(sgnorm / delta) + sqrt(numext::abs2(temp - delta / qnorm) + (1.-numext::abs2(delta / qnorm)) * (1.-numext::abs2(sgnorm / delta)));
0096     alpha = delta / qnorm * (1. - numext::abs2(sgnorm / delta)) / temp;
0097 algo_end:
0098
0099     /* form appropriate convex combination of the gauss-newton */
0100     /* direction and the scaled gradient direction. */
0101     temp = (1.-alpha) * (std::min)(sgnorm,delta);
0102     x = temp * wa1 + alpha * x;
0103 }
0104
0105 } // end namespace internal
0106
0107 } // end namespace Eigen