MOAB: Mesh Oriented datABase  (version 5.3.0)
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 ()
MESQUITE_EXPORT void set_lower_gradient_bound (double gradc)
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

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

Reimplemented from MBMesquite::PatchSetUser.

Definition at line 62 of file FeasibleNewton.cpp.

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

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

Member Data Documentation

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines