home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / LEASTSQ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  8.3 KB  |  281 lines
  1. function [x,OPTIONS,f,JOCOB] = leastsq(FUN,x,OPTIONS,GRADFUN,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
  2. %LEASTSQ Solves non-linear least squares problems.
  3. %     LEASTSQ solves problems of the form:
  4. %    min  sum {FUN(X).^2}    where FUN and X may be vectors of matrices.   
  5. %             x
  6. %
  7. %    X=LEASTSQ('FUN',X0) starts at the matrix X0 and finds a minimum to the
  8. %    sum of squares of the functions described in FUN. FUN is usually
  9. %    an M-file which returns a matrix of objective functions: F=FUN(X).
  10. %
  11. %    X=LEASTSQ('FUN',X0,OPTIONS) allows a vector of optional parameters to
  12. %    be defined. OPTIONS(2) is a measure of the precision required for the 
  13. %    values of X at the solution. OPTIONS(3) is a measure of the precision
  14. %    required of the objective function at the solution. See HELP FOPTIONS. 
  15. %
  16. %    X=LEASTSQ('FUN',X0,OPTIONS,'GRADFUN') enables a function'GRADFUN'
  17. %    to be entered which returns the partial derivatives of the functions,
  18. %    dF/dX, (stored in columns) at the point X: gf = GRADFUN(X).
  19. %    
  20.  
  21. %    [X,OPTIONS,F,J]=LEASTSQ('FUN',X0,...) returns, F, the residuals and
  22. %    J the Jacobian of the function FUN at the solution. 
  23. %    F and J can be used with to calculate the variance and
  24. %    confidence intervals using CONFINT(X,S,J).
  25.  
  26. %    Copyright (c) 1990 by the MathWorks, Inc.
  27. %    Andy Grace 7-9-90.
  28.  
  29. %    The default algorithm is the Levenberg-Marquardt method with a 
  30. %    mixed quadratic and cubic line search procedure.  A Gauss-Newton
  31. %    method is selected by setting  OPTIONS(5)=1. 
  32. %
  33. %    X=LEASTSQ('FUN',X,OPTIONS,GRADFUN,P1,P2,..) allows
  34. %    coefficients, P1, P2, ... to be passed directly to FUN:
  35. %    F=FUN(X,P1,P2,...). Empty arguments ([]) are ignored.
  36.  
  37. % ------------Initialization----------------
  38. XOUT=x(:);
  39. [nvars]=length(XOUT);
  40.  
  41. evalstr = [FUN];
  42. if ~any(FUN<48)
  43.     evalstr=[evalstr, '(x'];
  44.     for i=1:nargin - 4
  45.         evalstr = [evalstr,',P',int2str(i)];
  46.     end
  47.     evalstr = [evalstr, ')'];
  48. end
  49.  
  50.  
  51. if nargin < 3, OPTIONS=[]; end
  52. if nargin < 4, GRADFUN=[]; end
  53.  
  54. if length(GRADFUN)
  55.     evalstr2 = [GRADFUN];
  56.     if ~any(GRADFUN<48)
  57.         evalstr2 = [evalstr2, '(x'];
  58.         for i=1:nargin - 4
  59.             evalstr2 = [evalstr2,',P',int2str(i)];
  60.         end
  61.         evalstr2 = [evalstr2, ')'];
  62.     end
  63. end
  64.  
  65. f = eval(evalstr);
  66. f=f(:);
  67. nfun=length(f);
  68. GRAD=zeros(length(XOUT),nfun);
  69. OLDX=XOUT;
  70. MATX=zeros(3,1);
  71. MATL=[f'*f;0;0];
  72. OLDF=f'*f;
  73. FIRSTF=f'*f;
  74. [OLDX,OLDF,OPTIONS]=lsint(XOUT,f,OPTIONS);
  75. PCNT = 0;
  76. EstSum=0.5;
  77. GradFactor=0; 
  78. CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);
  79.  
  80. OPTIONS(10)=1;
  81. status=-1;
  82.  
  83. while status~=1
  84.  
  85. % Work Out Gradients
  86.     if ~length(GRADFUN) | OPTIONS(9)
  87.         OLDF=f;
  88.         CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
  89.         for gcnt=1:nvars
  90.             temp = XOUT(gcnt);
  91.             XOUT(gcnt) = temp +CHG(gcnt);
  92.             x(:) = XOUT;
  93.             f(:) = eval(evalstr);
  94.             GRAD(gcnt,:)=(f-OLDF)'/(CHG(gcnt));
  95.             XOUT(gcnt) = temp;
  96.         end
  97.         f = OLDF;
  98.         OPTIONS(10)=OPTIONS(10)+nvars;
  99. % Gradient check
  100.         if OPTIONS(9) == 1
  101.             GRADFD = GRAD;
  102.             x(:)=XOUT; GRAD = eval(evalstr2);
  103.             graderr(GRADFD, GRAD, evalstr2);
  104.             OPTIONS(9) = 0;
  105.         end
  106.     else
  107.         x(:) = XOUT;
  108.         OPTIONS(11)=OPTIONS(11)+1;
  109.         GRAD = eval(evalstr2);
  110.     end
  111.     % Try to set difference to 1e-8 for next iteration
  112.     if nfun==1
  113.         CHG = nfun*1e-8./GRAD;
  114.     else
  115.         CHG = nfun*1e-8./sum(abs(GRAD)')';
  116.     end
  117.  
  118.     gradf = 2*GRAD*f;
  119.     fnew = f'*f;
  120. %---------------Initialization of Search Direction------------------
  121.     if status==-1
  122.         if cond(GRAD)>1e8
  123.             SD=-(GRAD*GRAD'+(norm(GRAD)+1)*(eye(nvars,nvars)))\(GRAD*f);
  124.             if OPTIONS(5)==0, GradFactor=norm(GRAD)+1; end
  125.             how='COND';
  126.         else
  127. %        SD=GRAD'\(GRAD'*X-F)-X;
  128.             SD=-(GRAD*GRAD'+GradFactor*(eye(nvars,nvars)))\(GRAD*f);
  129.         end
  130.         FIRSTF=fnew;
  131.         OLDG=GRAD;
  132.         GDOLD=gradf'*SD;
  133.         % OPTIONS(18) controls the initial starting step-size.
  134.         % If OPTIONS(18) has been set externally then it will
  135.         % be non-zero, otherwise set to 1.
  136.         if OPTIONS(18) == 0, OPTIONS(18)=1; end
  137.         if OPTIONS(1)>0
  138.             disp([sprintf('%5.0f %12.6g %12.3g ',OPTIONS(10),fnew,OPTIONS(18)),sprintf('%12.3g  ',GDOLD)]);
  139.         end
  140.         XOUT=XOUT+OPTIONS(18)*SD;
  141.         if OPTIONS(5)==0
  142.             newf=GRAD'*SD+f;
  143.             GradFactor=newf'*newf;
  144.             SD=-(GRAD*GRAD'+GradFactor*(eye(nvars,nvars)))\(GRAD*f); 
  145.         end
  146.         newf=GRAD'*SD+f;
  147.         XOUT=XOUT+OPTIONS(18)*SD;
  148.         EstSum=newf'*newf;
  149.         status=0;
  150.         if OPTIONS(7)==0; PCNT=1; end
  151.         
  152.     else
  153. %-------------Direction Update------------------
  154.         gdnew=gradf'*SD;
  155.         if OPTIONS(1)>0, 
  156.             num=[sprintf('%5.0f %12.6g %12.3g ',OPTIONS(10),fnew,OPTIONS(18)),sprintf('%12.3g  ',gdnew)];
  157.         end
  158.         if gdnew>0 & fnew>FIRSTF
  159.  
  160. % Case 1: New function is bigger than last and gradient w.r.t. SD -ve
  161. % ... interpolate. 
  162.             how='inter';
  163.             [stepsize]=cubici1(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
  164.             OPTIONS(18)=0.9*stepsize;
  165.         elseif fnew<FIRSTF
  166.  
  167. %  New function less than old fun. and OK for updating 
  168. %         .... update and calculate new direction. 
  169.             [newstep,fbest] =cubici3(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
  170.             if fbest>fnew,fbest=0.9*fnew; end 
  171.             if gdnew<0
  172.                 how='incstep';
  173.                 if newstep<OPTIONS(18),  newstep=(2*OPTIONS(18)+1e-4); how=[how,'IF']; end
  174.                 OPTIONS(18)=abs(newstep);
  175.             else 
  176.                 if OPTIONS(18)>0.9
  177.                     how='int_step';
  178.                     OPTIONS(18)=min([1,abs(newstep)]);
  179.                 end
  180.             end
  181. % SET DIRECTION.
  182. % Gauss-Newton Method    
  183.             temp=1;
  184.             if OPTIONS(5)==1 
  185.                 if OPTIONS(18)>1e-8 & cond(GRAD)<1e8
  186.                     SD=GRAD'\(GRAD'*XOUT-f)-XOUT;
  187.                     if SD'*gradf>eps,how='ERROR- GN not descent direction',  end
  188.                     temp=0;
  189.                 else
  190.                     if OPTIONS(1) > 0
  191.                         disp('Conditioning of Gradient Poor - Switching To LM method')
  192.                     end
  193.                     how='CHG2LM';
  194.                     OPTIONS(5)=0;
  195.                     OPTIONS(18)=abs(OPTIONS(18));                
  196.                 end
  197.             end
  198.             
  199.             if (temp)      
  200. % Levenberg_marquardt Method N.B. EstSum is the estimated sum of squares.
  201. %                                 GradFactor is the value of lambda.
  202. % Estimated Residual:
  203.                 if EstSum>fbest
  204.                     GradFactor=GradFactor/(1+OPTIONS(18));
  205.                 else
  206.                     GradFactor=GradFactor+(fbest-EstSum)/(OPTIONS(18)+eps);
  207.                 end
  208.                 SD=-(GRAD*GRAD'+GradFactor*(eye(nvars,nvars)))\(GRAD*f); 
  209.                 OPTIONS(18)=1; 
  210.                 estf=GRAD'*SD+f;
  211.                 EstSum=estf'*estf;
  212.                 if OPTIONS(1)>0, num=[num,sprintf('%12.6g ',GradFactor)]; end
  213.             end
  214.             gdnew=gradf'*SD;
  215.  
  216.             OLDX=XOUT;
  217. % Save Variables
  218.             FIRSTF=fnew;
  219.             OLDG=gradf;
  220.             GDOLD=gdnew;    
  221.  
  222.             % If quadratic interpolation set PCNT
  223.             if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1);  MATL(1)=fnew; end
  224.         else 
  225. % Halve Step-length
  226.             how='Red_Step';
  227.             if fnew==FIRSTF,
  228.                 if OPTIONS(1)>0,
  229.                     disp('No improvement in search direction: Terminating')
  230.                 end
  231.                 status=1;
  232.             else
  233.                 OPTIONS(18)=OPTIONS(18)/8;
  234.                 if OPTIONS(18)<1e-8
  235.                     OPTIONS(18)=-OPTIONS(18);
  236.                 end
  237.             end
  238.         end
  239.         XOUT=OLDX+OPTIONS(18)*SD;
  240.         if OPTIONS(1)>0,disp([num,how]),end
  241.  
  242.     end %----------End of Direction Update-------------------
  243.     if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1);  MATL(1)=fnew; end
  244. % Check Termination 
  245.     if max(abs(SD))< OPTIONS(2) & (gradf'*SD) < OPTIONS(3) & max(abs(gradf)) < 10*(OPTIONS(3)+OPTIONS(2))
  246.         if OPTIONS(1) > 0
  247.             disp('Optimization Terminated Successfully')  
  248.         end
  249.         status=1; 
  250.     elseif OPTIONS(10)>OPTIONS(14)
  251.         disp('maximum number of iterations has been exceeded');
  252.         if OPTIONS(1)>0
  253.             disp('Increase OPTIONS(14)')
  254.         end
  255.         status=1;
  256.     else
  257.  
  258. % Line search using mixed polynomial interpolation and extrapolation.
  259.         if PCNT~=0
  260.             while PCNT > 0
  261.                 x(:) = XOUT; f(:)  = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
  262.                 fnew = f'*f;
  263.             % <= used in case when no improvement found.
  264.                 if fnew <= OLDF'*OLDF, OX = XOUT; OLDF=f; end
  265.                 [PCNT,MATL,MATX,steplen,fnew,how]=searchq(PCNT,fnew,OLDX,MATL,MATX,SD,GDOLD,OPTIONS(18),how);
  266.                 OPTIONS(18)=steplen;
  267.                 XOUT=OLDX+steplen*SD;
  268.                 if fnew==FIRSTF,  PCNT=0; end
  269.             end
  270.             XOUT = OX;
  271.             f=OLDF;
  272.         else
  273.             x(:)=XOUT; f(:) = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
  274.         end
  275.     end
  276. end
  277. OPTIONS(8) = fnew;
  278. XOUT=OLDX;
  279. x(:)=XOUT;
  280.  
  281.