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

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

Решаемые уравнения - $${}_{\mathrm{Fi}}(x)\mathrm{}$$= $$0,$$1≤i≤n. В примере $$используется\mathrm{}\mathrm{n}\mathrm{}\mathrm{=}$$1000.

Эта целевая функция достаточно проста, чтобы можно было вычислить ее Jacobian аналитически. Как поясняется в документе Writing Vector and Matrix Objective Functions, якобиановым $$J(x$$) системы уравнений $$F($$x) $${\mathrm{является}}_{}\mathrm{Jij}(\frac{x{)}_{}}{=\partial Fi{}_{(}}$$x) ∂xj. Предоставьте эту производную как второй вывод целевой функции. nlsf1 вспомогательная функция в конце этого примера создает целевую функцию $$F(x$$) и ее Jacobian $$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)

        7007

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)

     6

disp(output4.cgiterations)

     8

'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)

     6

disp(output5.cgiterations)

    19

Остаточная норма на этот раз немного ниже, а количество оценок функций не изменилось. Число итераций сопряженного градиента увеличивается с 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)

     6

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

Подведите итоги.

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

Связанные темы

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