Текст подпрограммы и версий
mnk3r_c.zip 		Тексты тестовых примеров
tmnk3r_c.zip

Подпрограмма:  mnk3r_c

Назначение

Решение задачи минимизации диффеpенциpуемой функции многих переменных по заданному направлению (задача одномерной минимизации) методами линейной и квадратичной аппроксимации производной по направлению.

Математическое описание

Для решения задачи 

    min  f(x0 - h p0) ,   0 ≤ h ≤ hL < +∞ ,   x0 ∈ EN ,   p0 ∈ EN ,
      h

x0 и  p0 заданы, используются последовательно методы линейной и квадратичной аппроксимации производной по направлению  ( f ' (x0 - h p0), p0 ) .

Минимум  f (x) по направлению  p0 считается найденным, если выполнено одно из следующих условий:

1.  | hK - hK + 1| < EPS1 , где  hK, hK + 1 - два последовательных шага по направлению  p0,  EPS1 - заданная абсолютная точность одномерной минимизации;
2.  | hK - hK + 1| / hK < EPS2, где EPS2 - заданная относительная точность одномерной минимизации;
3.  | fi' (x0 - hK p0)| < EPSGi для всех  i = 1, 2, ..., N, где EPSG - вектоp покомпонентной точности вычисления градиента в точке  x0 - hK p0;
4.  
| ( fi' (x0 - hK p0), p0) | < max   { max  EPSGi, 10 -1 -EM/2} ,
                                                  i = 1, 2, ..., N
   где  EM - максимальный порядок абсолютной величины
   вещественной константы;
5.  K ≥ Kmax ,  где  Kmax - максимально допустимое количество вычислений градиента.

Алгоритм является модификацией алгоритма одномерной минимизации, изложенного в работе:

B.A.Cкоков, Стандартная программа минимизации диффеpенциpуемой функции многих переменных методом сопряженных градиентов, серия "Стандартные программы решения задач математического программирования", вып.10, Изд - во МГУ, 1968, ротапринт.

Использование 

    int mnk3r_c (integer *n, integer *i0, real *p, real *hl,
            real *h, real *x, real *f, real *fe, real *g, real *ge,
            real *func, real *grad, real *up, integer *ip, real *rm)

   Параметры
n - размерность пространства переменных (тип: целый);
i0 - целый вектоp длины  n, задающий фиксированные на время минимизации компоненты вектоpа  x: если  i0 (i) = 0, то  i - ая компонента вектоpа  x остается равной своему начальному значению;
p - вещественный вектоp длины  n, задающий направление одномерной минимизации  p0;
hl - правый конец отрезка одномерной минимизации  hl (тип: вещественный);
h - вещественная переменная, содержащая на входе заданный начальный шаг одномерного поиска, а на выходе - расстояние от начальной точки поиска до вычисленной точки минимума;
x - вещественный вектоp длины  n, содержащий на входе начальную точку  x0 поиска минимума  f (x) по направлению  p0, а на выходе - точку минимального вычисленного значения  f (x);
f - вещественная переменная, содержащая на входе значение функции в начальной точке  x0, а на выходе - вычисленное минимальное значение функции вдоль направления  p0;
fe - вещественная переменная, содержащая на входе заданную точность вычисления значений функции, а на выходе - точность вычисления функции в точке минимума;
g - вещественный вектоp длины  n, содержащий на входе градиент функции в начальной точке  x0, а на выходе - градиент функции в вычисленной точке минимума;
ge - вещественный вектоp длины  n, содержащий на входе заданную покомпонентную точность вычисления градиента, а на выходе - точность вычисления градиента в точке минимума;
func - имя подпрограммы вычисления значения функции  f (x) (см. замечания по использованию);
grad - имя подпрограммы вычисления градиента функции (см. замечания по использованию);
up - вещественный вектоp длины 6, задающий упpавляющие параметры алгоритма:
up(1) - заданный признак; если  up (1) = 2, то для построения линейной аппроксимации используется алгоритм удвоения пробных шагов; если  up (1) = 1, то используется постоянный шаг (см. замечания по использованию);
up(2) - заданная точность одномерной минимизации (относительная при  up (3) = 1, абсолютная при  up (3) ≠ 1);
up(3) - признак:  up (3) = 1, если задана относительная точность одномерной минимизации;  up (3) ≠ 1, если задана абсолютная точность;
up(4) - заданное максимально допустимое количество вычислений градиента;
up(5) - максимальный порядок абсолютной величины вещественной константы;
up(6) - максимальное количество десятичных значащих цифр в вещественной константе;
ip - целый вектоp длины 4, содержащий на выходе:
ip(1) - количество вычислений функции при одномерной минимизации;
ip(2) - количество вычислений градиента при одномерной минимизации;
ip(3) - признак: если  ip (3) = 0, то вычисление градиента  g в точке минимума не производилось; если  ip (3) = 1, то в  g содержится вычисленное значение градиента в точке минимума;
ip(4) - признак: если  ip (4) = 0, то одномерный минимум функции найден и лежит внутpи отрезка [0, hl]; если  ip (4) = 1, то минимум достигается в точке  hl; если  ip (4) = 2, то функция  f (x) не убывает по направлению  p0 и алгоритм не может обеспечить сходимость (см. замечания по использованию);
rm - вещественный вектоp длины 2 * n, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

Подпрограммы func и grad составляются пользователем.

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

           int func(float *x, float *f, float *fe)

        Параметры
        x   - вещественный вектор длины  n, задающий точку
                пространства, в которой вычисляется значение функции;
        f   - вещественная переменная, содержащая
                вычисленное значение функции в точке  x;
        fe - вещественная переменная, содержащая на входе
                заданную точность вычисления значения функции
                в точке  x; на выходе - достигнутую точность;

Если значение достигнутой точности вычисления функции не известно, то в теле подпрограммы func параметр fe не должен переопределяться.

Первый оператор подпрограммы вычисления градиента функции  f (x) должен иметь вид: 

           int grad(float *x, float *g, float *ge, int *i0)

        Параметры
        x   - вещественный вектор длины  n, задающий точку
                пространства, в которой вычисляется градиент;
        g   - вещественный вектоp длины  n, содержащий
                градиент функции в точке  x;
        ge - вещественный вектоp длины  n, содержащий  на
                входе заданную покомпонентную точность
                вычисления градиента функции; на выходе -
                достигнутую точность вычисления градиента;
        i0  - целый вектоp фиксированных компонент,
                упpавляющий вычислением компонент градиента:
                если  i0(i) = 0, то полагается  g(i) = 0.

Если значение достигнутой точности ge (i) для некоторого  i не известно, то в теле подпрограммы grad параметр ge не должен переопределяться.

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

Если решается задача безусловной минимизации функции по направлению  p0, то следует задать  hl = c, где  c - достаточно большое представимое в машине вещественное число.

B общем случае наибольшее влияние на эффективность программы оказывает выбор начального шага одномерной минимизации  h и значения параметра  up (1).

При одномерной минимизации функции  f (x) на направлении  p0 вначале ищется отрезок для построения линейной аппроксимации производной по направлению с помощью постоянного шага (при  up (1) = 1) или алгоритма удвоения пробных шагов (при  up (1) = 2). Затем последовательно осуществляется линейная и квадратичные аппроксимации производной по направлению.

Если функция многоэкстремальна, то при некоторых значениях  h алгоритм может не обеспечить сходимость к ближайшей к  x0 точке локального одномерного минимума. В этом случае, а также в случае неубывания функции по направлению  p0 (ip (4) = 2) следует уменьшить начальный шаг  h и задать  up (1) = 1.

Используются служебные подпрограммы: mnku1_c, mnku3_c, mnku5_c, mnku6_c, mnko2_c, mnko3_c.

Пример использования 

    min  f(x) ,  x ∈ E1 ,  x0 = 0.0 ,  p0 = -1.0 ,
    f(x)  =  100 e -x + x .
    Точка минимума  f(x) вдоль  p0 pавна   x*  =  - ln(0.01) ,
    f(x*)  =  1 - ln(0.01)

int main(void)
{
    /* Initialized data */
    static int n = 1;
    static int i0[1] = { 1 };
    static float p[1] = { -1.f };
    static float x[1] = { 0.f };
    static float h__ = .5f;
    static float hl = 1e10f;
    static float fe = 1e-7f;
    static float ge[1] = { 1e-7f };
    static float up[6] = { 2.f,1e-6f,2.f,50.f,18.f,12.f };

    /* Local variables */
    extern int grad_c(float *, float *, float *, int *),
            func_c(float *, float *, float *),
            mnk3r_c(int *, int *, float *, float *, float *, float *,
            float *, float *, float *, float *, S_fp, S_fp, float *,
            int *, float *);
    static float f, g[1];
    static int ip[4];
    static float rm[2];

    func_c(x, &f, &fe);
    grad_c(x, g, ge, i0);
    mnk3r_c(&n, i0, p, &hl, &h__, x, &f, &fe, g, ge, (S_fp)func_c,
            (S_fp)grad_c, up, ip, rm);

/* l1: */
    printf("\n %5i %5i \n", ip[0], ip[1]);
    printf("\n %16.7e \n", x[0]);
    printf("\n %16.7e \n", f);
    return 0;
} /* main */

int func_c(float *x, float *f, float *fe)
{
    /* Builtin functions */
    double exp(double);

    /* Parameter adjustments */
    --x;

    /* Function Body */
    *f = (float)exp((float)(-x[1])) * 100.f + x[1];
    return 0;
} /* func_c */

int grad_c(float *x, float *g, float *ge, int *i0)
{
    /* Builtin functions */
    double exp(double);

    /* Parameter adjustments */
    --i0;
    --ge;
    --g;
    --x;

    /* Function Body */
    g[1] = (float)exp((float)(-x[1])) * -100.f + 1.f;
    return 0;
} /* grad_c */


Результаты счета:

     вычислений функции  1        вычислений градиента  10

     x  =  4.6051700 + 00              f  =  5.6051700 + 00