Gauss Newton Method, How to fix Singular Matrix (or Prevent it).

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 

Answer 1

[Nipun]

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