Exponenta Banner

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

Скрипт Open Live Script

Общие рекомендации для настройки задачи наименьших квадратов должны сформулировать проблему способом, которая позволяет 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')

Проблема состоит в том, чтобы соответствовать функции

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,lam)c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t);

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

diffexpr = sum((fcn2optimexpr(diffun,c,lam) - 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.8852e-06
       iterations: 6
        funcCount: 35
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0096
          message: '...'
           solver: 'lsqnonlin'

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

resp = diffun(sol.c,sol.lam);
hold on
plot(t,resp)
hold off

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

Подход решения Используя 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: 185
         stepsize: 0.0017
     lssteplength: 1
    firstorderopt: 2.9562e-05
        algorithm: 'quasi-newton'
          message: '...'
           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 185 function evaluations using fminunc, and 35 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 функционируйте к выражению оптимизации. Чтобы избежать предупреждения, дайте выходной размер получившегося выражения. Создайте новую задачу оптимизации с целью как сумма различий в квадрате между конвертированным 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.4071e-06
       iterations: 10
        funcCount: 33
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0080
          message: '...'
           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 35 using the 4-d formulation.

Разделите проблема Более устойчива к Исходному предположению

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

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.0057
       iterations: 32
        funcCount: 165
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0021
          message: '...'
           solver: 'lsqnonlin'


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

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.

Смотрите также

|

Похожие темы

Документация Optimization Toolbox
Поддержка