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

Этот пример демонстрирует, что 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

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