| #include "./InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| namespace internal { |
| |
| template <typename Scalar> |
| void lmpar( |
| Matrix< Scalar, Dynamic, Dynamic > &r, |
| const VectorXi &ipvt, |
| const Matrix< Scalar, Dynamic, 1 > &diag, |
| const Matrix< Scalar, Dynamic, 1 > &qtb, |
| Scalar delta, |
| Scalar &par, |
| Matrix< Scalar, Dynamic, 1 > &x) |
| { |
| using std::abs; |
| using std::sqrt; |
| typedef DenseIndex Index; |
| |
| /* Local variables */ |
| Index i, j, l; |
| Scalar fp; |
| Scalar parc, parl; |
| Index iter; |
| Scalar temp, paru; |
| Scalar gnorm; |
| Scalar dxnorm; |
| |
| |
| /* Function Body */ |
| const Scalar dwarf = (std::numeric_limits<Scalar>::min)(); |
| const Index n = r.cols(); |
| eigen_assert(n==diag.size()); |
| eigen_assert(n==qtb.size()); |
| eigen_assert(n==x.size()); |
| |
| Matrix< Scalar, Dynamic, 1 > wa1, wa2; |
| |
| /* compute and store in x the gauss-newton direction. if the */ |
| /* jacobian is rank-deficient, obtain a least squares solution. */ |
| Index nsing = n-1; |
| wa1 = qtb; |
| for (j = 0; j < n; ++j) { |
| if (r(j,j) == 0. && nsing == n-1) |
| nsing = j - 1; |
| if (nsing < n-1) |
| wa1[j] = 0.; |
| } |
| for (j = nsing; j>=0; --j) { |
| wa1[j] /= r(j,j); |
| temp = wa1[j]; |
| for (i = 0; i < j ; ++i) |
| wa1[i] -= r(i,j) * temp; |
| } |
| |
| for (j = 0; j < n; ++j) |
| x[ipvt[j]] = wa1[j]; |
| |
| /* initialize the iteration counter. */ |
| /* evaluate the function at the origin, and test */ |
| /* for acceptance of the gauss-newton direction. */ |
| iter = 0; |
| wa2 = diag.cwiseProduct(x); |
| dxnorm = wa2.blueNorm(); |
| fp = dxnorm - delta; |
| if (fp <= Scalar(0.1) * delta) { |
| par = 0; |
| return; |
| } |
| |
| /* if the jacobian is not rank deficient, the newton */ |
| /* step provides a lower bound, parl, for the zero of */ |
| /* the function. otherwise set this bound to zero. */ |
| parl = 0.; |
| if (nsing >= n-1) { |
| for (j = 0; j < n; ++j) { |
| l = ipvt[j]; |
| wa1[j] = diag[l] * (wa2[l] / dxnorm); |
| } |
| // it's actually a triangularView.solveInplace(), though in a weird |
| // way: |
| for (j = 0; j < n; ++j) { |
| Scalar sum = 0.; |
| for (i = 0; i < j; ++i) |
| sum += r(i,j) * wa1[i]; |
| wa1[j] = (wa1[j] - sum) / r(j,j); |
| } |
| temp = wa1.blueNorm(); |
| parl = fp / delta / temp / temp; |
| } |
| |
| /* calculate an upper bound, paru, for the zero of the function. */ |
| for (j = 0; j < n; ++j) |
| wa1[j] = r.col(j).head(j+1).dot(qtb.head(j+1)) / diag[ipvt[j]]; |
| |
| gnorm = wa1.stableNorm(); |
| paru = gnorm / delta; |
| if (paru == 0.) |
| paru = dwarf / (std::min)(delta,Scalar(0.1)); |
| |
| /* if the input par lies outside of the interval (parl,paru), */ |
| /* set par to the closer endpoint. */ |
| par = (std::max)(par,parl); |
| par = (std::min)(par,paru); |
| if (par == 0.) |
| par = gnorm / dxnorm; |
| |
| /* beginning of an iteration. */ |
| while (true) { |
| ++iter; |
| |
| /* evaluate the function at the current value of par. */ |
| if (par == 0.) |
| par = (std::max)(dwarf,Scalar(.001) * paru); /* Computing MAX */ |
| wa1 = sqrt(par)* diag; |
| |
| Matrix< Scalar, Dynamic, 1 > sdiag(n); |
| qrsolv<Scalar>(r, ipvt, wa1, qtb, x, sdiag); |
| |
| wa2 = diag.cwiseProduct(x); |
| dxnorm = wa2.blueNorm(); |
| temp = fp; |
| fp = dxnorm - delta; |
| |
| /* if the function is small enough, accept the current value */ |
| /* of par. also test for the exceptional cases where parl */ |
| /* is zero or the number of iterations has reached 10. */ |
| if (abs(fp) <= Scalar(0.1) * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10) |
| break; |
| |
| /* compute the newton correction. */ |
| for (j = 0; j < n; ++j) { |
| l = ipvt[j]; |
| wa1[j] = diag[l] * (wa2[l] / dxnorm); |
| } |
| for (j = 0; j < n; ++j) { |
| wa1[j] /= sdiag[j]; |
| temp = wa1[j]; |
| for (i = j+1; i < n; ++i) |
| wa1[i] -= r(i,j) * temp; |
| } |
| temp = wa1.blueNorm(); |
| parc = fp / delta / temp / temp; |
| |
| /* depending on the sign of the function, update parl or paru. */ |
| if (fp > 0.) |
| parl = (std::max)(parl,par); |
| if (fp < 0.) |
| paru = (std::min)(paru,par); |
| |
| /* compute an improved estimate for par. */ |
| /* Computing MAX */ |
| par = (std::max)(parl,par+parc); |
| |
| /* end of an iteration. */ |
| } |
| |
| /* termination. */ |
| if (iter == 0) |
| par = 0.; |
| return; |
| } |
| |
| template <typename Scalar> |
| void lmpar2( |
| const ColPivHouseholderQR<Matrix< Scalar, Dynamic, Dynamic> > &qr, |
| const Matrix< Scalar, Dynamic, 1 > &diag, |
| const Matrix< Scalar, Dynamic, 1 > &qtb, |
| Scalar delta, |
| Scalar &par, |
| Matrix< Scalar, Dynamic, 1 > &x) |
| |
| { |
| using std::sqrt; |
| using std::abs; |
| typedef DenseIndex Index; |
| |
| /* Local variables */ |
| Index j; |
| Scalar fp; |
| Scalar parc, parl; |
| Index iter; |
| Scalar temp, paru; |
| Scalar gnorm; |
| Scalar dxnorm; |
| |
| |
| /* Function Body */ |
| const Scalar dwarf = (std::numeric_limits<Scalar>::min)(); |
| const Index n = qr.matrixQR().cols(); |
| eigen_assert(n==diag.size()); |
| eigen_assert(n==qtb.size()); |
| |
| Matrix< Scalar, Dynamic, 1 > wa1, wa2; |
| |
| /* compute and store in x the gauss-newton direction. if the */ |
| /* jacobian is rank-deficient, obtain a least squares solution. */ |
| |
| // const Index rank = qr.nonzeroPivots(); // exactly double(0.) |
| const Index rank = qr.rank(); // use a threshold |
| wa1 = qtb; |
| wa1.tail(n-rank).setZero(); |
| qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solveInPlace(wa1.head(rank)); |
| |
| x = qr.colsPermutation()*wa1; |
| |
| /* initialize the iteration counter. */ |
| /* evaluate the function at the origin, and test */ |
| /* for acceptance of the gauss-newton direction. */ |
| iter = 0; |
| wa2 = diag.cwiseProduct(x); |
| dxnorm = wa2.blueNorm(); |
| fp = dxnorm - delta; |
| if (fp <= Scalar(0.1) * delta) { |
| par = 0; |
| return; |
| } |
| |
| /* if the jacobian is not rank deficient, the newton */ |
| /* step provides a lower bound, parl, for the zero of */ |
| /* the function. otherwise set this bound to zero. */ |
| parl = 0.; |
| if (rank==n) { |
| wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2)/dxnorm; |
| qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1); |
| temp = wa1.blueNorm(); |
| parl = fp / delta / temp / temp; |
| } |
| |
| /* calculate an upper bound, paru, for the zero of the function. */ |
| for (j = 0; j < n; ++j) |
| wa1[j] = qr.matrixQR().col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)]; |
| |
| gnorm = wa1.stableNorm(); |
| paru = gnorm / delta; |
| if (paru == 0.) |
| paru = dwarf / (std::min)(delta,Scalar(0.1)); |
| |
| /* if the input par lies outside of the interval (parl,paru), */ |
| /* set par to the closer endpoint. */ |
| par = (std::max)(par,parl); |
| par = (std::min)(par,paru); |
| if (par == 0.) |
| par = gnorm / dxnorm; |
| |
| /* beginning of an iteration. */ |
| Matrix< Scalar, Dynamic, Dynamic > s = qr.matrixQR(); |
| while (true) { |
| ++iter; |
| |
| /* evaluate the function at the current value of par. */ |
| if (par == 0.) |
| par = (std::max)(dwarf,Scalar(.001) * paru); /* Computing MAX */ |
| wa1 = sqrt(par)* diag; |
| |
| Matrix< Scalar, Dynamic, 1 > sdiag(n); |
| qrsolv<Scalar>(s, qr.colsPermutation().indices(), wa1, qtb, x, sdiag); |
| |
| wa2 = diag.cwiseProduct(x); |
| dxnorm = wa2.blueNorm(); |
| temp = fp; |
| fp = dxnorm - delta; |
| |
| /* if the function is small enough, accept the current value */ |
| /* of par. also test for the exceptional cases where parl */ |
| /* is zero or the number of iterations has reached 10. */ |
| if (abs(fp) <= Scalar(0.1) * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10) |
| break; |
| |
| /* compute the newton correction. */ |
| wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2/dxnorm); |
| // we could almost use this here, but the diagonal is outside qr, in sdiag[] |
| // qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1); |
| for (j = 0; j < n; ++j) { |
| wa1[j] /= sdiag[j]; |
| temp = wa1[j]; |
| for (Index i = j+1; i < n; ++i) |
| wa1[i] -= s(i,j) * temp; |
| } |
| temp = wa1.blueNorm(); |
| parc = fp / delta / temp / temp; |
| |
| /* depending on the sign of the function, update parl or paru. */ |
| if (fp > 0.) |
| parl = (std::max)(parl,par); |
| if (fp < 0.) |
| paru = (std::min)(paru,par); |
| |
| /* compute an improved estimate for par. */ |
| par = (std::max)(parl,par+parc); |
| } |
| if (iter == 0) |
| par = 0.; |
| return; |
| } |
| |
| } // end namespace internal |
| |
| } // end namespace Eigen |
