|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
High Performance implementation of the Feasible Newton algorythm. More...
#include <FeasibleNewton.hpp>
Inheritance diagram for MBMesquite::FeasibleNewton: Collaboration diagram for MBMesquite::FeasibleNewton:Public Member Functions | |
| MESQUITE_EXPORT | FeasibleNewton (ObjectiveFunction *of) |
| virtual MESQUITE_EXPORT | ~FeasibleNewton () |
| MESQUITE_EXPORT void | set_lower_gradient_bound (double gradc) |
| PatchSet * | get_patch_set () |
| MESQUITE_EXPORT std::string | get_name () const |
| Get string name for use in diagnostic and status output. | |
Protected Member Functions | |
| virtual void | initialize (PatchData &pd, MsqError &err) |
| virtual void | optimize_vertex_positions (PatchData &pd, MsqError &err) |
| virtual void | initialize_mesh_iteration (PatchData &pd, MsqError &err) |
| virtual void | terminate_mesh_iteration (PatchData &pd, MsqError &err) |
| virtual void | cleanup () |
Private Attributes | |
| double | convTol |
| MsqHessian | mHessian |
| PatchDataVerticesMemento * | coordsMem |
| bool | havePrintedDirectionMessage |
High Performance implementation of the Feasible Newton algorythm.
Consider our non-linear objective function \( f: I\!\!R^{3N} \rightarrow I\!\!R \) where \( N \) is the number of vertices of the mesh, and \( 3N \) is therefore the number of degrees of freedom of the mesh. The Taylor expansion of \( f \) around the point \( x_0 \) is
\[ f(x_0+d) = f(x_0) + \nabla f(x_0)d + \frac{1}{2} d^T\nabla^2 f(x_0)d + ... \;\;\; .\]
Each iteration of the Newton algorithm tries to find a descent vector that minimizes the above quadratic approximation, i.e. it looks for
\[ \min_{d} q(d;x_0) = f(x_0) + \nabla f(x_0)d + \frac{1}{2} d^T\nabla^2 f(x_0)d \;\; . \]
We know that if a quadratic function has a finite minimum, it is reached at the point where the function gradient is null and that the function Hessian is then positive definite. Therefore we are looking for \( d \) such that \( \nabla q(d;x_0) =0 \). We have
\[ \nabla q(d;x_0) = \nabla f(x_0) + \nabla^2 f(x_0)d \;\;, \]
therefore we must solve for \( d \) the system
\[ \nabla^2 f(x_0)d = -\nabla f(x_0) \;\; . \]
We assume that the Hessian is positive definite and we use the conjugate gradient algebraic solver to solve the above system. If the conjugate gradient solver finds a direction of negative curvature, the Hessian was not positive definite and we take a step in that direction of negative curvature, which is a descent direction.
Definition at line 95 of file FeasibleNewton.hpp.
Definition at line 67 of file FeasibleNewton.cpp.
References MBMesquite::TerminationCriterion::add_absolute_gradient_L2_norm(), and MBMesquite::QualityImprover::get_inner_termination_criterion().
: VertexMover( of ), PatchSetUser( true ), convTol( 1e-6 ), coordsMem( 0 ), havePrintedDirectionMessage( false ) { TerminationCriterion* default_crit = get_inner_termination_criterion(); default_crit->add_absolute_gradient_L2_norm( 5e-5 ); }
| virtual MESQUITE_EXPORT MBMesquite::FeasibleNewton::~FeasibleNewton | ( | ) | [inline, virtual] |
| void FeasibleNewton::cleanup | ( | ) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 430 of file FeasibleNewton.cpp.
References coordsMem.
| std::string FeasibleNewton::get_name | ( | ) | const [virtual] |
Get string name for use in diagnostic and status output.
Implements MBMesquite::Instruction.
Definition at line 57 of file FeasibleNewton.cpp.
{
return "FeasibleNewton";
}
| PatchSet * FeasibleNewton::get_patch_set | ( | ) | [virtual] |
Reimplemented from MBMesquite::PatchSetUser.
Definition at line 62 of file FeasibleNewton.cpp.
{
return PatchSetUser::get_patch_set();
}
| void FeasibleNewton::initialize | ( | PatchData & | pd, |
| MsqError & | err | ||
| ) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 74 of file FeasibleNewton.cpp.
References coordsMem, MBMesquite::PatchData::create_vertices_memento(), havePrintedDirectionMessage, and MSQ_CHKERR.
{
// Cannot do anything. Variable sizes with maximum size dependent
// upon the entire MeshSet.
coordsMem = pd.create_vertices_memento( err );MSQ_CHKERR( err );
havePrintedDirectionMessage = false;
}
| void FeasibleNewton::initialize_mesh_iteration | ( | PatchData & | pd, |
| MsqError & | err | ||
| ) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 82 of file FeasibleNewton.cpp.
References MBMesquite::PatchData::reorder().
{
pd.reorder();
}
| void FeasibleNewton::optimize_vertex_positions | ( | PatchData & | pd, |
| MsqError & | err | ||
| ) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 87 of file FeasibleNewton.cpp.
References MBMesquite::TerminationCriterion::accumulate_inner(), MBMesquite::TerminationCriterion::accumulate_patch(), MBMesquite::arrptr(), MBMesquite::MsqError::BARRIER_VIOLATED, beta, MBMesquite::MsqHessian::cg_solver(), MBMesquite::MsqError::clear(), convTol, coordsMem, MBMesquite::PatchData::domain_set(), epsilon, MBMesquite::MsqError::error_code(), MBMesquite::OFEvaluator::evaluate(), MBMesquite::PatchData::get_domain(), MBMesquite::QualityImprover::get_inner_termination_criterion(), MBMesquite::VertexMover::get_objective_function_evaluator(), MBMesquite::grad(), havePrintedDirectionMessage, MBMesquite::MsqHessian::initialize(), MBMesquite::inner(), INTERNAL_ERROR, MBMesquite::length(), mHessian, MBMesquite::PatchData::move_free_vertices_constrained(), MSQ_DBG, MSQ_DBGOUT, MSQ_ERRRTN, MSQ_FUNCTION_TIMER, MSQ_PRINT, MSQ_SETERR, MBMesquite::MsqHessian::norm(), MBMesquite::PatchData::num_free_vertices(), MBMesquite::PatchData::recreate_vertices_memento(), MBMesquite::PatchData::set_to_vertices_memento(), MBMesquite::TerminationCriterion::terminate(), and MBMesquite::OFEvaluator::update().
{
MSQ_FUNCTION_TIMER( "FeasibleNewton::optimize_vertex_positions" );
MSQ_DBGOUT( 2 ) << "\no Performing Feasible Newton optimization.\n";
//
// the only valid 2D meshes that FeasibleNewton works for are truly planar which
// lie in the X-Y coordinate plane.
//
XYPlanarDomain* xyPlanarDomainPtr = dynamic_cast< XYPlanarDomain* >( pd.get_domain() );
// only optimize if input mesh is a volume or an XYPlanarDomain
if( !pd.domain_set() || xyPlanarDomainPtr != NULL )
{
const double sigma = 1e-4;
const double beta0 = 0.25;
const double beta1 = 0.80;
const double tol1 = 1e-8;
const double tol2 = 1e-12;
const double epsilon = 1e-10;
double original_value, new_value;
double beta;
int nv = pd.num_free_vertices();
std::vector< Vector3D > grad( nv ), d( nv );
bool fn_bool = true; // bool used for determining validity of patch
OFEvaluator& objFunc = get_objective_function_evaluator();
int i;
// TODD -- Don't blame the code for bad things happening when using a
// bad termination test or requesting more accuracy than is
// possible.
//
// Also,
// 1. Allocate a hessian and calculate the sparsity pattern.
mHessian.initialize( pd, err );MSQ_ERRRTN( err );
// does the Feasible Newton iteration until stopping is required.
// Terminate when inner termination criterion signals.
/* Computes the value of the stopping criterion*/
TerminationCriterion* term_crit = get_inner_termination_criterion();
while( !term_crit->terminate() )
{
fn_bool = objFunc.update( pd, original_value, grad, mHessian, err );MSQ_ERRRTN( err );
if( !fn_bool )
{
MSQ_SETERR( err )
( "invalid patch for hessian calculation", MsqError::INTERNAL_ERROR );
return;
}
if( MSQ_DBG( 3 ) )
{ // avoid expensive norm calculations if debug flag is off
MSQ_DBGOUT( 3 ) << " o objective function: " << original_value << std::endl;
MSQ_DBGOUT( 3 ) << " o gradient norm: " << length( grad ) << std::endl;
MSQ_DBGOUT( 3 ) << " o Hessian norm: " << mHessian.norm() << std::endl;
}
// Prints out free vertices coordinates.
//
// Comment out the following because it is way to verbose for larger
// problems. Consider doing:
// inner_term_crit->write_mesh_steps( "filename.vtk" );
// instead.
// - j.kraftcheck 2010-11-17
// if (MSQ_DBG(3)) {
// MSQ_DBGOUT(3) << "\n o Free vertices ("<< pd.num_free_vertices()
// <<")original coordinates:\n ";
// MSQ_ERRRTN(err);
// const MsqVertex* toto1 = pd.get_vertex_array(err); MSQ_ERRRTN(err);
// MsqFreeVertexIndexIterator ind1(pd, err); MSQ_ERRRTN(err);
// ind1.reset();
// while (ind1.next()) {
// MSQ_DBGOUT(3) << "\t\t\t" << toto1[ind1.value()];
// }
// }
// 4. Calculate a direction using preconditionned conjugate gradients
// to find a zero of the Newton system of equations (H*d = -g)
// (a) stop if conjugate iteration limit reached
// (b) stop if relative residual is small
// (c) stop if direction of negative curvature is obtained
mHessian.cg_solver( arrptr( d ), arrptr( grad ), err );MSQ_ERRRTN( err );
// 5. Check for descent direction (inner produce of gradient and
// direction is negative.
double alpha = inner( grad, d );
// TODD -- Add back in if you encounter problems -- do a gradient
// step if the direction from the conjugate gradient solver
// is not a descent direction for the objective function. We
// SHOULD always get a descent direction from the conjugate
// method though, unless the preconditioner is not positive
// definite.
// If direction is positive, does a gradient (steepest descent) step.
if( alpha > -epsilon )
{
MSQ_DBGOUT( 3 ) << " o alpha = " << alpha << " (rejected)" << std::endl;
if( !havePrintedDirectionMessage )
{
MSQ_PRINT( 1 )
( "Newton direction not guaranteed descent. Ensure preconditioner is positive "
"definite.\n" );
havePrintedDirectionMessage = true;
}
// TODD: removed performing gradient step here since we will use
// gradient if step does not produce descent. Instead we set
// alpha to a small negative value.
alpha = -epsilon;
// alpha = inner(grad, grad, nv); // compute norm squared of gradient
// if (alpha < 1) alpha = 1; // take max with constant
// for (i = 0; i < nv; ++i) {
// d[i] = -grad[i] / alpha; // compute scaled gradient
// }
// alpha = inner(grad, d, nv); // recompute alpha
// // equal to one for large gradient
}
else
{
MSQ_DBGOUT( 3 ) << " o alpha = " << alpha << std::endl;
}
alpha *= sigma;
beta = 1.0;
pd.recreate_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
// TODD: Unrolling the linesearch loop. We do a function and
// gradient evaluation when beta = 1. Otherwise, we end up
// in the linesearch regime. We expect that several
// evaluations will be made, so we only do a function evaluation
// and finish with a gradient evaluation. When beta = 1, we also
// check the gradient for stability.
// TODD -- the Armijo linesearch is based on the objective function,
// so theoretically we only need to evaluate the objective
// function. However, near a very accurate solution, say with
// the two norm of the gradient of the objective function less
// than 1e-5, the numerical error in the objective function
// calculation is enough that the Armijo linesearch will
// fail. To correct this situation, the iterate is accepted
// when the norm of the gradient is also small. If you need
// high accuracy and have a large mesh, talk with Todd about
// the numerical issues so that we can fix it.
// TODD -- the Armijo linesearch here is only good for differentiable
// functions. When projections are involved, you should change
// to a form of the linesearch meant for nondifferentiable
// functions.
pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
fn_bool = objFunc.evaluate( pd, new_value, grad, err );
if( err.error_code() == err.BARRIER_VIOLATED )
err.clear(); // barrier violated does not represent an actual error here
MSQ_ERRRTN( err );
if( ( fn_bool && ( original_value - new_value >= -alpha * beta - epsilon ) ) ||
( fn_bool && ( length( arrptr( grad ), nv ) < 100 * convTol ) ) )
{
// Armijo linesearch rules passed.
MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (accepted without line search)" << std::endl;
}
else
{
if( !fn_bool )
{
// Function undefined. Use the higher decrease rate.
beta *= beta0;
MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (invalid step)" << std::endl;
}
else
{
// Function defined, but not sufficient decrease
// Use the lower decrease rate.
beta *= beta1;
MSQ_DBGOUT( 3 ) << " o beta = " << beta << " (insufficient decrease)" << std::endl;
}
pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
// Standard Armijo linesearch rules
MSQ_DBGOUT( 3 ) << " o Doing line search" << std::endl;
while( beta >= tol1 )
{
// 6. Search along the direction
// (a) trial = x + beta*d
pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
// (b) function evaluation
fn_bool = objFunc.evaluate( pd, new_value, err );
if( err.error_code() == err.BARRIER_VIOLATED )
err.clear(); // barrier violated does not represent an actual error here
MSQ_ERRRTN( err );
// (c) check for sufficient decrease and stop
if( !fn_bool )
{
// function not defined at trial point
beta *= beta0;
}
else if( original_value - new_value >= -alpha * beta - epsilon )
{
// iterate is acceptable.
break;
}
else
{
// iterate is not acceptable -- shrink beta
beta *= beta1;
}
pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
}
if( beta < tol1 )
{
// assert(pd.set_to_vertices_memento called last)
// TODD -- Lower limit on steplength reached. Direction does not
// appear to make sufficient progress decreasing the
// objective function. This can happen when you are
// very near a solution due to numerical errors in
// computing the objective function. It can also happen
// when the direction is not a descent direction and when
// you are projecting the iterates onto a surface.
//
// The latter cases require the use of a linesearch on
// a gradient step. If this linesearch terminate with
// insufficient decrease, then you are at a critical
// point and should stop!
//
// The numerical errors with the objective function cannot
// be overcome. Eventually, the gradient step will
// fail to compute a new direction and you will stop.
MSQ_PRINT( 1 )
( "Sufficient decrease not obtained in linesearch; switching to gradient.\n" );
alpha = inner( arrptr( grad ), arrptr( grad ),
nv ); // compute norm squared of gradient
if( alpha < 1 ) alpha = 1; // take max with constant
for( i = 0; i < nv; ++i )
{
d[i] = -grad[i] / alpha; // compute scaled gradient
}
alpha = inner( arrptr( grad ), arrptr( d ), nv ); // recompute alpha
alpha *= sigma; // equal to one for large gradient
beta = 1.0;
// Standard Armijo linesearch rules
while( beta >= tol2 )
{
// 6. Search along the direction
// (a) trial = x + beta*d
pd.move_free_vertices_constrained( arrptr( d ), nv, beta, err );MSQ_ERRRTN( err );
// (b) function evaluation
fn_bool = objFunc.evaluate( pd, new_value, err );
if( err.error_code() == err.BARRIER_VIOLATED )
err.clear(); // barrier violated does not represent an actual error
// here
MSQ_ERRRTN( err );
// (c) check for sufficient decrease and stop
if( !fn_bool )
{
// function not defined at trial point
beta *= beta0;
}
else if( original_value - new_value >= -alpha * beta - epsilon )
{
// iterate is acceptable.
break;
}
else
{
// iterate is not acceptable -- shrink beta
beta *= beta1;
}
pd.set_to_vertices_memento( coordsMem, err );MSQ_ERRRTN( err );
}
if( beta < tol2 )
{
// assert(pd.set_to_vertices_memento called last)
// TODD -- Lower limit on steplength reached. Gradient does not
// appear to make sufficient progress decreasing the
// objective function. This can happen when you are
// very near a solution due to numerical errors in
// computing the objective function. Most likely you
// are at a critical point for the problem.
MSQ_PRINT( 1 )
( "Sufficient decrease not obtained with gradient; critical point likely "
"found.\n" );
break;
}
}
// Compute the gradient at the new point -- needed by termination check
fn_bool = objFunc.update( pd, new_value, grad, err );MSQ_ERRRTN( err );
}
// Prints out free vertices coordinates.
// if (MSQ_DBG(3)) {
// MSQ_DBGOUT(3) << " o Free vertices new coordinates: \n";
// const MsqVertex* toto1 = pd.get_vertex_array(err); MSQ_ERRRTN(err);
// MsqFreeVertexIndexIterator ind(pd, err); MSQ_ERRRTN(err);
// ind.reset();
// while (ind.next()) {
// MSQ_DBGOUT(3) << "\t\t\t" << toto1[ind.value()];
// }
// }
// checks stopping criterion
term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
term_crit->accumulate_inner( pd, new_value, arrptr( grad ), err );MSQ_ERRRTN( err );
}
MSQ_PRINT( 2 )( "FINISHED\n" );
}
else
{
std::cout << "WARNING: Feasible Newton optimization only supported for volume meshes" << std::endl
<< " and XYPlanarDomain surface meshes." << std::endl
<< std::endl
<< "Try a different solver such as Steepest Descent." << std::endl;
}
}
| MESQUITE_EXPORT void MBMesquite::FeasibleNewton::set_lower_gradient_bound | ( | double | gradc | ) | [inline] |
Sets a minimum value for the gradient. If the gradient is below that value, we stop iterating.
Definition at line 107 of file FeasibleNewton.hpp.
References convTol.
{
convTol = gradc;
}
| void FeasibleNewton::terminate_mesh_iteration | ( | PatchData & | pd, |
| MsqError & | err | ||
| ) | [protected, virtual] |
Implements MBMesquite::VertexMover.
Definition at line 423 of file FeasibleNewton.cpp.
{
// Michael:: Should the vertices memento be delete here???
// cout << "- Executing FeasibleNewton::iteration_complete()\n";
}
double MBMesquite::FeasibleNewton::convTol [private] |
Definition at line 124 of file FeasibleNewton.hpp.
Referenced by optimize_vertex_positions(), and set_lower_gradient_bound().
Definition at line 126 of file FeasibleNewton.hpp.
Referenced by cleanup(), initialize(), optimize_vertex_positions(), and ~FeasibleNewton().
bool MBMesquite::FeasibleNewton::havePrintedDirectionMessage [private] |
Definition at line 127 of file FeasibleNewton.hpp.
Referenced by initialize(), and optimize_vertex_positions().
Definition at line 125 of file FeasibleNewton.hpp.
Referenced by optimize_vertex_positions().
1.7.6.1