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