"occasional" division by zero in optim-1.0.4/leasqr.m

[Sergei Steshenko]

--- On Fri, 7/10/09, Rob Mahurin <rob at utk.edu> wrote:

> From: Rob Mahurin <rob at utk.edu>
> Subject: Re: "occasional" division by zero in optim-1.0.4/leasqr.m
> To: "Sergei Steshenko" <sergstesh at yahoo.com>
> Cc: bug-octave at octave.org
> Date: Friday, July 10, 2009, 12:47 PM
> On Jul 10, 2009, at 2:07 PM, Sergei
> Steshenko wrote:
> > Looking into the code I found that division by zero
> occurs at lines ##
> > 301, 302 in the following piece of code:
> > [...]
> >     300 resid=y-f;   
>                
>                
> %un-weighted residuals
> >     301
> covr=resid'*Qinv*resid*Q/(m-n);       
>          %covariance of
> residuals
> >     302 Vy=1/(1-n/m)*covr;  %
> Eq. 7-13-22, Bard     
>    %covariance of the data
> >     303
> >     304
> jtgjinv=inv(jac'*Qinv*jac);       
>          
>    %argument of inv may be singular
> >     305
> covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13,
> Bard %cov of parm est
> > 
> > I have verified that division by zero occurs because n
> == m; in my case
> > they both were 10, but they can be anything, of
> course.
> > 
> > Earlier in the code m, n are set this way (line #161)
> > 
> >     159 y=y(:); wt=wt(:);
> pin=pin(:); dp=dp(:); %change all vectors to columns
> >     160 % check data vectors- same
> length?
> >     161 m=length(y);
> n=length(pin); p=pin;[m1,m2]=size(x);
> >     162 if m1~=m
> ,error('input(x)/output(y) data must have same number of
> rows ') ,end;
> > .
> > 
> > In English 'm' is number of elements in 'y' vector
> which is the vector
> > of observed function values, and 'n' is number of
> elements in parameters
> > ('p') which are supposed to be approximated by
> 'leasqr'.
> > 
> > Though I know nothing about covariance matrix (the
> piece of code where
> > division by zero occurs) and a little about
> Levenberg-Marquardt algorithm
> > implemented by 'leasqr', I do not think there should
> be a hard requirement
> > for number of elements in 'y' and 'p' to be
> different.
> 
> If there are as many observations as parameters, it's
> possible in principle to find a set of parameters so that
> all the residuals are zero.  You must have more
> observations than parameters to get error estimates on the
> parameters, which is where all the covariance matrices come
> in.
> 
> What does leasqr do if you feed it, say, three observations
> and four parameters?  I think in that case it should
> just fail ("underdetermined problem" or something).  I
> don't know whether the should try to fit for the case m==n
> or not.
> 
> Cheers,
> Rob
> 
> --Rob Mahurin
> University of Manitoba, Dept. of Physics & Astronomy
> at:    Oak Ridge National
> Laboratory    865 207 2594
>     Oak Ridge, Tennessee   
>     rob at utk.edu


I think the algorithm should try to fit - n == m is "exactly" determined
problem.

I think that the algorithm should try to fit in any case - because it fits
in the least squares sense.

The code which causes division by zero is _after_ the fitting part, i.e.
it is not necessary for the algorithm to work. In my calculations I simply
ignore results produced by the sometimes failing piece of code - not 
because it sometimes fails, but because I simply do no need them.

Thanks,
  Sergei.