Одномерный поиск корня с минимальным количеством вычислений функции

Я ищу алгоритмы поиска корней, которые используют очень мало вычислений функций (цель - минимум). Задача поиска корней имеет следующие характеристики:

f(x) = 0, R -> R

  • оценка функции (f(.)) чрезвычайно затратна*;
  • для начала доступен ограничивающий интервал ([a,b]) (относительно хорошее приближение, а не дикое предположение);
  • f(.) непрерывен;
  • f(.) дифференцируема (аналитическая производная недоступна);
  • известно, что только один корень лежит в пределах начального интервала ([a,b]);
  • плавно меняющееся f(.) (от функции не ожидается ничего экстремального);
  • допустимые критерии остановки, т.е. |f(x)| < 1e-2 достаточно.

* Мы можем с уверенностью предположить, что любое вычисление, выполненное алгоритмом, ничтожно мало по сравнению с единственной оценкой f(.). Таким образом, экономия даже одной оценки функции является значительным преимуществом.

Учитывая это, какие алгоритмы наиболее эффективны для поиска корня с наименьшим количеством вычислений функции?

На основе fzero Matlab и root-finding functions, метод Брента кажется наиболее популярным, хотя может существовать более эффективный алгоритм для конкретной задачи. описано выше.

Также приветствуются ссылки на книги и обзорные статьи.


python matlab numerical-methods
person rozsasarpi    schedule 17.11.2016    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Что ж, если вы хотите попробовать метод Брента, вы можете попробовать мой перевод оригинального алгоритма ниже. Это всего лишь мой перевод исходного кода C в MATLAB. Я убедился, что все тестовые примеры для исходного кода C дают идентичные результаты в моем переводе.

В приведенном ниже коде This — это дескриптор функции, а a и b — границы поиска.

function x = brent(This,a,b) 

tol = 1e-2 ;
maxit = 1e+3 ;
displayIter = true ;

Fa = This(a) ;
Fb = This(b) ;
c = a ;
Fc = Fa ;

it = 0 ;
if displayIter
    fprintf(1,'Searching for a root in the interval [%g,%g] using Brent''s method...\n',a,b) ;
end
while it<maxit
    it = it + 1 ;

    prevStep = b - a ;

    if abs(Fc) < abs(Fb)
        % swap for b to be best approximation
        a = b ;
        b = c ;
        c = a ;
        Fa = Fb ;
        Fb = Fc ;
        Fc = Fa ;
    end

    tolAct =  2*eps*abs(b) + tol/2 ;

    newStep = (c-b)/2 ;

    if abs(newStep) <= tolAct || abs(Fb)<eps
        % acceptable solution found
        x = b ;
        return 
    end

    if abs(prevStep) >= tolAct && Fa == Fb
        % previous step was large enough and in the right direction, try
        % interpolation
        cb = c-b ;
        if abs(a-c)<eps
            % if only two distinct points, interpolate linearly
            t1 = Fb/Fa ;
            p = cb*t1 ;
            q = 1 - t1 ;
        else
            % three points, do quadratic inverse  interpolation
            a = Fa/Fc ;
            t1 = Fb/Fc ;
            t2 = Fb/Fa ;
            p = t2*( cb*q*(q-t1) - (b-a)*(t1-1) ) ;
            q = (q-1)*(t1-1)*(t2-1) ;
        end
        if p>0
            q = -q ;
        else
            p = -p ;
        end
        if p < ( 0.75*cb*q-abs(tolAct*q)/2 ) && p < abs(prevStep*q/2)
            newStep = p/q ;
        end
    end
    % step must be at least as large as tolerance
    if abs(newStep)<tolAct
        if newStep>0
            newStep = tolAct ;
        else
            newStep = -tolAct ;
        end
    end

    a = b ;
    Fa = Fb ;
    b = b + newStep ;
    Fb = This(b) ;
    if ( Fb > 0 && Fc > 0 ) || ( Fb < 0 && Fc < 0 )
        c = a ;
        Fc = Fa ;
    end
    if displayIter
        fprintf(1,'[%g] Gap = %g, dx = %g\n',it,abs(Fb),newStep) ;
    end
end
fprintf(1,'[%g] Gap = %g, dx = %g\n',it,abs(Fb),newStep) ;

end
person transversality condition    schedule 05.12.2016
comment
Спасибо за ваш код (3). Я уже протестировал метод Брента, кроме того, что fzero (1) я использую эта(2) реализация Matlab. Быстрый тест с использованием fun = @(x) exp(-exp(-x)) - x и x0 = [0, 1] с tolx = 1e-2 показывает следующие номера оценки функции: (1) 5; (2) 4; (3) 7. Я ищу что-то более эффективное. :) - person rozsasarpi; 05.12.2016