| Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 1 | namespace Eigen { |
| 2 | |
| 3 | namespace internal { |
| 4 | |
| 5 | template <typename Scalar> |
| 6 | void dogleg( |
| 7 | const Matrix< Scalar, Dynamic, Dynamic > &qrfac, |
| 8 | const Matrix< Scalar, Dynamic, 1 > &diag, |
| 9 | const Matrix< Scalar, Dynamic, 1 > &qtb, |
| 10 | Scalar delta, |
| 11 | Matrix< Scalar, Dynamic, 1 > &x) |
| 12 | { |
| 13 | using std::abs; |
| 14 | using std::sqrt; |
| 15 | |
| 16 | typedef DenseIndex Index; |
| 17 | |
| 18 | /* Local variables */ |
| 19 | Index i, j; |
| 20 | Scalar sum, temp, alpha, bnorm; |
| 21 | Scalar gnorm, qnorm; |
| 22 | Scalar sgnorm; |
| 23 | |
| 24 | /* Function Body */ |
| 25 | const Scalar epsmch = NumTraits<Scalar>::epsilon(); |
| 26 | const Index n = qrfac.cols(); |
| 27 | eigen_assert(n==qtb.size()); |
| 28 | eigen_assert(n==x.size()); |
| 29 | eigen_assert(n==diag.size()); |
| 30 | Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n); |
| 31 | |
| 32 | /* first, calculate the gauss-newton direction. */ |
| 33 | for (j = n-1; j >=0; --j) { |
| 34 | temp = qrfac(j,j); |
| 35 | if (temp == 0.) { |
| 36 | temp = epsmch * qrfac.col(j).head(j+1).maxCoeff(); |
| 37 | if (temp == 0.) |
| 38 | temp = epsmch; |
| 39 | } |
| 40 | if (j==n-1) |
| 41 | x[j] = qtb[j] / temp; |
| 42 | else |
| 43 | x[j] = (qtb[j] - qrfac.row(j).tail(n-j-1).dot(x.tail(n-j-1))) / temp; |
| 44 | } |
| 45 | |
| 46 | /* test whether the gauss-newton direction is acceptable. */ |
| 47 | qnorm = diag.cwiseProduct(x).stableNorm(); |
| 48 | if (qnorm <= delta) |
| 49 | return; |
| 50 | |
| 51 | // TODO : this path is not tested by Eigen unit tests |
| 52 | |
| 53 | /* the gauss-newton direction is not acceptable. */ |
| 54 | /* next, calculate the scaled gradient direction. */ |
| 55 | |
| 56 | wa1.fill(0.); |
| 57 | for (j = 0; j < n; ++j) { |
| 58 | wa1.tail(n-j) += qrfac.row(j).tail(n-j) * qtb[j]; |
| 59 | wa1[j] /= diag[j]; |
| 60 | } |
| 61 | |
| 62 | /* calculate the norm of the scaled gradient and test for */ |
| 63 | /* the special case in which the scaled gradient is zero. */ |
| 64 | gnorm = wa1.stableNorm(); |
| 65 | sgnorm = 0.; |
| 66 | alpha = delta / qnorm; |
| 67 | if (gnorm == 0.) |
| 68 | goto algo_end; |
| 69 | |
| 70 | /* calculate the point along the scaled gradient */ |
| 71 | /* at which the quadratic is minimized. */ |
| 72 | wa1.array() /= (diag*gnorm).array(); |
| 73 | // TODO : once unit tests cover this part,: |
| 74 | // wa2 = qrfac.template triangularView<Upper>() * wa1; |
| 75 | for (j = 0; j < n; ++j) { |
| 76 | sum = 0.; |
| 77 | for (i = j; i < n; ++i) { |
| 78 | sum += qrfac(j,i) * wa1[i]; |
| 79 | } |
| 80 | wa2[j] = sum; |
| 81 | } |
| 82 | temp = wa2.stableNorm(); |
| 83 | sgnorm = gnorm / temp / temp; |
| 84 | |
| 85 | /* test whether the scaled gradient direction is acceptable. */ |
| 86 | alpha = 0.; |
| 87 | if (sgnorm >= delta) |
| 88 | goto algo_end; |
| 89 | |
| 90 | /* the scaled gradient direction is not acceptable. */ |
| 91 | /* finally, calculate the point along the dogleg */ |
| 92 | /* at which the quadratic is minimized. */ |
| 93 | bnorm = qtb.stableNorm(); |
| 94 | temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta); |
| 95 | temp = temp - delta / qnorm * numext::abs2(sgnorm / delta) + sqrt(numext::abs2(temp - delta / qnorm) + (1.-numext::abs2(delta / qnorm)) * (1.-numext::abs2(sgnorm / delta))); |
| 96 | alpha = delta / qnorm * (1. - numext::abs2(sgnorm / delta)) / temp; |
| 97 | algo_end: |
| 98 | |
| 99 | /* form appropriate convex combination of the gauss-newton */ |
| 100 | /* direction and the scaled gradient direction. */ |
| 101 | temp = (1.-alpha) * (std::min)(sgnorm,delta); |
| 102 | x = temp * wa1 + alpha * x; |
| 103 | } |
| 104 | |
| 105 | } // end namespace internal |
| 106 | |
| 107 | } // end namespace Eigen |