Оптимизация оптимизации в MatLab: nested и anonymous functions

Добрый день!
Я занимаюсь научными исследованиями в области систем управления, и Matlab — мой основной рабочий инструмент. Одна из возможностей в MatLab — численная оптимизация. Оптимизировать (минимизировать) можно любую функцию, которая принимает на вход вектор варьируемых параметров и возвращает значение минимизируемого критерия. Естественно, в процессе оптимизации целевая функция вызывается множество раз и ее быстродействие существенно. В матлабе есть хорошие программные средства, которые часто позволяют существенно улучшить быстродействие, сохранив при этом читаемость и удобство сопровождения кода. Я приведу пример задачи, покажу на нём, что такое anonymous functions и nested functions, а потом покажу, как можно совместить эти два инструмента для заметного повышения быстродействия.


Постановка задачи и оптимизация

Я долго думал над примером для оптимизации, пытался выбрать что-то реалистичное – поиск оптимальной траектории для кинематически избыточного манипулятора, аппроксимация экспериментальных данных, идентификация параметров. Но все это как-то уводило от сути, так что решил остановится на непрактичном, но иллюстративном примере. Итак, заданы параметры a и b. Надо найти x в диапазоне [0;1], максимизирующий критерий:
image
где image — некоторая функция, не зависящая явно от x, но нужная для вычисления критерия. Пусть это будет модуль суммы первого миллиона псевдослучайных чисел, полученных при установка seed в z. Напишем матлабовскую функцию, вычисляющую значение критерия:
function J = GetJ(a,b,x) 
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
J=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end

Обратите внимание, что возвращаем значение со знаком минус, так как мы хотим найти максимум, а оптимизация ищет минимум.
Теперь для оптимизации при конкретных значениях a и b нам надо сделать анонимную функцию. Делается это так:
a=121233;
b=31235;
fhnd_GetJ = @(x) GetJ(a,b,x);

Теперь переменная fhnd_GetJ содержит хендл анонимной функции, которая принимает один параметр и позволяет вычислить GetJ(a,b, принятый параметр) для значений a и b, указанных при создании этой анонимной функции. Можно перейти непосредственно к минимизации. Скажем, мы хотим найти оптимальное значение x с точностью до 1 милионной:
opt=optimset('TolX',1e-6);
optimal_x=fminbnd(fhnd_GetJ,0,1,opt);

Функция fminbnd(fun,x_min,x_max) минимизирует скалярную функцию скалярного аргумента на интервале [x_min; x_max]. Здесь fun — хендл оптимизируемой функции. При выполнении оптимизации функция GetJ вызывается 12 раз, это можно увидеть, задав в опциях параметр ‘Display’ как ‘iter’ – покажет все итерации. На моем компе оптимизация занимает, в среднем, 10 мс.

Повышение быстродействия

Посмотрим, как это можно оптимизировать. Очевидно, что каждый раз при вызове функции GetJ мы совершенно зря считаем значения phi_a и phi_b, так как они при вариации x не меняются и зависят только от заданных значений a и b. Как тут можно сэкономить? Самый частый вариант, который мне встречается – вынести предварительные расчеты из целевой функции. Сделать вот такую функцию
function J = GetJ2(phi_a,phi_b,x)
J=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end

и вот такой код:
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
fhnd_GetJ2 = @(x) GetJ2(phi_a,phi_b,x);

optimal_x=fminbnd(fhnd_GetJ2,0,1,opt);

Конечно, это считается намного быстрее. Оптимизация проходит, в среднем, за 0.8 мс. Но платой за это является ухудшение кода. Разрывается целостность функционала – вычисление phi_a и phi_b самостоятельной ценности не имеет и в отрыве от вычисления J не нужно. Если где-то в другом месте опять потребуется считать J(a,b,x), уже не для оптимизации, а просто так, то придется вместо простого вызова GetJ за собой таскать еще и вычисление phi_a и phi_b или делать отдельно функции для оптимизации, отдельно для вычислений. Ну и просто не очень красиво.
Хорошо, что у матлаба есть способ обойти эту ситуацию. Для этого давайте создадим nested function внутри функции GetJ:
function J = GetJ3(a,b,x)
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));

J=nf_GetJ(x);

    function out_val = nf_GetJ(x)
        out_val=-(phi_a*x^2+phi_b*sqrt(1-x^2));
    end
end 

Nested function nf_GetJ видит все переменные в области видимости родительской функции и, в принципе, понятно, что и как она делает. Пока что мы никакого выигрыша по быстродействию не получили – все те же 10 мс на оптимизацию.

А теперь самое интересное – матлаб умеет работать с anonymous nested function. Вместо конкретного значения наша функция GetJ может вернуть хендл на собственную вложенную функцию. И при вызове функции по этому хендлу родительская функция выполняться не будет, но ее область видимости с посчитанными параметрами останется! Запишем:
function fhnd_J = GetJ4(a,b,x)
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));

fhnd_J=@(x) nf_GetJ(x);

    function out_val = nf_GetJ(x)
        out_val=-(phi_a*x^2+phi_b*sqrt(1-x^2));
    end
end

Теперь в возвращаемую переменную fhnd_J записывается хендл функции, позволяющей вычислить значение J для использованных параметров a и b, не вычисляя при этом переменные phi_a и phi_b, а используя значения, посчитанные при создании хендла. Далее наша оптимизация выглядит вот так:
fhnd_GetJ4=GetJ4(a,b,x);
optimal_x=fminbnd(fhnd_GetJ4,0,1,opt);

И выполняется такая оптимизация за 0,8 мс.

Результат
Мы получили такое же быстродействие, как при явном выносе вычислений в примере с GetJ2, но сохранили целостность функции и удобство ее использования.
В дальнейшем, если предполагается использовать функцию и для оптимизации, и просто для вычислений, то к ней можно добавить один необязательный параметр-флаг. Если он не задан, то функция возвращает значение. Если задан, то возвращается хендл на nested function.
Поделиться публикацией

Комментарии 12

    0
    Иными словами, вы описали использование замыканий в MATLAB. Наверное, стоило об этом упомянуть в самом начале для большей ясности.
      0
      Да, действительно. Но до вашего комментария я и не знал, что это так называется. Спасибо.
        0
        Это ни в коем случае не критика и не придирка :) Просто штука эта достаточно известная, а именно она и раскрывается в статье.
        Если для вас это новое, то советую почитать теорию про замыкания, функции высшего порядка и плюшки Lexical Scope — это в ряде случаев мощнейшие инструменты, позволяющие получать простые и красивые решения сложных задач. Думаю, это может помочь в вашей области деятельности.
      0
      Это все здорово, но не стоит забывать про высокие накладные расходы при использовании указателей на функции и замыкания в матлабе. Не зная этого, можно долго искать узкое место. Просто для примера попробуйте сравнить вызов в цикле встроенной функции, скажем sin и её же через указатель @sin. Разница будет раз в 20. Могу найти свой старый бенчмарк, в котором оцениваются накладные расходы на вызов функций и методов при разных условиях, в том числе при использовании замыканий. Например, при использовании ООП накладные расходы просто огромны, о чём я писал в комментарии к недавней статье про матлаб.
        0
        Да, конечно, накладные расходы высоки. Но задач, когда экономия в разы выше — множество. И не только таких надуманных, как в статье, но и вполне жизненных. В частности, эта статья навеяна рабочей ситуацией — занимались планированием оптимального движения избыточного манипулятора вдоль заданного пути, и я разбирал чужой матлабовский код. За счет введения таких замыканий существенно повысили быстродействие.
          0
          А между вызовом встроенной и пользовательской функцией разница заметна? Сейчас проверить не могу, но предпологаю, что пользовательские функции будут проигрывать замыканиям не столь сильно.
            0
            Примерно как-то так получается. Функции вызывались 10000 раз в цикле. Это время и замерялось.

            По функциям:
            0 - вызов функции sin как есть sin(x)
            1 - вызов функции sin через указатель @ sin
            2 - вызов функции sin через обёртку из анонимной функции @(x) sin(x)
            3 - вызов функции sin через обёртку из анонимной и nested функций @(x) sin1(x)
            4 - вызов функции sin через обёртку из анонимной и nested функции с передачей параметра, доступного в анонимной функции @(x) sin2(x, p)
            5 - вызов функции sin через обёртку из анонимной и nested функции, в которой используется параметр, доступный в nested функции @(x) sin3(x)
            



              0
              Thanks!

              Странно, что 1 медленнее чем 2.

              То, что 5 так тормозит, неприятно, но будем учитывать.
                0
                Странно, что 1 медленнее чем 2.

                Это меня тоже сильно удивляет, но факт остаётся фактом: прямой указатель на встроенную функцию гораздо медленнее, чем указатель на обёртку из анонимной функции.
                0
                Хм. На моем i7-2600@3.4GHz, 16GB RAM, Matlab R2012a это выполняется почти в 10 раз быстрее. Может, там цикл на 100'000 раз? Кстати, в R2013a еще примерно на 5% быстрее.
                Можно посмотреть код sin2 и sin3?
                  0
                  Да, цикл выполнялся 100000, потерялся один нолик.

                  Код функций простейший:
                  function y = sin2(x, a)
                      y = sin(x);
                      a;
                  end
                  
                  function y = sin3(x)
                       % nested-функция, параметр p определён во внешней функции
                       y = sin(x);
                       p;
                  end
                  

                  Как видно, параметры даже явно не используются, просто присутствуют в теле функции, что значительно замедляет их работу.
            0
            Спасибо, прекрасное решение (я тоже не знал слова «замыкание» — спасибо и за него)!

            Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

            Самое читаемое