Quelle dogleg.h Sprache: C

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();
    eigen_assert(n==qtb.size());
    eigen_assert(n==x.size());
    eigen_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 * numext::abs2(sgnorm / delta) + sqrt(numext::abs2(temp - delta / qnorm) + (1.-numext::abs2(delta / qnorm)) * (1.-numext::abs2(sgnorm / delta)));
    alpha = delta / qnorm * (1. - numext::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

quality80%

¤ Dauer der Verarbeitung: 0.13 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.