Улучшение точности и производительности в оптимизации

Этот пример демонстрирует, что Symbolic Math Toolbox помогает минимизировать ошибки при решении нелинейной системы уравнений.

Этот пример использует следующую поддержку Symbolic Math Toolbox:

Вычисление аналитического градиента с помощью gradient
Вычисление аналитического якобиана с помощью jacobian
Преобразование символьных результатов в числовые функции с помощью matlabFunction
Аналитическая визуализация с fplot

Цель состоит в том, чтобы найти минимум Функции Розенброка или так называемую "банановую" функцию.

$f (x) = \sum_{i = 1}^{n / 2} 100 (x_{2 i} - x_{2 i - 1}^{2})^{2} + (1 - x_{2 i - 1})^{2}$

Используйте fplot создать быструю визуализацию Функции Розенброка двух переменных

syms x y
a=1;    b=100;
f(x,y)=(a-x)^2+b*(y-x^2)^2

f(x, y) =  ${(x - 1)}^{2} + 100 {(y - x^{2})}^{2}$

fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
caxis([0,1000])
colorbar

Мотивация

“Нелинейные уравнения с Аналитическим якобиевским” примером в Optimization Toolbox решают ту же задачу при помощи fsolve функция. Рабочий процесс, показанный в том примере, имеет два потенциальных источника ошибок. Вы должны были бы

Переведите уравнения для Функции Розенброка и якобиана от математической формы в тексте к цифровому коду.
Явным образом вычислите якобиан. Для сложных уравнений эта задача подвержена ошибкам.

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

Перепишите функцию Розенброка как систему нелинейных уравнений

Во-первых, преобразуйте Функцию Розенброка f к системе нелинейных уравнений F, где F градиент f .

clear
n = 64;
x = sym('x',[n,1]);
f = sum(100*(x(2:2:n)-x(1:2:n).^2).^2 + (1-x(1:2:n)).^2);
F = gradient(f,x);

Покажите первые 10 уравнений:

F(1:10)

ans = 
 $(\begin{matrix} 2 x_{1} - 400 x_{1} (x_{2} - {x_{1}}^{2}) - 2 \\ 200 x_{2} - 200 {x_{1}}^{2} \\ 2 x_{3} - 400 x_{3} (x_{4} - {x_{3}}^{2}) - 2 \\ 200 x_{4} - 200 {x_{3}}^{2} \\ 2 x_{5} - 400 x_{5} (x_{6} - {x_{5}}^{2}) - 2 \\ 200 x_{6} - 200 {x_{5}}^{2} \\ 2 x_{7} - 400 x_{7} (x_{8} - {x_{7}}^{2}) - 2 \\ 200 x_{8} - 200 {x_{7}}^{2} \\ 2 x_{9} - 400 x_{9} (x_{10} - {x_{9}}^{2}) - 2 \\ 200 x_{10} - 200 {x_{9}}^{2} \end{matrix})$

Чтобы совпадать с промежуточными результатами, показанными в примере Optimization Toolbox, используйте ту же форму F.

F(1:2:n) = simplify(F(1:2:n) + 2*x(1:2:n).*F(2:2:n));
F(1:2:n) = -F(1:2:n)/2;
F(2:2:n) = F(2:2:n)/20;

Теперь системы уравнений идентичны:

F(1:10)

ans = 
 $(\begin{matrix} 1 - x_{1} \\ 10 x_{2} - 10 {x_{1}}^{2} \\ 1 - x_{3} \\ 10 x_{4} - 10 {x_{3}}^{2} \\ 1 - x_{5} \\ 10 x_{6} - 10 {x_{5}}^{2} \\ 1 - x_{7} \\ 10 x_{8} - 10 {x_{7}}^{2} \\ 1 - x_{9} \\ 10 x_{10} - 10 {x_{9}}^{2} \end{matrix})$

Вычислите якобиан матрицы, представляющей систему уравнений

Используйте jacobian вычислить якобиан F . Эта функция вычисляет якобиан символически, таким образом избегая ошибок, сопоставленных числовыми приближениями производных.

JF = jacobian(F,x);

Покажите первые 10 строк и столбцов якобиевской матрицы.

JF(1:10,1:10)

ans = 
 $(\begin{matrix} - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - 20 x_{1} & 10 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 20 x_{3} & 10 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 20 x_{5} & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{7} & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{9} & 10 \end{matrix})$

Большинство элементов якобиевского матричного JF нули. Тем не менее, когда вы преобразуете эту матрицу в функцию MATLAB, результатом является плотная числовая матрица. Часто, операции на разреженных матрицах более эффективны, чем те же операции на плотных матрицах.

Поэтому, чтобы создать эффективный код, преобразуйте только ненулевые элементы JF к символьному вектору прежде, чем вызвать matlabFunction. is и js представляйте шаблон разреженности JF.

[is,js] = find(JF);
JF = JF(JF~=0);

Преобразуйте Систему уравнений и якобиан к функции MATLAB

Система уравнений F, представление Функции Розенброка, символьная матрица, которая состоит из символьных выражений. Смочь решить его с fsolve функция, преобразуйте эту систему в функцию MATLAB. На этом шаге удобно преобразовать оба F и его якобиан, JF, к одной основанной на файле функции MATLAB, FJFfun. В принципе это позволяет F и JF к переменным повторного использования. В качестве альтернативы можно преобразовать их в отдельные функции MATLAB.

matlabFunction([F;JF],'var',x,'file','FJFfun');

FJFfun принимает переменные как список, но fsolve ожидает вектор переменных. Функциональный genericObj преобразует вектор в список. Анонимная функция bananaObj задан, чтобы зафиксировать значения аргументов для genericObj. Отметьте использование шаблона разреженности в функциональном genericObj. Далее обратите внимание что этот подход использования FJFfun и genericObj вместе не связывается к конкретному выражению, представленному F и может быть использован как есть для любого F.

bananaobj = @(y) genericObj(y,@FJFfun,is,js)

bananaobj = function_handle with value:
    @(y)genericObj(y,@FJFfun,is,js)

Используйте fsolve, чтобы Численно Решить Систему Нелинейных уравнений

Используйте fsolve для системы уравнений, преобразованной в функцию MATLAB. Используйте начальную точку x (i) = –1.9 для нечетных индексов и x (i) = 2 для ровных индексов. Установите 'Display' to'iter' видеть прогресс решателя. Установите 'Jacobian' к 'on' использовать якобиан задало in bananaobj.

x0(1:n,1) = -1.9;
x0(2:2:n,1) = 2;
options = optimoptions(@fsolve,'Display','iter','Jacobian','on');
[sol1,Fsol,exitflag,output,JAC] = fsolve(bananaobj,x0,options);

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          1         8563.84                           615               1
     1          2         3093.71              1            329               1
     2          3         225.104            2.5           34.8             2.5
     3          4          212.48           6.25           34.1            6.25
     4          5          212.48           6.25           34.1            6.25
     5          6         102.771         1.5625           6.39            1.56
     6          7         102.771        3.90625           6.39            3.91
     7          8         87.7443       0.976563           2.19           0.977
     8          9         74.1426        2.44141           6.27            2.44
     9         10         74.1426        2.44141           6.27            2.44
    10         11          52.497       0.610352           1.52            0.61
    11         12         41.3297        1.52588           4.63            1.53
    12         13         34.5115        1.52588           6.97            1.53
    13         14         16.9716        1.52588           4.69            1.53
    14         15         8.16797        1.52588           3.77            1.53
    15         16         3.55178        1.52588           3.56            1.53
    16         17         1.38476        1.52588           3.31            1.53
    17         18        0.219553        1.16206           1.66            1.53
    18         19               0      0.0468565              0            1.53

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Используйте vpasolve, чтобы Численно Решить Систему Нелинейных уравнений

Система нелинейных уравнений, F, может быть альтернативно решен с vpasolve функция, которая является числовым решателем, доступным в Symbolic Math Toolbox. Основное отличие between vpasolve и fsolve тот fsolve использование арифметика с двойной точностью, чтобы решить уравнения, тогда как vpasolve арифметика переменной точности использования и поэтому, позволяет вам управлять точностью расчетов.

sol2 = vpasolve(F);

Если вы присваиваете решение системы уравнений к одной переменной, sol2, затем vpasolve возвращает решения в форме массива структур. Можно получить доступ к каждому решению (каждое поле массива структур):

sol2.x1

ans =  $1.0$

Также можно преобразовать sol2 к символьному вектору, и затем получают доступ к области значений решений. Например, отобразите 5-е к 25-м решениям:

for k=(1:64)
  sol2_array(k) = sol2.(char(x(k)));
end
sol2_array(5:25)

ans =  $(\begin{matrix} 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{matrix})$

Функции помощника

function [F,J] = genericObj(x,FJFfun,is,js)
% genericObj takes as input, vector x, Function and Jacobian evaluation 
% function handle FJFfun, and sparsity pattern is,js and returns as output,
% the numeric values of the Function and Jacobian: F and Jfunction 
% FJs(1:length(x)) is the numeric vector for the Function
% FJs(length(x)+1:end) is the numeric vector for the non-zero entries of the Jacobian

xcell = num2cell(x);
FJs = FJFfun(xcell{:});
F = FJs(1:length(x));
J = sparse(is,js,FJs(length(x)+1:end));

end

Документация