Documentation of gn

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[g,R,grad]=gn(z,th,lim,maxsize,index)

Help text

GN     Computes the Gauss-Newton direction for PE-criteria

   [G,R] = gn(Z,TH,LIM,MAXSIZE,INDEX)

   G : The Gauss-Newton direction
   R : The Hessian of the criterion (using the GN-approximation)

   The routine is intended primarily as a subroutine to MINCRIT.
   See PEM for an explanation of arguments.

Cross-Reference Information

Listing of function gn

function [g,R,grad]=gn(z,th,lim,maxsize,index)

%   L. Ljung 10-1-86
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.3 $  $Date: 1997/12/02 03:43:35 $

% *** Set up the model orders ***

nu=th(1,3); na=th(1,4); nb=th(1,5:4+nu); nc=th(1,5+nu);
nd=th(1,6+nu); nf=th(1,7+nu:6+2*nu); nk=th(1,7+2*nu:6+3*nu);
[Ncap,nz]=size(z); n=na+sum(nb)+nc+nd+sum(nf);
nmax=max([na nb+nk-ones(1,nu) nc nd nf]);

% *** Prepare for gradients calculations ***

[e,v,w,a,b,c,d,f]=pe(z,th);
[me,ne]=size(e);
ll=lim*ones(me,ne);
if lim==0,el=e;else la=abs(e)+eps*ll;el=e.*(min(la,ll)./la);end
%el=min(e,ll); el=max(el,-ll); clear ll;
if sum(nf)==0 , clear w, end
if na>0, yf=filter(-d,c,z(:,1));end
if nc>0, ef=filter(1,c,e); end
if nd>0, vf=filter([-1],c,v); end
for k=1:nu
      gg=conv(c,f(k,:));
      uf(:,k)=filter(d,gg,z(:,k+1));
      if nf(k)>0, wf(:,k)=filter(-d,gg,w(:,k));end
end

% *** Compute the gradient PSI. If N>M do it in portions ***

M=floor(maxsize/n);n1=length(index);
R=zeros(n1);Fcap=zeros(n1,1);
for k=nmax:M:Ncap-1
      jj=(k+1:min(Ncap,k+M));
      psi=zeros(length(jj),n);
      for kl=1:na, psi(:,kl)=yf(jj-kl);,end
      ss=na;ss1=na+sum(nb)+nc+nd;
      for ku=1:nu
             for kl=1:nb(ku), psi(:,ss+kl)=uf(jj-kl-nk(ku)+1,ku);end
             for kl=1:nf(ku), psi(:,ss1+kl)=wf(jj-kl,ku);end
             ss=ss+nb(ku);ss1=ss1+nf(ku);
      end
      for kl=1:nc, psi(:,ss+kl)=ef(jj-kl);end
      ss=ss+nc;
      for kl=1:nd, psi(:,ss+kl)=vf(jj-kl);,end
      psi=psi(:,index);
      R=R+psi'*psi;  Fcap=Fcap+ psi'*el(jj);
end

% *** Compute the Gauss-Newton direction g ***

if Ncap>M, g=R\Fcap; else g=psi\el(jj);end
g1=zeros(n,1);g1(index)=g;g=g1;grad=zeros(n,1);grad(index)=Fcap;