Нелинейная подгонка данных с использованием нескольких основанных на проблеме подходов

Общий совет для настройки задачи наименьших квадратов состоит в том, чтобы сформулировать задачу таким образом, чтобы это позволяло solve чтобы распознать, что задача имеет форму наименьших квадратов. Когда ты это делаешь, solve внутренние вызовы lsqnonlin, что эффективно при решении задач наименьших квадратов. Смотрите Запись Целевой Функции для основанных на проблеме наименьших квадратов.

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

Setup задачи

Примите во внимание следующие данные:

Data = ...
  [0.0000    5.8955
   0.1000    3.5639
   0.2000    2.5173
   0.3000    1.9790
   0.4000    1.8990
   0.5000    1.3938
   0.6000    1.1359
   0.7000    1.0096
   0.8000    1.0343
   0.9000    0.8435
   1.0000    0.6856
   1.1000    0.6100
   1.2000    0.5392
   1.3000    0.3946
   1.4000    0.3903
   1.5000    0.5474
   1.6000    0.3459
   1.7000    0.1370
   1.8000    0.2211
   1.9000    0.1704
   2.0000    0.2636];

Постройте график точек данных.

t = Data(:,1);
y = Data(:,2);
plot(t,y,'ro')
title('Data points')

Figure contains an axes. The axes with title Data points contains an object of type line.

Задача состоит в том, чтобы соответствовать функции

y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)

к данным.

Подход к решению с использованием решателя по умолчанию

Для начала задайте переменные оптимизации, соответствующие уравнению.

c = optimvar('c',2);
lam = optimvar('lam',2);

Произвольно установите начальную точку x0 следующим образом: c(1) = 1, c(2) = 1, lam(1) = 1, и lam(2) = 0:

x0.c = [1,1];
x0.lam = [1,0];

Создайте функцию, которая вычисляет значение отклика в разы t когда параметры c и lam.

diffun = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t);

Преобразование diffun к выражению оптимизации, которое суммирует квадраты различий между функцией и данными y.

diffexpr = sum((diffun - y).^2);

Создайте задачу оптимизации, имеющую diffexpr как целевая функция.

ssqprob = optimproblem('Objective',diffexpr);

Решите задачу с помощью решателя по умолчанию.

[sol,fval,exitflag,output] = solve(ssqprob,x0)
Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
sol = struct with fields:
      c: [2x1 double]
    lam: [2x1 double]

fval = 0.1477
exitflag = 
    FunctionChangeBelowTolerance

output = struct with fields:
          firstorderopt: 7.8870e-06
             iterations: 6
              funcCount: 7
           cgiterations: 0
              algorithm: 'trust-region-reflective'
               stepsize: 0.0096
                message: '...'
    objectivederivative: "forward-AD"
                 solver: 'lsqnonlin'

Постройте график получившейся кривой на основе возвращенных значений решения sol.c и sol.lam.

resp = evaluate(diffun,sol);
hold on
plot(t,resp)
hold off

Figure contains an axes. The axes with title Data points contains 2 objects of type line.

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

Подход к решению с использованием fminunc

Чтобы решить проблему с помощью fminunc решатель, установите 'Solver' опция для 'fminunc' при вызове solve.

[xunc,fvalunc,exitflagunc,outputunc] = solve(ssqprob,x0,'Solver',"fminunc")
Solving problem using fminunc.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
xunc = struct with fields:
      c: [2x1 double]
    lam: [2x1 double]

fvalunc = 0.1477
exitflagunc = 
    OptimalSolution

outputunc = struct with fields:
             iterations: 30
              funcCount: 37
               stepsize: 0.0017
           lssteplength: 1
          firstorderopt: 2.9454e-05
              algorithm: 'quasi-newton'
                message: '...'
    objectivederivative: "forward-AD"
                 solver: 'fminunc'

Заметьте, что fminunc найдено то же решение, что и lsqcurvefitно для этого потребовалось еще много вычисления функции. Параметры для fminunc находятся в противоположном порядке, как для lsqcurvefit; больший lam является lam(2), не lam(1). Это неудивительно, порядок переменных произволен.

fprintf(['There were %d iterations using fminunc,' ...
    ' and %d using lsqcurvefit.\n'], ...
    outputunc.iterations,output.iterations)
There were 30 iterations using fminunc, and 6 using lsqcurvefit.
fprintf(['There were %d function evaluations using fminunc,' ...
    ' and %d using lsqcurvefit.'], ...
    outputunc.funcCount,output.funcCount)
There were 37 function evaluations using fminunc, and 7 using lsqcurvefit.

Разделение линейной и нелинейной задач

Заметьте, что задача аппроксимации линейна в параметрах c(1) и c(2). Это означает для любых значений lam(1) и lam(2)можно использовать оператор обратной косой черты, чтобы найти значения c(1) и c(2) которые решают задачу наименьших квадратов.

Переделайте задачу как двумерную задачу, ища лучшие значения lam(1) и lam(2). Значения c(1) и c(2) вычисляются на каждом этапе с использованием оператора обратной косой черты, как описано выше. Для этого используйте fitvector функция, которая выполняет операцию обратной косой черты для получения c(1) и c(2) при каждой итерации решателя.

type fitvector
function yEst = fitvector(lam,xdata,ydata)
%FITVECTOR  Used by DATDEMO to return value of fitting function.
%   yEst = FITVECTOR(lam,xdata) returns the value of the fitting function, y
%   (defined below), at the data points xdata with parameters set to lam.
%   yEst is returned as a N-by-1 column vector, where N is the number of
%   data points.
%
%   FITVECTOR assumes the fitting function, y, takes the form
%
%     y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t)
%
%   with n linear parameters c, and n nonlinear parameters lam.
%
%   To solve for the linear parameters c, we build a matrix A
%   where the j-th column of A is exp(-lam(j)*xdata) (xdata is a vector).
%   Then we solve A*c = ydata for the linear least-squares solution c,
%   where ydata is the observed values of y.

A = zeros(length(xdata),length(lam));  % build A matrix
for j = 1:length(lam)
   A(:,j) = exp(-lam(j)*xdata);
end
c = A\ydata; % solve A*c = y for linear parameters c
yEst = A*c; % return the estimated response based on c

Решите задачу, используя solve начиная с двумерной начальной точки x02.lam = [1,0]. Для этого сначала преобразуйте fitvector функцию в выражение оптимизации с использованием fcn2optimexpr. См. Преобразование нелинейной функции в выражение оптимизации. Чтобы избежать предупреждения, задайте выходной размер получившегося выражения. Создайте новую задачу оптимизации с целью как суммой квадратов различий между преобразованными fitvector функция и данные y.

x02.lam = x0.lam;
F2 = fcn2optimexpr(@(x) fitvector(x,t,y),lam,'OutputSize',[length(t),1]);
ssqprob2 = optimproblem('Objective',sum((F2 - y).^2));
[sol2,fval2,exitflag2,output2] = solve(ssqprob2,x02)
Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
sol2 = struct with fields:
    lam: [2x1 double]

fval2 = 0.1477
exitflag2 = 
    FunctionChangeBelowTolerance

output2 = struct with fields:
          firstorderopt: 4.4018e-06
             iterations: 10
              funcCount: 33
           cgiterations: 0
              algorithm: 'trust-region-reflective'
               stepsize: 0.0080
                message: '...'
    objectivederivative: "finite-differences"
                 solver: 'lsqnonlin'

Эффективность двумерного решения аналогична эффективности четырехмерного решения:

fprintf(['There were %d function evaluations using the 2-d ' ...
    'formulation, and %d using the 4-d formulation.'], ...
    output2.funcCount,output.funcCount)
There were 33 function evaluations using the 2-d formulation, and 7 using the 4-d formulation.

Задача разделения более устойчива к первоначальным предположениям

Выбор плохой начальной точки для исходной четырехпараметрической задачи приводит к локальному решению, которое не является глобальным. Выбор начальной точки с такими же плохими lam(1) и lam(2) значения для разделенной двухпараметрической задачи приводят к глобальному решению. Чтобы показать это, перезапустите исходную задачу с начальной точкой, которая приводит к относительно плохому локальному решению, и сравните полученную подгонку с глобальным решением.

x0bad.c = [5 1];
x0bad.lam = [1 0];
[solbad,fvalbad,exitflagbad,outputbad] = solve(ssqprob,x0bad)
Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
solbad = struct with fields:
      c: [2x1 double]
    lam: [2x1 double]

fvalbad = 2.2173
exitflagbad = 
    FunctionChangeBelowTolerance

outputbad = struct with fields:
          firstorderopt: 0.0036
             iterations: 31
              funcCount: 32
           cgiterations: 0
              algorithm: 'trust-region-reflective'
               stepsize: 0.0012
                message: '...'
    objectivederivative: "forward-AD"
                 solver: 'lsqnonlin'

respbad = evaluate(diffun,solbad);
hold on
plot(t,respbad,'g')
legend('Data','Global fit','Bad local fit','Location','NE')
hold off

Figure contains an axes. The axes with title Data points contains 3 objects of type line. These objects represent Data, Global fit, Bad local fit.

fprintf(['The residual norm at the good ending point is %f,' ...
   ' and the residual norm at the bad ending point is %f.'], ...
   fval,fvalbad)
The residual norm at the good ending point is 0.147723, and the residual norm at the bad ending point is 2.217300.

См. также

|

Похожие темы

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