MOAB: Mesh Oriented datABase  (version 5.4.1)
MBMesquite::FeasibleNewton Class Reference

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 ()
PatchSetget_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
PatchDataVerticesMementocoordsMem
bool havePrintedDirectionMessage

## Detailed Description

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.

## Constructor & Destructor Documentation

 FeasibleNewton::FeasibleNewton ( ObjectiveFunction * of )

Definition at line 67 of file FeasibleNewton.cpp.

    : VertexMover( of ), PatchSetUser( true ), convTol( 1e-6 ), coordsMem( 0 ), havePrintedDirectionMessage( false )
{
TerminationCriterion* default_crit = get_inner_termination_criterion();
}

 virtual MESQUITE_EXPORT MBMesquite::FeasibleNewton::~FeasibleNewton ( )  [inline, virtual]

Definition at line 100 of file FeasibleNewton.hpp.

References coordsMem.

    {
delete coordsMem;
}


## Member Function Documentation

 void FeasibleNewton::cleanup ( )  [protected, virtual]

Implements MBMesquite::VertexMover.

Definition at line 430 of file FeasibleNewton.cpp.

References coordsMem.

{
delete coordsMem;
coordsMem = NULL;
}

 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.

{
// 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.

{
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" );
// - 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.

    {
}

 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";
}


## Member Data Documentation

 double MBMesquite::FeasibleNewton::convTol [private]

Definition at line 124 of file FeasibleNewton.hpp.

Referenced by optimize_vertex_positions(), and set_lower_gradient_bound().

 PatchDataVerticesMemento* MBMesquite::FeasibleNewton::coordsMem [private]
 bool MBMesquite::FeasibleNewton::havePrintedDirectionMessage [private]

Definition at line 127 of file FeasibleNewton.hpp.

Referenced by initialize(), and optimize_vertex_positions().

 MsqHessian MBMesquite::FeasibleNewton::mHessian [private]

Definition at line 125 of file FeasibleNewton.hpp.

Referenced by optimize_vertex_positions().

List of all members.

The documentation for this class was generated from the following files: