01-05-11 - Golden 1d Searches

GoldenSearch1d finds the minimum of some function if you know the finite range to look in. Search1d_ExpandingThenGoldenDown looks in the infinite interval >= 0 . Pretty trivial thing, but handy.

The Golden ratio arises from doing a four point search and trying to reuse evaluations. If your four points are { 0, rho, 1-rho ,1 } and you shrink to the lower three { 0, rho, 1-rho }, then you want to reuse the rho evaluation as your new high interior point, so you require (1-rho)^2 = rho , hence 1-rho = (sqrt(5)-1)/2

BTW because you have various points you could obviously use interpolation search (either linear or quadratic) at various places here and reduce the number of evaluations needed to converge. See Brent's method or Dekker's method .

Also obviously this kind of stuff only works on functions with simple minima.

Also obviously if func is analytic and you can take derivatives of it there are better ways. I use this for evaluating things that aren't simple functions, but have nice shapes (such as running my image compressor with different quantizers).

---

template< typename t_functor >  
double GoldenSearch1d( t_functor func, double lo, double v_lo, double hi, double v_hi, double minstep )
{
    const double rho = 0.381966;
    const double irho = 1.0 - rho; // = (sqrt(5)-1)/2 

    // four points :
    // [lo,m1,m2,hi]

    double m1 = irho*lo + rho*hi;
    double m2 = irho*hi + rho*lo;

    double v_m1 = func( m1 );
    double v_m2 = func( m2 );
    
    while( (m1-lo) > minstep )
    {
        // step to [lo,m1,m2] or [m1,m2,hi]
        // only one func eval per iteration :
        if ( MIN(v_lo,v_m1) < MIN(v_hi,v_m2) )
        {
            hi = m2; v_hi = v_m2;
            m2 = m1; v_m2 = v_m1;
            m1 = irho*lo + rho*hi;
            v_m1 = func( m1 );
        }
        else
        {
            lo = m1; v_lo = v_m1;
            m1 = m2; v_m1 = v_m2;
            m2 = irho*hi + rho*lo;
            v_m2 = func( m2 );
        }
        
        ASSERT( fequal(m2, irho*hi + rho*lo) );
        ASSERT( fequal(m1, irho*lo + rho*hi) );
    }
    
    // could do a cubic fit with the 4 samples I have now
    //  but they're close together so would be numerically unstable
    
    /*
    //return (lo+hi)/2.0;
    
    if ( v_m1 < v_m2 ) return m1;
    else return m2;
    
    /*/
    // return best of the 4 samples :
    if ( v_lo < v_m1 ) { v_m1 = v_lo; m1 = lo; }
    if ( v_hi < v_m2 ) { v_m2 = v_hi; m2 = hi; }
    
    if ( v_m1 < v_m2 ) return m1;
    else return m2;
    /**/
}

template< typename t_functor >  
double GoldenSearch1d( t_functor func, double lo, double hi, double minstep )
{
    double v_lo = func( lo );
    double v_hi = func( hi );
    
    return GoldenSearch1d( func, lo, v_lo, hi, v_hi, minstep );
}

template< typename t_functor >  
double Search1d_ExpandingThenGoldenDown( t_functor func, double start, double step, double minstep , const int min_steps = 8)
{
    struct Triple
    {
        double t0,f0;
        double t1,f1;
        double t2,f2;
    };
    
    Triple cur;
    cur.t2 = cur.t1 = cur.t0 = start;
    cur.f2 = cur.f1 = cur.f0 = func( start );
        
    int steps = 0;
    
    Triple best = cur;
    
    for(;;)
    {
        cur.t0 = cur.t1; 
        cur.f0 = cur.f1;
        cur.t1 = cur.t2; 
        cur.f1 = cur.f2;
        
        cur.t2 = cur.t1 + step;
        cur.f2 = func( cur.t2 );
        
        if ( cur.f1 <= best.f1 )
        {
            best = cur;
        }
        
        // if we got worst and we're past min_steps :
        if ( cur.f2 > cur.f1 && steps > min_steps )
            break;
                
        const double golden_growth = 1.618034; // 1/rho - 1
        step *= golden_growth; // grow step by some amount
        steps++;
    }
        
    // best is at t1 bracketed in [t0,t2]   
    // could save one function eval by passing in t1,f1 as well
    
    return GoldenSearch1d(func,best.t0,best.f0,best.t2,best.f2,minstep);
}

---

Usage example :

---

#define MAKE_FUNCTOR(type,func) \
struct STRING_JOIN(func,_functor) { \
  type operator() (type x) { return func(x); } \
};

double TFunc( double x )
{
    return 100 / ( x + 1) + x;
}

MAKE_FUNCTOR(double,TFunc);

int main(int argc,char *argv[])
{

double t = Search1d_ExpandingThenGoldenDown(TFunc_functor(),0.0,1.0,0.0001);

lprintfvar(t);
lprintf(TFunc(t) , "\n");

return 0;
}