Gauss Newton Method, How to fix Singular Matrix (or Prevent it).
16 次查看(过去 30 天)
显示 更早的评论
I am trying to fit a curve to a set of data using a summation of gaussians. I first used the matlab fit function to see the solution I am trying to get. I fitted a summation of 8 gaussians. I then went on to use the Gauss Netwon Method to try to obtain the same solution but I am ending up with a singular matrix since the partial derivatives mostly go to 0 or very close to it. Is there any way to fix this?
close all; clear all; clc
% DATA LOADING
D=dlmread('FinalProject.csv',',',1,0);
x=D(:,1);
y=D(:,2);
figure
plot(x,y,'.')
f=fit(x,y,'gauss8')
hold on
plot(f,x,y)
Tol = .1;
L2Norm = 1;
a = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
while L2Norm > Tol
for i = 1:411
F(i) = a(1)*exp(-((x(i)-a(2))/a(3))^2) + a(4)*exp(-((x(i)-a(5))/a(6))^2) + a(7)*exp(-((x(i)-a(8))/a(9))^2) + a(10)*exp(-((x(i)-a(11))/a(12))^2) + a(13)*exp(-((x(i)-a(14))/a(15))^2) + a(16)*exp(-((x(i)-a(17))/a(18))^2) + a(19)*exp(-((x(i)-a(20))/a(21))^2) + a(22)*exp(-((x(i)-a(23))/a(24))^2);
r(i) = F(i) - y(i);
Dr(i,1) = exp(-((a(2) - x(i)))^2/(a(3)^2));
Dr(i,2) = exp(-(a(5) - x(i))^2/a(6)^2);
Dr(i,3) = exp(-(a(8) - x(i))^2/a(9)^2);
Dr(i,4) = exp(-(a(11) - x(i))^2/a(12)^2);
Dr(i,5) = exp(-(a(14) - x(i))^2/a(15)^2);
Dr(i,6) = exp(-(a(17) - x(i))^2/a(18)^2);
Dr(i,7) = exp(-(a(20) - x(i))^2/a(21)^2);
Dr(i,8) = exp(-(a(23) - x(i))^2/a(24)^2);
Dr(i,9) = -(a(1)*exp(-(a(2) - x(i))^2/a(3)^2)*(2*a(2) - 2*x(i)))/a(3)^2;
Dr(i,10) = -(a(4)*exp(-(a(5) - x(i))^2/a(6)^2)*(2*a(5) - 2*x(i)))/a(6)^2;
Dr(i,11) = -(a(7)*exp(-(a(8) - x(i))^2/a(9)^2)*(2*a(8) - 2*x(i)))/a(9)^2;
Dr(i,12) = -(a(10)*exp(-(a(11) - x(i))^2/a(12)^2)*(2*a(11) - 2*x(i)))/a(12)^2;
Dr(i,13) = -(a(13)*exp(-(a(14) - x(i))^2/a(15)^2)*(2*a(14) - 2*x(i)))/a(15)^2;
Dr(i,14) = -(a(16)*exp(-(a(17) - x(i))^2/a(18)^2)*(2*a(17) - 2*x(i)))/a(18)^2;
Dr(i,15) = -(a(19)*exp(-(a(20) - x(i))^2/a(21)^2)*(2*a(20) - 2*x(i)))/a(21)^2;
Dr(i,16) = -(a(22)*exp(-(a(23) - x(i))^2/a(24)^2)*(2*a(23) - 2*x(i)))/a(24)^2;
Dr(i,17) = (2*a(1)*exp(-(a(2) - x(i))^2/a(3)^2)*(a(2) - x(i))^2)/a(3)^3;
Dr(i,18) = (2*a(4)*exp(-(a(5) - x(i))^2/a(6)^2)*(a(5) - x(i))^2)/a(6)^3;
Dr(i,19) = (2*a(7)*exp(-(a(8) - x(i))^2/a(9)^2)*(a(8) - x(i))^2)/a(9)^3;
Dr(i,20) = (2*a(10)*exp(-(a(11) - x(i))^2/a(12)^2)*(a(11) - x(i))^2)/a(12)^3;
Dr(i,21) = (2*a(13)*exp(-(a(14) - x(i))^2/a(15)^2)*(a(14) - x(i))^2)/a(15)^3;
Dr(i,22) = (2*a(16)*exp(-(a(17) - x(i))^2/a(18)^2)*(a(17) - x(i))^2)/a(18)^3;
Dr(i,23) = (2*a(19)*exp(-(a(20) - x(i))^2/a(21)^2)*(a(20) - x(i))^2)/a(21)^3;
Dr(i,24) = (2*a(22)*exp(-(a(23) - x(i))^2/a(24)^2)*(a(23) - x(i))^2)/a(24)^3;
end
LHS = transpose(Dr) * Dr;
RHS = transpose(Dr) * transpose(r);
da = - (inv(LHS) * RHS);
L2Norm = 0;
for j = 1:24
S = da(j)^2;
L2Norm = S + L2Norm;
end
L2Norm = sqrt(L2Norm);
i = 1;
j = 1;
a = a + transpose(da);
end
0 个评论
回答(1 个)
Nipun
2024-5-22
Hi James,
I understand that you are trying to fit a curve to your data using a summation of Gaussians and are encountering issues with a singular matrix when applying the Gauss-Newton method due to partial derivatives approaching zero. This can indeed happen when the Jacobian matrix (which you're computing as Dr) becomes ill-conditioned or nearly singular, making its inverse either impossible to compute accurately or leading to very large errors in the computed parameters.
I recommend "regularizing" your code. Add a regularization term to the diagonal of your Jacobian or Hessian matrix to improve its condition. This is a common approach to deal with ill-conditioned matrices. A simple form of regularization involves adding a small constant ( \lambda ) to the diagonal elements of the matrix ( LHS ) before inverting it.
Here is a suggestion for incorporating a simple form of regularization into your existing code. Note that this is a basic approach and might need adjustments based on your specific problem:
lambda = 1e-3; % Regularization parameter, adjust based on your problem
while L2Norm > Tol
% Your existing code to compute LHS, RHS, and r goes here
% Modify LHS with regularization
LHS_reg = LHS + lambda * eye(size(LHS));
% Use LHS_reg instead of LHS for computing da
da = - (inv(LHS_reg) * RHS);
% Your existing code to update a and compute L2Norm goes here
end
This modification adds a small value to the diagonal elements of the LHS matrix, which can help to make it invertible and stabilize the Gauss-Newton iteration. The regularization parameter ( \lambda ) may need to be adjusted based on the scale of your problem and the degree of ill-conditioning.
Hope this helps.
Regards,
Nipun
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!