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

[Rob Mahurin]

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