Warning: Matrix is singular, close to singular or badly scaled: Help with finding error in code

Greetings, I am working on a code for my numerical methods class and am having some trouble understanding why my code is not operating properly. I keep getting the warning "Matrix is singular, close to singular, or badly scaled.

The goal of my code is to vary Y to fit experimental data using a non-linear numerical model. I believe it will be a simple fix with the line: dY = (Z'*Z)\(Z'*D); but I am not sure what fix it would require.

I have checked the sizes of the matrixes, they are as follow: Z 9x3 D 9x1 dY 3x3

In my pseudo code dY should only be a 3x1. I have a basic understanding of matrix math but not enough to work through this with google searches and my know knowledge. Any ideas would be of help. Code is below.

Thanks in advance! Vinny Morasko

clear; clc;
clear all;

    %data
    T=60;
    Y=[0.0021;0.2628;0.5769];

Activity =
[1;0.951693618101212;0.787050622169703;0.703301583776796;0.611130056863694;0.430688946871704;0.217509576417032;0.058196064746889;0.007090117454854];
t=      [0;30;120;180;240;360;720;1440;2880];

%plot data
figure(1); clf(1)
plot(t,Activity,'o')

%iterate 
for iter=1:10;

      %build matricies 
      n=length(t);
      m=3;
      Z=zeros(n,m);
      D=zeros(n,1);

      for k=1:n

          Z(k,1)= ((Y(1)^2*t(k)-Y(1)*t(k)*Y(2)+Y(2)*exp(-Y(1)*t(k))*Y(3))/...
              ((Y(1)-Y(2))^2))-((exp(-t(k)*Y(2))*Y(2)*Y(3))/((Y(1)-Y(2))^2));

          Z(k,2)= (((Y(2)*t(k)-t(k)*Y(1)+1)*exp(-Y(2)*t(k))*Y(1)*Y(3))/...
              ((Y(2)-Y(1))^2))-((exp(-t(k)*Y(1))*Y(1)*Y(3))/((Y(2)-Y(1))^2));

          Z(k,3)= ((exp(-t(k)*Y(2))*Y(1))/(Y(1)-Y(2)))-((exp(-t(k)*Y(1))*...
              Y(1))/(Y(1)-Y(2)));

          D(k,1)= Activity(k)-(1+(Y(3)*Y(1))/(Y(2)-Y(1))*exp(-Y(1)*t(k)))-...
              ((Y(1)*Y(3))/(Y(2)-Y(1))*exp(-Y(2)*t(k)));
      end

      %solve for dY
      dY = (Z'*Z)\(Z'*D);

      %update Y's

      Y=Y+dY;

      %plot fit
      xf=linspace(min(t)-0.5,max(t+0.5),100);
      yf=(1+(Y(3)*Y(1))/(Y(2)-Y(1))*exp(-Y(1)*xf))-...
              ((Y(1)*Y(3))/(Y(2)-Y(1))*exp(-Y(2)*xf));

      hold on 
      plot (xf,yf)
      pause (0.5)
  end

Answer 1

[John D'Errico]

One hint:

If you knowhow to use backslash, then why are you using it as

      %solve for dY
      dY = (Z'*Z)\(Z'*D);

Instead, use backslash as it is designed to be used!

      %solve for dY
      dY = Z\D;

The solution that you are using squares the condition number. So if a linear system is already poorly scaled, then it is scaled far more poorly using the form you used.

Next, any time you have exponentials in a problem, expect things to be poorly scaled. Now, reread what I said above.

To learn more about your problem, try this:

cond(Z)

what does it tell you? You can learn even more about the problem by this:

svd(Z)

What does it show?

I cannot go further, since I cannot even test your code, lacking the data you are reading in.

[Vincent Morasko]

Thank you for the quick response. I have edited the code so you could run it yourself (added data). I am using backslash in that manner to create non-linear least squares code. A weighted residual analysis so to speak.

I used the cond(Z) and svd(Z) codes like you requested, the results are below. From what I gathered from the help section of MATLab, my matrix is poorly conditioned... Maybe a better initial guess for Y is needed to develop a better condition for matrix Z?

>> cond(Z)

ans =

     Inf

    >> svd(Z)

ans =

     0.104798737310912
                     0
                     0

[John D'Errico]

Sigh. My point is that the form you use for a linear least squares is a POOR one. backslash is already able to solve the problem directly as I wrote, which works for non-square systems.

The reason you use backslash in the form

%solve for dY
dY = Z\D;

is because backslash was written using high quality numerical methods to AVOID squaring the condition number of the problem. REREAD MY ANSWER.

That the condition number is infinite and SVD returns ZERO for some of the singular values means that your matrix is singular. You can NEVER solve that linear system, at least at that point.

It MAY reflect that you have chosen poor stating values. But as well, it begs the question, WHY ARE YOU WRITING YOUR OWN NONLINEAR ESTIMATION CODE? If you don't know how to effectively solve a simple linear system well, then why do you think you can write high quality nonlinear estimation code? (Sorry, but that is a completely valid question from my point of view.)

Use high quality tools written by professionals in the art rather than writing your own code from scratch.