HybridNonLinearSolver.h
Go to the documentation of this file.
1 // -*- coding: utf-8
2 // vim: set fileencoding=utf-8
3 
4 // This file is part of Eigen, a lightweight C++ template library
5 // for linear algebra.
6 //
7 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
8 //
9 // This Source Code Form is subject to the terms of the Mozilla
10 // Public License v. 2.0. If a copy of the MPL was not distributed
11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
12 
13 #ifndef EIGEN_HYBRIDNONLINEARSOLVER_H
14 #define EIGEN_HYBRIDNONLINEARSOLVER_H
15 
16 #include "./InternalHeaderCheck.h"
17 
18 namespace Eigen {
19 
20 namespace HybridNonLinearSolverSpace {
21  enum Status {
22  Running = -1,
29  UserAsked = 6
30  };
31 }
32 
44 template<typename FunctorType, typename Scalar=double>
46 {
47 public:
48  typedef DenseIndex Index;
49 
50  HybridNonLinearSolver(FunctorType &_functor)
51  : functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; useExternalScaling=false;}
52 
53  struct Parameters {
55  : factor(Scalar(100.))
56  , maxfev(1000)
57  , xtol(numext::sqrt(NumTraits<Scalar>::epsilon()))
58  , nb_of_subdiagonals(-1)
60  , epsfcn(Scalar(0.)) {}
61  Scalar factor;
62  Index maxfev; // maximum number of function evaluation
63  Scalar xtol;
66  Scalar epsfcn;
67  };
70  /* TODO: if eigen provides a triangular storage, use it here */
72 
74  FVectorType &x,
75  const Scalar tol = numext::sqrt(NumTraits<Scalar>::epsilon())
76  );
77 
81 
83  FVectorType &x,
84  const Scalar tol = numext::sqrt(NumTraits<Scalar>::epsilon())
85  );
86 
90 
99  Scalar fnorm;
101 private:
102  FunctorType &functor;
104  Scalar sum;
105  bool sing;
106  Scalar temp;
107  Scalar delta;
108  bool jeval;
110  Scalar ratio;
111  Scalar pnorm, xnorm, fnorm1;
114  Scalar actred, prered;
116 
118 };
119 
120 
121 
122 template<typename FunctorType, typename Scalar>
125  FVectorType &x,
126  const Scalar tol
127  )
128 {
129  n = x.size();
130 
131  /* check the input parameters for errors. */
132  if (n <= 0 || tol < 0.)
134 
135  resetParameters();
136  parameters.maxfev = 100*(n+1);
137  parameters.xtol = tol;
138  diag.setConstant(n, 1.);
139  useExternalScaling = true;
140  return solve(x);
141 }
142 
143 template<typename FunctorType, typename Scalar>
146 {
147  n = x.size();
148 
149  wa1.resize(n); wa2.resize(n); wa3.resize(n); wa4.resize(n);
150  fvec.resize(n);
151  qtf.resize(n);
152  fjac.resize(n, n);
153  if (!useExternalScaling)
154  diag.resize(n);
155  eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
156 
157  /* Function Body */
158  nfev = 0;
159  njev = 0;
160 
161  /* check the input parameters for errors. */
162  if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
164  if (useExternalScaling)
165  for (Index j = 0; j < n; ++j)
166  if (diag[j] <= 0.)
168 
169  /* evaluate the function at the starting point */
170  /* and calculate its norm. */
171  nfev = 1;
172  if ( functor(x, fvec) < 0)
174  fnorm = fvec.stableNorm();
175 
176  /* initialize iteration counter and monitors. */
177  iter = 1;
178  ncsuc = 0;
179  ncfail = 0;
180  nslow1 = 0;
181  nslow2 = 0;
182 
184 }
185 
186 template<typename FunctorType, typename Scalar>
189 {
190  using std::abs;
191 
192  eigen_assert(x.size()==n); // check the caller is not cheating us
193 
194  Index j;
195  std::vector<JacobiRotation<Scalar> > v_givens(n), w_givens(n);
196 
197  jeval = true;
198 
199  /* calculate the jacobian matrix. */
200  if ( functor.df(x, fjac) < 0)
202  ++njev;
203 
204  wa2 = fjac.colwise().blueNorm();
205 
206  /* on the first iteration and if external scaling is not used, scale according */
207  /* to the norms of the columns of the initial jacobian. */
208  if (iter == 1) {
209  if (!useExternalScaling)
210  for (j = 0; j < n; ++j)
211  diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
212 
213  /* on the first iteration, calculate the norm of the scaled x */
214  /* and initialize the step bound delta. */
215  xnorm = diag.cwiseProduct(x).stableNorm();
216  delta = parameters.factor * xnorm;
217  if (delta == 0.)
218  delta = parameters.factor;
219  }
220 
221  /* compute the qr factorization of the jacobian. */
222  HouseholderQR<JacobianType> qrfac(fjac); // no pivoting:
223 
224  /* copy the triangular factor of the qr factorization into r. */
225  R = qrfac.matrixQR();
226 
227  /* accumulate the orthogonal factor in fjac. */
228  fjac = qrfac.householderQ();
229 
230  /* form (q transpose)*fvec and store in qtf. */
231  qtf = fjac.transpose() * fvec;
232 
233  /* rescale if necessary. */
234  if (!useExternalScaling)
235  diag = diag.cwiseMax(wa2);
236 
237  while (true) {
238  /* determine the direction p. */
239  internal::dogleg<Scalar>(R, diag, qtf, delta, wa1);
240 
241  /* store the direction p and x + p. calculate the norm of p. */
242  wa1 = -wa1;
243  wa2 = x + wa1;
244  pnorm = diag.cwiseProduct(wa1).stableNorm();
245 
246  /* on the first iteration, adjust the initial step bound. */
247  if (iter == 1)
248  delta = (std::min)(delta,pnorm);
249 
250  /* evaluate the function at x + p and calculate its norm. */
251  if ( functor(wa2, wa4) < 0)
253  ++nfev;
254  fnorm1 = wa4.stableNorm();
255 
256  /* compute the scaled actual reduction. */
257  actred = -1.;
258  if (fnorm1 < fnorm) /* Computing 2nd power */
259  actred = 1. - numext::abs2(fnorm1 / fnorm);
260 
261  /* compute the scaled predicted reduction. */
262  wa3 = R.template triangularView<Upper>()*wa1 + qtf;
263  temp = wa3.stableNorm();
264  prered = 0.;
265  if (temp < fnorm) /* Computing 2nd power */
266  prered = 1. - numext::abs2(temp / fnorm);
267 
268  /* compute the ratio of the actual to the predicted reduction. */
269  ratio = 0.;
270  if (prered > 0.)
271  ratio = actred / prered;
272 
273  /* update the step bound. */
274  if (ratio < Scalar(.1)) {
275  ncsuc = 0;
276  ++ncfail;
277  delta = Scalar(.5) * delta;
278  } else {
279  ncfail = 0;
280  ++ncsuc;
281  if (ratio >= Scalar(.5) || ncsuc > 1)
282  delta = (std::max)(delta, pnorm / Scalar(.5));
283  if (abs(ratio - 1.) <= Scalar(.1)) {
284  delta = pnorm / Scalar(.5);
285  }
286  }
287 
288  /* test for successful iteration. */
289  if (ratio >= Scalar(1e-4)) {
290  /* successful iteration. update x, fvec, and their norms. */
291  x = wa2;
292  wa2 = diag.cwiseProduct(x);
293  fvec = wa4;
294  xnorm = wa2.stableNorm();
295  fnorm = fnorm1;
296  ++iter;
297  }
298 
299  /* determine the progress of the iteration. */
300  ++nslow1;
301  if (actred >= Scalar(.001))
302  nslow1 = 0;
303  if (jeval)
304  ++nslow2;
305  if (actred >= Scalar(.1))
306  nslow2 = 0;
307 
308  /* test for convergence. */
309  if (delta <= parameters.xtol * xnorm || fnorm == 0.)
311 
312  /* tests for termination and stringent tolerances. */
313  if (nfev >= parameters.maxfev)
315  if (Scalar(.1) * (std::max)(Scalar(.1) * delta, pnorm) <= NumTraits<Scalar>::epsilon() * xnorm)
317  if (nslow2 == 5)
319  if (nslow1 == 10)
321 
322  /* criterion for recalculating jacobian. */
323  if (ncfail == 2)
324  break; // leave inner loop and go for the next outer loop iteration
325 
326  /* calculate the rank one modification to the jacobian */
327  /* and update qtf if necessary. */
328  wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm );
329  wa2 = fjac.transpose() * wa4;
330  if (ratio >= Scalar(1e-4))
331  qtf = wa2;
332  wa2 = (wa2-wa3)/pnorm;
333 
334  /* compute the qr factorization of the updated jacobian. */
335  internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing);
336  internal::r1mpyq<Scalar>(n, n, fjac.data(), v_givens, w_givens);
337  internal::r1mpyq<Scalar>(1, n, qtf.data(), v_givens, w_givens);
338 
339  jeval = false;
340  }
342 }
343 
344 template<typename FunctorType, typename Scalar>
347 {
348  HybridNonLinearSolverSpace::Status status = solveInit(x);
350  return status;
352  status = solveOneStep(x);
353  return status;
354 }
355 
356 
357 
358 template<typename FunctorType, typename Scalar>
361  FVectorType &x,
362  const Scalar tol
363  )
364 {
365  n = x.size();
366 
367  /* check the input parameters for errors. */
368  if (n <= 0 || tol < 0.)
370 
371  resetParameters();
372  parameters.maxfev = 200*(n+1);
373  parameters.xtol = tol;
374 
375  diag.setConstant(n, 1.);
376  useExternalScaling = true;
377  return solveNumericalDiff(x);
378 }
379 
380 template<typename FunctorType, typename Scalar>
383 {
384  n = x.size();
385 
386  if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals= n-1;
387  if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals= n-1;
388 
389  wa1.resize(n); wa2.resize(n); wa3.resize(n); wa4.resize(n);
390  qtf.resize(n);
391  fjac.resize(n, n);
392  fvec.resize(n);
393  if (!useExternalScaling)
394  diag.resize(n);
395  eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
396 
397  /* Function Body */
398  nfev = 0;
399  njev = 0;
400 
401  /* check the input parameters for errors. */
402  if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
404  if (useExternalScaling)
405  for (Index j = 0; j < n; ++j)
406  if (diag[j] <= 0.)
408 
409  /* evaluate the function at the starting point */
410  /* and calculate its norm. */
411  nfev = 1;
412  if ( functor(x, fvec) < 0)
414  fnorm = fvec.stableNorm();
415 
416  /* initialize iteration counter and monitors. */
417  iter = 1;
418  ncsuc = 0;
419  ncfail = 0;
420  nslow1 = 0;
421  nslow2 = 0;
422 
424 }
425 
426 template<typename FunctorType, typename Scalar>
429 {
430  using std::sqrt;
431  using std::abs;
432 
433  eigen_assert(x.size()==n); // check the caller is not cheating us
434 
435  Index j;
436  std::vector<JacobiRotation<Scalar> > v_givens(n), w_givens(n);
437 
438  jeval = true;
439  if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals= n-1;
440  if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals= n-1;
441 
442  /* calculate the jacobian matrix. */
443  if (internal::fdjac1(functor, x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0)
445  nfev += (std::min)(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n);
446 
447  wa2 = fjac.colwise().blueNorm();
448 
449  /* on the first iteration and if external scaling is not used, scale according */
450  /* to the norms of the columns of the initial jacobian. */
451  if (iter == 1) {
452  if (!useExternalScaling)
453  for (j = 0; j < n; ++j)
454  diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
455 
456  /* on the first iteration, calculate the norm of the scaled x */
457  /* and initialize the step bound delta. */
458  xnorm = diag.cwiseProduct(x).stableNorm();
459  delta = parameters.factor * xnorm;
460  if (delta == 0.)
461  delta = parameters.factor;
462  }
463 
464  /* compute the qr factorization of the jacobian. */
465  HouseholderQR<JacobianType> qrfac(fjac); // no pivoting:
466 
467  /* copy the triangular factor of the qr factorization into r. */
468  R = qrfac.matrixQR();
469 
470  /* accumulate the orthogonal factor in fjac. */
471  fjac = qrfac.householderQ();
472 
473  /* form (q transpose)*fvec and store in qtf. */
474  qtf = fjac.transpose() * fvec;
475 
476  /* rescale if necessary. */
477  if (!useExternalScaling)
478  diag = diag.cwiseMax(wa2);
479 
480  while (true) {
481  /* determine the direction p. */
482  internal::dogleg<Scalar>(R, diag, qtf, delta, wa1);
483 
484  /* store the direction p and x + p. calculate the norm of p. */
485  wa1 = -wa1;
486  wa2 = x + wa1;
487  pnorm = diag.cwiseProduct(wa1).stableNorm();
488 
489  /* on the first iteration, adjust the initial step bound. */
490  if (iter == 1)
491  delta = (std::min)(delta,pnorm);
492 
493  /* evaluate the function at x + p and calculate its norm. */
494  if ( functor(wa2, wa4) < 0)
496  ++nfev;
497  fnorm1 = wa4.stableNorm();
498 
499  /* compute the scaled actual reduction. */
500  actred = -1.;
501  if (fnorm1 < fnorm) /* Computing 2nd power */
502  actred = 1. - numext::abs2(fnorm1 / fnorm);
503 
504  /* compute the scaled predicted reduction. */
505  wa3 = R.template triangularView<Upper>()*wa1 + qtf;
506  temp = wa3.stableNorm();
507  prered = 0.;
508  if (temp < fnorm) /* Computing 2nd power */
509  prered = 1. - numext::abs2(temp / fnorm);
510 
511  /* compute the ratio of the actual to the predicted reduction. */
512  ratio = 0.;
513  if (prered > 0.)
514  ratio = actred / prered;
515 
516  /* update the step bound. */
517  if (ratio < Scalar(.1)) {
518  ncsuc = 0;
519  ++ncfail;
520  delta = Scalar(.5) * delta;
521  } else {
522  ncfail = 0;
523  ++ncsuc;
524  if (ratio >= Scalar(.5) || ncsuc > 1)
525  delta = (std::max)(delta, pnorm / Scalar(.5));
526  if (abs(ratio - 1.) <= Scalar(.1)) {
527  delta = pnorm / Scalar(.5);
528  }
529  }
530 
531  /* test for successful iteration. */
532  if (ratio >= Scalar(1e-4)) {
533  /* successful iteration. update x, fvec, and their norms. */
534  x = wa2;
535  wa2 = diag.cwiseProduct(x);
536  fvec = wa4;
537  xnorm = wa2.stableNorm();
538  fnorm = fnorm1;
539  ++iter;
540  }
541 
542  /* determine the progress of the iteration. */
543  ++nslow1;
544  if (actred >= Scalar(.001))
545  nslow1 = 0;
546  if (jeval)
547  ++nslow2;
548  if (actred >= Scalar(.1))
549  nslow2 = 0;
550 
551  /* test for convergence. */
552  if (delta <= parameters.xtol * xnorm || fnorm == 0.)
554 
555  /* tests for termination and stringent tolerances. */
556  if (nfev >= parameters.maxfev)
558  if (Scalar(.1) * (std::max)(Scalar(.1) * delta, pnorm) <= NumTraits<Scalar>::epsilon() * xnorm)
560  if (nslow2 == 5)
562  if (nslow1 == 10)
564 
565  /* criterion for recalculating jacobian. */
566  if (ncfail == 2)
567  break; // leave inner loop and go for the next outer loop iteration
568 
569  /* calculate the rank one modification to the jacobian */
570  /* and update qtf if necessary. */
571  wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm );
572  wa2 = fjac.transpose() * wa4;
573  if (ratio >= Scalar(1e-4))
574  qtf = wa2;
575  wa2 = (wa2-wa3)/pnorm;
576 
577  /* compute the qr factorization of the updated jacobian. */
578  internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing);
579  internal::r1mpyq<Scalar>(n, n, fjac.data(), v_givens, w_givens);
580  internal::r1mpyq<Scalar>(1, n, qtf.data(), v_givens, w_givens);
581 
582  jeval = false;
583  }
585 }
586 
587 template<typename FunctorType, typename Scalar>
590 {
591  HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x);
593  return status;
595  status = solveNumericalDiffOneStep(x);
596  return status;
597 }
598 
599 } // end namespace Eigen
600 
601 #endif // EIGEN_HYBRIDNONLINEARSOLVER_H
602 
603 //vim: ai ts=4 sts=4 et sw=4
int n
Array< double, 1, 3 > e(1./3., 0.5, 2.)
#define eigen_assert(x)
HouseholderSequenceType householderQ() const
const MatrixType & matrixQR() const
TransposeReturnType transpose() const
Finds a zero of a system of n nonlinear functions in n variables by a modification of the Powell hybr...
HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x)
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x)
HybridNonLinearSolverSpace::Status solveInit(FVectorType &x)
HybridNonLinearSolver & operator=(const HybridNonLinearSolver &)
HybridNonLinearSolverSpace::Status solve(FVectorType &x)
Matrix< Scalar, Dynamic, Dynamic > JacobianType
HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x)
Matrix< Scalar, Dynamic, 1 > FVectorType
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x)
HybridNonLinearSolverSpace::Status hybrj1(FVectorType &x, const Scalar tol=numext::sqrt(NumTraits< Scalar >::epsilon()))
HybridNonLinearSolver(FunctorType &_functor)
HybridNonLinearSolverSpace::Status hybrd1(FVectorType &x, const Scalar tol=numext::sqrt(NumTraits< Scalar >::epsilon()))
Matrix< Scalar, Dynamic, Dynamic > UpperTriangularType
DenseIndex fdjac1(const FunctorType &Functor, Matrix< Scalar, Dynamic, 1 > &x, Matrix< Scalar, Dynamic, 1 > &fvec, Matrix< Scalar, Dynamic, Dynamic > &fjac, DenseIndex ml, DenseIndex mu, Scalar epsfcn)
Definition: fdjac1.h:8
bool abs2(bool x)
EIGEN_ALWAYS_INLINE double sqrt(const double &x)
: TensorContractionSycl.h, provides various tensor contraction kernel for SYCL backend
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Eigen::internal::remove_all_t< DerType >, typename Eigen::internal::traits< Eigen::internal::remove_all_t< DerType >>::Scalar, product) > sqrt(const Eigen::AutoDiffScalar< DerType > &x)
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)
CleanedUpDerType< DerType >::type() max(const AutoDiffScalar< DerType > &x, const T &y)
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
adouble abs(const adouble &x)
Definition: AdolcForward:74
std::ptrdiff_t j