From 9a876806e1db0cb32e5e63ed697513773c5f45d3 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Thu, 20 Aug 2009 21:04:38 +0200 Subject: [PATCH] use eigen stableNorm() instead of cminpack 'enorm'. The results are mostly slightly better in tests (one test needs 15 iterations intead of 16, for the same result). Some numerical results have improved slightly, too. If one uses blueNorm() instead, an assert for 'overflow' is raised from blueNorm() --- unsupported/Eigen/NonLinear | 7 ++++++ unsupported/Eigen/src/NonLinear/hybrd.h | 12 +++++----- unsupported/Eigen/src/NonLinear/hybrj.h | 12 +++++----- unsupported/Eigen/src/NonLinear/lmder.h | 12 +++++----- unsupported/Eigen/src/NonLinear/lmdif.h | 12 +++++----- unsupported/Eigen/src/NonLinear/lmstr.h | 14 +++++------ unsupported/test/NonLinear.cpp | 32 ++++++++++++++++--------- 7 files changed, 59 insertions(+), 42 deletions(-) diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear index 5606b7b42..a01d705cd 100644 --- a/unsupported/Eigen/NonLinear +++ b/unsupported/Eigen/NonLinear @@ -116,9 +116,16 @@ int lmstr ( minpack_funcderstr_mn fcn, void *p, int m, void chkder ( int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, double *xp, double *fvecp, int mode, double *err ); + +void covar(int n, double *r__, int ldr, const int *ipvt, double tol, double *wa); #endif +template +Scalar ei_enorm ( int n, const Scalar *x ){ + return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).stableNorm(); +// return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).blueNorm(); +} #include "src/NonLinear/lmder1.h" #include "src/NonLinear/lmder.h" diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h index 53019a003..6a28b63e4 100644 --- a/unsupported/Eigen/src/NonLinear/hybrd.h +++ b/unsupported/Eigen/src/NonLinear/hybrd.h @@ -77,7 +77,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(n, &fvec[1]); + fnorm = ei_enorm(n, &fvec[1]); /* determine the number of calls to fcn needed to compute */ /* the jacobian matrix. */ @@ -140,7 +140,7 @@ L50: wa3[j] = diag[j] * x[j]; /* L60: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -250,7 +250,7 @@ L190: wa3[j] = diag[j] * wa1[j]; /* L200: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -265,7 +265,7 @@ L190: if (iflag < 0) { goto L300; } - fnorm1 = enorm(n, &wa4[1]); + fnorm1 = ei_enorm(n, &wa4[1]); /* compute the scaled actual reduction. */ @@ -291,7 +291,7 @@ L190: wa3[i__] = qtf[i__] + sum; /* L220: */ } - temp = enorm(n, &wa3[1]); + temp = ei_enorm(n, &wa3[1]); prered = 0.; if (temp < fnorm) { /* Computing 2nd power */ @@ -344,7 +344,7 @@ L240: fvec[j] = wa4[j]; /* L250: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L260: diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h index 30712e5db..013757ae7 100644 --- a/unsupported/Eigen/src/NonLinear/hybrj.h +++ b/unsupported/Eigen/src/NonLinear/hybrj.h @@ -78,7 +78,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(n, &fvec[1]); + fnorm = ei_enorm(n, &fvec[1]); /* initialize iteration counter and monitors. */ @@ -133,7 +133,7 @@ L50: wa3[j] = diag[j] * x[j]; /* L60: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -243,7 +243,7 @@ L190: wa3[j] = diag[j] * wa1[j]; /* L200: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -258,7 +258,7 @@ L190: if (iflag < 0) { goto L300; } - fnorm1 = enorm(n, &wa4[1]); + fnorm1 = ei_enorm(n, &wa4[1]); /* compute the scaled actual reduction. */ @@ -284,7 +284,7 @@ L190: wa3[i__] = qtf[i__] + sum; /* L220: */ } - temp = enorm(n, &wa3[1]); + temp = ei_enorm(n, &wa3[1]); prered = 0.; if (temp < fnorm) { /* Computing 2nd power */ @@ -337,7 +337,7 @@ L240: fvec[j] = wa4[j]; /* L250: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L260: diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h index 6df4dcbb2..99b912f4b 100644 --- a/unsupported/Eigen/src/NonLinear/lmder.h +++ b/unsupported/Eigen/src/NonLinear/lmder.h @@ -71,7 +71,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -136,7 +136,7 @@ L60: wa3[j] = diag[j] * x[j]; /* L70: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -242,7 +242,7 @@ L200: wa3[j] = diag[j] * wa1[j]; /* L210: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -257,7 +257,7 @@ L200: if (iflag < 0) { goto L300; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -283,7 +283,7 @@ L200: } /* L230: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -351,7 +351,7 @@ L260: fvec[i__] = wa4[i__]; /* L280: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L290: diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h index 136451a09..da04caacf 100644 --- a/unsupported/Eigen/src/NonLinear/lmdif.h +++ b/unsupported/Eigen/src/NonLinear/lmdif.h @@ -73,7 +73,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -139,7 +139,7 @@ L60: wa3[j] = diag[j] * x[j]; /* L70: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -245,7 +245,7 @@ L200: wa3[j] = diag[j] * wa1[j]; /* L210: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -260,7 +260,7 @@ L200: if (iflag < 0) { goto L300; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -286,7 +286,7 @@ L200: } /* L230: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -354,7 +354,7 @@ L260: fvec[i__] = wa4[i__]; /* L280: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L290: diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h index 17aba6865..015e0bd13 100644 --- a/unsupported/Eigen/src/NonLinear/lmstr.h +++ b/unsupported/Eigen/src/NonLinear/lmstr.h @@ -72,7 +72,7 @@ L20: if (iflag < 0) { goto L340; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -136,7 +136,7 @@ L40: sing = TRUE_; } ipvt[j] = j; - wa2[j] = enorm(j, &fjac[j * fjac_dim1 + 1]); + wa2[j] = ei_enorm(j, &fjac[j * fjac_dim1 + 1]); /* L80: */ } if (! sing) { @@ -194,7 +194,7 @@ L150: wa3[j] = diag[j] * x[j]; /* L160: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -269,7 +269,7 @@ L240: wa3[j] = diag[j] * wa1[j]; /* L250: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -284,7 +284,7 @@ L240: if (iflag < 0) { goto L340; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -310,7 +310,7 @@ L240: } /* L270: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -378,7 +378,7 @@ L300: fvec[i__] = wa4[i__]; /* L320: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L330: diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 00ac38795..ff94198c2 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -1142,7 +1142,7 @@ void testNistLanczos1(void) VERIFY( 79 == nfev); VERIFY( 72 == njev); // check norm^2 - VERIFY_IS_APPROX(fvec.squaredNorm(), 1.4292388868910E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats + VERIFY_IS_APPROX(fvec.squaredNorm(), 1.429604433690E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats // check x VERIFY_IS_APPROX(x[0], 9.5100000027E-02 ); VERIFY_IS_APPROX(x[1], 1.0000000001E+00 ); @@ -1294,7 +1294,7 @@ void testNistMGH10(void) info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); // check return value - VERIFY( 3 == info); + VERIFY( 2 == info); VERIFY( 285 == nfev); VERIFY( 250 == njev); // check norm^2 @@ -1312,7 +1312,7 @@ void testNistMGH10(void) info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); // check return value - VERIFY( 3 == info); + VERIFY( 2 == info); VERIFY( 126 == nfev); VERIFY( 116 == njev); // check norm^2 @@ -1388,8 +1388,8 @@ void testNistBoxBOD(void) // check return value VERIFY( 1 == info); - VERIFY( 16 == nfev); - VERIFY( 15 == njev); + VERIFY( 15 == nfev); + VERIFY( 14 == njev); // check norm^2 VERIFY_IS_APPROX(fvec.squaredNorm(), 1.1680088766E+03); // check x @@ -1531,15 +1531,15 @@ void testNistMGH09(void) // check return value VERIFY( 1 == info); - VERIFY( 487== nfev); - VERIFY( 378 == njev); + VERIFY( 503== nfev); + VERIFY( 385 == njev); // check norm^2 VERIFY_IS_APPROX(fvec.squaredNorm(), 3.0750560385E-04); // check x - VERIFY_IS_APPROX(x[0], 0.19280590); // should be 1.9280693458E-01 - VERIFY_IS_APPROX(x[1], 0.19130543); // should be 1.9128232873E-01 - VERIFY_IS_APPROX(x[2], 0.12306085); // should be 1.2305650693E-01 - VERIFY_IS_APPROX(x[3], 0.13607303); // should be 1.3606233068E-01 + VERIFY_IS_APPROX(x[0], 0.19280624); // should be 1.9280693458E-01 + VERIFY_IS_APPROX(x[1], 0.19129774); // should be 1.9128232873E-01 + VERIFY_IS_APPROX(x[2], 0.12305940); // should be 1.2305650693E-01 + VERIFY_IS_APPROX(x[3], 0.13606946); // should be 1.3606233068E-01 /* * Second try @@ -1924,3 +1924,13 @@ void test_NonLinear() CALL_SUBTEST(testNistEckerle4()); } +/* + * Can be useful for debugging... + printf("info, nfev, jev : %d, %d, %d\n", info, nfev, njev); + printf("x[0] : %.32g\n", x[0]); + printf("x[1] : %.32g\n", x[1]); + printf("x[2] : %.32g\n", x[2]); + printf("x[3] : %.32g\n", x[3]); + printf("fvec.squaredNorm() : %.32g\n", fvec.squaredNorm()); +*/ +