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

Открыть сценарий в реальном времени

В этом примере показано, что инструментарий символьной математики помогает минимизировать ошибки при решении нелинейной системы уравнений.

В этом примере используются следующие возможности инструментария символьной математики:

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

Цель - найти минимум функции Розенброка, или так называемой «банановой» функции.

$f (x)_{}^{} {=∑i=1n/2100}_{(}_{}^{}^{} x2i-x2i-12) 2_{+ (}^{}$ 1-x2i-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

Figure contains an axes. The axes contains an object of type functionsurface.

Мотивация

Пример «Нелинейные уравнения с аналитическим якобианом» в панели инструментов оптимизации решает ту же задачу с помощью 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{array}{c} 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{array})$

Чтобы сопоставить промежуточные результаты, показанные в примере панели инструментов оптимизации, используйте ту же форму 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{array}{c} 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{array})$

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

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

JF = jacobian(F,x);

Отображение первых 10 строк и столбцов матрицы якобиана.

JF(1:10,1:10)

ans = 
 $(\begin{array}{c} - 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{array})$

Большинство элементов матрицы Якобиана 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' кому'iter' для просмотра хода выполнения решателя. Набор 'Jacobian' кому 'on' для использования якобиана, определенного в 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          212.48         1.5625           34.1            1.56
     6          7         116.793       0.390625           5.79           0.391
     7          8         109.947       0.390625          0.753           0.391
     8          9         99.4696       0.976563            1.2           0.977
     9         10         83.6416        2.44141           7.13            2.44
    10         11         77.7663        2.44141           9.94            2.44
    11         12         77.7663        2.44141           9.94            2.44
    12         13          43.013       0.610352           1.38            0.61
    13         14         36.4334       0.610352           1.58            0.61
    14         15         34.1448        1.52588           6.71            1.53
    15         16         18.0108        1.52588           4.91            1.53
    16         17         8.48336        1.52588           3.74            1.53
    17         18         3.74566        1.52588           3.58            1.53
    18         19         1.46166        1.52588           3.32            1.53
    19         20         0.29997        1.24265           1.94            1.53
    20         21               0      0.0547695              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 функция, которая является числовым решателем, доступным в инструментарии символьной математики. Ключевое различие между 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{array}{c} 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{array})$

Вспомогательные функции

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

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

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

Мотивация

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

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

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

Использовать `fsolve` Численное решение системы нелинейных уравнений

Использовать vpasolve для численного решения системы нелинейных уравнений

Вспомогательные функции

Документация по инструментам символьной математики

Поддержка

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

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

Мотивация

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

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

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

Использовать fsolve Численное решение системы нелинейных уравнений

Использовать vpasolve для численного решения системы нелинейных уравнений

Вспомогательные функции

Документация по инструментам символьной математики

Поддержка

Использовать `fsolve` Численное решение системы нелинейных уравнений