optimise.h

Go to the documentation of this file.

#pragma once
#include "dg/algorithm.h"
#include "tridiaginv.h" // lapack wrapper
 
namespace dg{
namespace mat{
template<class Gradient, class InvHessian, class ContainerType>
unsigned newton( Gradient grad, InvHessian invhess,
    ContainerType& x0, double tol = 1e-5, unsigned max_iter = 1000)
{
    ContainerType jj(x0), p(x0), test(x0);
    for ( unsigned i=0; i<max_iter; i++)
    {
        grad( x0, jj);
        invhess( x0, jj, p);
        double alpha = 1.;
        dg::blas1::axpby( -alpha, p, 1., x0);
        double err = sqrt(dg::blas1::dot( jj, jj));
        if (err < tol)
            return i;
    }
    return max_iter;
}
 
template<class Func, class Jacobian, class ContainerType0, class ContainerType1>
unsigned levenberg_marquardt( Func fun, Jacobian jac,
    ContainerType0& x0,
    const ContainerType1& copyable, // size of return of fun
    double tol = 1e-8, unsigned max_iter = 1000)
{
    unsigned num_p = x0.size();
    auto x1 = x0, W = x0;
    auto rs(copyable), rs1(rs);
    std::vector<ContainerType1> jacs(num_p, copyable);
    dg::SquareMatrix<double> HH(num_p, 0.), evHH( HH), evHH_T(HH), WW(HH), syWW(WW);
    thrust::host_vector<double> evs( num_p), grad(num_p), gradbar(num_p),
        pk(num_p), pkbar(num_p);
    thrust::host_vector<double> work( 3*num_p-1);
    // init loop
    fun( x0, rs);
    jac( x0, jacs);
    double delta = 0;
    for( unsigned p=0; p<num_p; p++)
    {
        WW(p,p) = W[p] = dg::blas1::dot( jacs[p], jacs[p]);
    }
    for( unsigned p=0; p<num_p; p++)
    {
        grad[p] = dg::blas1::dot( jacs[p], rs);
        delta += grad[p]*grad[p]/W[p];
    }
    //std::cout << "Norm gT g/W^2 " << sqrt(delta)<<"\n";
    double f0 = dg::blas1::dot( rs, rs);
    delta = 0.25*f0/sqrt(delta);
    //std::cout << "initial delta " << delta<<"\n";
    double normx0 = sqrt(dg::blas1::dot( x0, x0));
    for ( unsigned k=0; k<max_iter; k++)
    {
        // 1. Solve (J^T J + lambda W)p = -J^T r with lambda : ||p||_W leq Delta
        // W[p] = max{ W_{k-1}, H_pp}
        // In: jacs, rs; Out: pk, normpk, lambda
        for( unsigned l=0; l<num_p; l++)
        {
            for( unsigned j=l; j<num_p; j++)
                HH(j,l) = HH(l,j) = dg::blas1::dot( jacs[l], jacs[j]);
            WW(l,l) = std::max( WW(l,l), HH(l,l)); // After Mor´e the weights must not decrease
        }
        // !!! sygv destroys WW on output (so use copy syWW instead) !!!
        syWW = WW;
        lapack::sygv( 1, 'V', 'U', num_p, HH.data(), num_p, syWW.data(), num_p, evs, work);
        evHH_T = HH;
        evHH = HH.transpose();
        //std::cout << "#########Iteration "<<k<<"\n";
        //std::cout << "Eigenvalues are \n";
        //for( unsigned p=0; p<num_p; p++)
        //    std::cout << "p EV "<<evs[p]<<"\n";
        //std::cout << "Weights are \n";
        //for( unsigned p=0; p<num_p; p++)
        //    std::cout << "W "<<WW(p,p)<<"\n";
        dg::blas2::gemv( evHH_T, grad, gradbar); // !! sygv gives A = WE_A Lambda E_A^TW
        double normp=0;
        auto target = [&]( double lambda)
        {
            // safeguard against 0 Eigenvalue!
            for( unsigned p=0; p<num_p; p++)
                pkbar[p]= -gradbar[p]/(evs[p] +lambda == 0 ? 1e-16 : evs[p]+lambda);
            normp = sqrt(dg::blas1::dot( pkbar, pkbar));
            return 1./delta - 1./normp;
        };
        auto dtarget = [&]( double lambda)
        {
            // safeguard against 0 Eigenvalue!
            double dnorm =0;
            for( unsigned p=0; p<num_p; p++)
                dnorm += pkbar[p]*pkbar[p]/(evs[p]+lambda == 0 ? 1e-16 : evs[p]+lambda);
            return - dnorm/normp/normp/normp;
        };
        double lambda = 0;
        const double sigma = 1e-4; // tolerance for Newton algorithm
        const unsigned max_newton = 100;
        double phi = target(lambda);
        if( phi > 0) // positive lambda only exist if phi(0) > 0 since dlambda < 0 for all lambda >= 0
        {
            for( unsigned i=0; i<max_newton; i++) // Safeguard
            {
                // if ||p|| leq delta ( 1+sigma)
                if ( fabs(phi) <= sigma/normp || i == max_newton-1)
                    break;
                // The first step is always to the right ... (phi > 0 , lambda < 0)
                lambda += - phi/dtarget(lambda);
                phi = target(lambda);
            }
        }
        dg::blas2::gemv(evHH, pkbar, pk);
        // target(lambda) updates grad and normp
        // 2. Check termination
        //std::cout << "Real Norm p "<<sqrt(dg::blas1::dot( pk, pk))<<" normx0 "<<normx0<<"\n";
        if( sqrt(dg::blas1::dot( pk,pk)) <= tol*(normx0 + 1.))
        {
            dg::blas1::axpby( 1., pk , 1., x0);
            return k;
        }
        // 3. Compute Ratio rhok
        // don't overwrite (x0, rs) because we might reject
        dg::blas1::axpby( 1., pk , 1., x0, x1);
        //for( unsigned l=0; l<num_p; l++)
        //    std::cout << "x0 "<<l<<" "<<x0[l]<<"\n";
        //for( unsigned l=0; l<num_p; l++)
        //    std::cout << "x1 "<<l<<" "<<x1[l]<<"\n";
        double f1 = 0;
        try{
            fun( x1, rs1);
            f1 = dg::blas1::dot( rs1, rs1);
        }
        catch( dg::Error& e)
        {
            // The only reason to throw an Error here is if there is a NaN or Inf
            // In this case will rs1 with huge values such the step is rejected
            dg::blas1::copy( 1e300, rs1);
            f1 = 1e300;
            //std::cout << "############################# STEP CONTAINS NAN OR INF! REJECT!\n";
        }
        // compute (Jp)^2
        double jp = dg::blas2::dot( pkbar, evs, pkbar);
        //std::cout<< "f0 "<<f0<<" f1 "<<f1<<"\n";
        double rhok = (1.-f1/f0)/( jp/f0 + 2*lambda*normp*normp/f0);
        //std::cout << "Actual    reduction "<<f0-f1<<"\n";
        //std::cout << "Predicted reduction "<<( jp + 2*lambda*normp*normp)<<"\n";
        //std::cout << "Ratio "<<rhok<<"\n";
        // Algorithm 4.1 from Nocedal & Wright
        if( rhok < 0.25)
            delta = 0.25*delta;
        else
        {
            // Mor´e has slightly different conditions
            if( rhok > 0.75 && lambda > 0) // steps lies on the trust region boundary
                delta = 2*delta;
            // else delta remains unchanged
        }
        const double eta = 1e-4; // when step is rejected
        if ( rhok > eta)
        {
            x0 = x1;
            // update all quantities
            f0 = f1;
            using std::swap;
            swap( rs, rs1);
            jac( x0, jacs);
            normx0 = sqrt(dg::blas1::dot( x0, x0));
            for( unsigned l=0; l<num_p; l++)
            {
                grad[l] = dg::blas1::dot( jacs[l], rs);
                //std::cout << "grad "<<l<<" "<<grad[l]<<"\n";
            }
        }
        else
        {
        //    std::cout << "REJECTED\n";
        // else step is rejected
        }
    }
    return max_iter;
}
} //namespace mat
} //namespace dg