Большая разреженная система нелинейных уравнений с якобианом

Открыть Live Script

В этом примере показано, как использовать функции fsolve решатель для эффективного решения больших разреженных систем уравнений. В примере используется целевая функция, заданная для системы $n$ уравнения,

$\begin{array}{l} F (1) = 3 x_{1} - 2 x_{1}^{2} - 2 x_{2} + 1, \\ F (i) = 3 x_{i} - 2 x_{i}^{2} - x_{i - 1} - 2 x_{i + 1} + 1, \\ F (n) = 3 x_{n} - 2 x_{n}^{2} - x_{n - 1} + 1 . \end{array}$

Уравнения, которые нужно решить, $F_{i} (x) = 0$ , $1 \leq i \leq n$ . Пример использует $n = 1000$ .

Эта целевая функция достаточно проста, чтобы можно было вычислить ее якобиан аналитически. Как объяснено в Writing Vector and Matrix Objective Functions, якобиан $J (x)$ системы уравнений $F (x)$ является $J_{i j} (x) = \frac{\partial F_{i} (x)}{\partial x_{j}}$ . Предоставьте эту производную как второй выход вашей целевой функции. The nlsf1 helper функция в конце этого примера создает целевую функцию $F (x)$ и его якобиан $J (x)$ .

Решить систему уравнений можно используя опции по умолчанию, которые называют 'trust-region-dogleg' алгоритм. Начните с точки xstart(i) = -1.

n = 1000;
xstart = -ones(n,1);
fun = @nlsf1; 
[x,fval,exitflag,output] = fsolve(fun,xstart);

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.

Отображение качества решения и количества выполненных вычислений функции.

disp(norm(fval))

   2.8577e-13

disp(output.funcCount)

fsolve решает уравнение точно, но для этого требуется тысячи вычислений функции.

Решить уравнение используя якобиан как по умолчанию, так и по 'trust-region' алгоритмы.

options = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,fval2,exitflag2,output2] = fsolve(fun,xstart,options);

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.

options.Algorithm = 'trust-region';
[x3,fval3,exitflag3,output3] = fsolve(fun,xstart,options);

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.

disp([norm(fval2),norm(fval3)])

   1.0e-08 *

    0.0000    0.1065

disp([output2.funcCount,output3.funcCount])

     7     5

Оба алгоритма берут крошечную часть количества вычислений функции по сравнению с количеством, не используя якобиан. Алгоритм по умолчанию принимает несколько больше вычислений функции, чем 'trust-region' алгоритм, но алгоритм по умолчанию достигает более точного ответа.

Посмотрите, задает ли 'PrecondBandWidth' опция 1 меняет 'trust-region' ответ или эффективность. Этот параметр вызывает fsolve использовать тридиагональный предварительный кондиционер, который должен быть эффективным для этой тридиагональной системы уравнений.

options.PrecondBandWidth = 1;
[x4,fval4,exitflag4,output4] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

disp(norm(fval4))

   3.1185e-05

disp(output4.funcCount)

disp(output4.cgiterations)

The 'PrecondBandWidth' установка опции приводит к fsolve дать несколько менее точный ответ, измеренный нормой невязки. Количество вычислений функции несколько увеличивается с 5 до 6. Решатель имеет менее десяти сопряженных градиентных итераций как часть процесса решения.

Посмотрите, как хорошо fsolve выполняет с диагональным предварительным кондиционером.

options.PrecondBandWidth = 0;
[x5,fval5,exitflag5,output5] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

disp(norm(fval5))

   2.0057e-05

disp(output5.funcCount)

disp(output5.cgiterations)

На этот раз норма невязки немного ниже, и количество вычислений функции остается неизменным. Количество сопряженных итераций градиента увеличивается с 8 до 19, что указывает на то, что это 'PrecondBandWidth' установка приводит к тому, что решатель делает больше работы.

Решить уравнения используя 'levenberg-marquardt' алгоритм.

options = optimoptions('fsolve','SpecifyObjectiveGradient',true,'Algorithm','levenberg-marquardt');
[x6,fval6,exitflag6,output6] = fsolve(fun,xstart,options);

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.

disp(norm(fval6))

   7.5854e-15

disp(output6.funcCount)

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

Результирующие результаты.

t = table([norm(fval);norm(fval2);norm(fval3);norm(fval4);norm(fval5);norm(fval6)],...
    [output.funcCount;output2.funcCount;output3.funcCount;output4.funcCount;output5.funcCount;output6.funcCount],...
    'VariableNames',["Residual" "Fevals"],...
    'RowNames',["Default" "Default+Jacobian" "Trust-Region+Jacobian" "Trust-Region+Jacobian,BW=1" "Trust-Region+Jacobian,BW=0" "Levenberg-Marquardt+Jacobian"])

t=6×2 table
                                     Residual     Fevals
                                    __________    ______

    Default                         2.8577e-13     7007 
    Default+Jacobian                2.5886e-13        7 
    Trust-Region+Jacobian           1.0646e-09        5 
    Trust-Region+Jacobian,BW=1      3.1185e-05        6 
    Trust-Region+Jacobian,BW=0      2.0057e-05        6 
    Levenberg-Marquardt+Jacobian    7.5854e-15        6

Этот код создает nlsf1 функция.

function [F,J] = nlsf1(x)
% Evaluate the vector function
n = length(x);
F = zeros(n,1);
i = 2:(n-1);
F(i) = (3-2*x(i)).*x(i)-x(i-1)-2*x(i+1) + 1;
F(n) = (3-2*x(n)).*x(n)-x(n-1) + 1;
F(1) = (3-2*x(1)).*x(1)-2*x(2) + 1;
% Evaluate the Jacobian if nargout > 1
if nargout > 1
   d = -4*x + 3*ones(n,1); D = sparse(1:n,1:n,d,n,n);
   c = -2*ones(n-1,1); C = sparse(1:n-1,2:n,c,n,n);
   e = -ones(n-1,1); E = sparse(2:n,1:n-1,e,n,n);
   J = C + D + E;
end
end

См. также

fsolve

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

Большая разреженная система нелинейных уравнений с якобианом

См. также

Похожие темы

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

Поддержка