Code shows matrix is singular to working precision. How can I fix it?
2 views (last 30 days)
Show older comments
x=[-2 0 1 2 3]
y=[2.5,1,2.5,5.7,11.6]
A=[x.^0 x.^2]
alin=A'*A\A'*log(y)
a=[exp(alin(1));alin(2)]
r=a(1)*exp(a(2)*x.^2)-y
RMSE=norm(r,2)/sqrt(5)
0 Comments
Answers (1)
Jon
on 13 Oct 2023
Consider your matrix A'A. If you take a column matrix, in your case A' and multiply it by a row, in your case, A, then every column of A'A will be just a scaled version of the column A', so the columns will be linearly dependent and A'A is therefore singular
A'A\A' tries to solve the equation A'A x = A' and fails because A'A is singular.
I'm not really sure what you are trying to accomplish, and whether it makes any sense, but if you want to get a solution to A'Ax=A' you can use x = pinv(A'*A)*A'. Alternatively you could just recognize that by the construction of A'A, a solution of A'Ax = A' would be [1;0;0;0;0;0;0;0;0;0]
2 Comments
Askic V
on 13 Oct 2023
Edited: Askic V
on 13 Oct 2023
If r means residual and you calculate root mean square error (RMSE), it seems to me that you're trying to find coefficients a1 and a2 such that minimize error fit of the following model a(1)*exp(a(2)*x.^2) to your data.
If that is the case, I would determine coefficients a1 and a2 in the following way:
x = [-2 0 1 2 3];
y = [2.5,1,2.5,5.7,11.6];
% Define the model function as a function handle
model = @(a, x) a(1) * exp(a(2) * x.^2);
% Initial guess for the coefficients a1 and a2
initialGuess = [2, 0.5];
% Perform nonlinear curve fitting using lsqcurvefit
options = optimoptions('lsqcurvefit', 'Display', 'off');
A = lsqcurvefit(model, initialGuess, x, y, [], [], options);
a1 = A(1)
a2 = A(2)
r = A(1)*exp(A(2)*x.^2) - y;
RMSE = norm(r,2)/sqrt(numel(r))
Jon
on 13 Oct 2023
You can solve this using linear least squares as follows, this looks like what you were attempting to do
% Data as vectors
x = [-2 0 1 2 3]'
y = [2.5,1,2.5,5.7,11.6]'
% solve for coefficients, to best
% fit y = a(1)*exp(a(2)*x.^2)
% taking log of both side
% log(y) = log(a(1)) + a(2)*x.^2
% log(y) = c(1) + c(2)*x.^2, a(1) = exp(c(1))
% find least squares solution
A = [x.^0 x.^2]
c = A\log(y);
a = [exp(c(1));c(2)]
% compute residual
r = a(1)*exp(a(2)*x.^2) - y
RMSE =norm(r,2)/sqrt(5)
Note that this gives a larger residual, than the non-linear approach listed above. I think this is because we are minimizing the error between log(yfit) - log(y), rather than yfit - y, on the other hand the linear solution is very reliable and doesn't require an initial guess
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!