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

[Sergei Steshenko]

Hello,

running my calculations based on optim-1.0.4/leasqr.m I've noticed
occasional "division by zero" message.

Looking into the code I found that division by zero occurs at lines ##
301, 302 in the following piece of code:

    281 % CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
    282 % re-evaluate the Jacobian at optimal values
    283 jac=feval(dFdp,x,f,p,dp,F);
    284 msk = dp ~= 0;
    285 n = sum(msk);           % reduce n to equal number of estimated parameters
    286 jac = jac(:, msk);      % use only fitted parameters
    287
    288 %% following section is Ray Muzic's estimate for covariance and correlation
    289 %% assuming covariance of data is a diagonal matrix proportional to
    290 %% diag(1/wt.^2).
    291 %% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1
    292
    293 if exist('sparse')  % save memory
    294   Q=sparse(1:m,1:m,1./wt.^2);
    295   Qinv=sparse(1:m,1:m,wt.^2);
    296 else
    297   Q=diag((0*wt+1)./(wt.^2));
    298   Qinv=diag(wt.*wt);
    299 end
    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
    306 d=sqrt(abs(diag(covp)));
    307 corp=covp./(d*d');
    308
    309 if exist('sparse')
    310   covr=spdiags(covr,0);
    311   stdresid=resid./sqrt(spdiags(Vy,0));
    312 else
    313   covr=diag(covr);                 % convert returned values to compact storage
    314   stdresid=resid./sqrt(diag(Vy));  % compute then convert for compact storage
    315 end
.

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.

I have worked around the problem by changing number of observed points in
my code, so the division by zero does not occur, but I think somebody who
in details understands the algorithm should properly fix 'leasqr'.

Thanks,
  Sergei.