Общий совет для настройки задачи наименьших квадратов состоит в том, чтобы сформулировать проблему таким образом, чтобы это позволяло solve признать, что проблема имеет форму наименьших квадратов. Когда вы это делаете, solve внутренние вызовы lsqnonlin, что эффективно при решении задач наименьших квадратов. См. раздел Целевая функция записи для наименьших квадратов, основанных на проблемах.
Этот пример показывает эффективность решателя наименьших квадратов, сравнивая производительность lsqnonlin с тем, что fminunc по той же проблеме. Кроме того, в примере показаны дополнительные преимущества, которые можно получить, явно распознав и обрабатывая отдельно линейные части задачи.
Рассмотрим следующие данные:
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(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
Посадка выглядит как можно лучше.
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 fitvectorfunction 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
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.