aprod3aug.m

书籍“Regularization tools for training large feed-forward neural networks using Automatic Differentiat

function [Af] =aprod3aug(mode,m,n,f,mu,rw, ...
	  w1,b1,Output_Layer2,Deriv_Layer2,df1,w2,b2,Output_Layer3,...
	  Deriv_Layer3,df2,w3,b3,Output_Layer4,Deriv_Layer4,df3,P,T)
%Call: [Af] =aprod3aug(mode,m,n,f,mu,rw, ...
%	  w1,b1,Output_Layer2,Deriv_Layer2,df1,w2,b2,Output_Layer3,Deriv_Layer3,
%	  ..., df2,w3,b3,Output_Layer4,Deriv_Layer4,df3,P,T)
%This routine computes Jaug*f or Jaug'*f, where J is the Jacobian corresponding
%to a FFANN with two hidden layers.
%Jaug = ( J  )
%       (mu*I)
%It should normally be used with the routine truncLS3 which is a
%modification of LSQR written by Paige and Saunders
%
% If mode == 1 aprod3aug returns y=Jaug*f
% If mode == 2 aprod3aug returns y=Jaug'*f
%
[Fst,q]=size(P);	% # nodes in Input layer (first layer)
[Snd,r]=size(w1);	% # nodes in 1'st hidden layer (second layer)
if Fst ~= r		% q is # samples presented to the inputs
   error('In aprod3aug: Fst .ne. r');
end
[Thrd,r]=size(w2);	% # nodes in 2'nd hidden layer (third layer)
if Snd ~= r
   error('In aprod3aug: Snd .ne. r')
end
[Fth,r]=size(w3);	% # nodes in output layer (forth layer)

Om=Fst*Snd+Snd;		% # weights in first hidden layer
Nm=Om+Snd*Thrd+Thrd;	% # weights up to and including second hidden layer
Mm=Nm+Thrd*Fth+Fth;	% # total weights
lenf=max(size(f));
if mode == 1
  if Mm ~= lenf
     error('In aprod3aug: Mm .ne. lenf (mode == 1)');
  else
     Af=zeros(q*Fth+Mm,1);
  end
elseif (q*Fth+Mm) ~= lenf
      error('In aprod3aug: q*Fth+Mm .ne. lenf (mode == 2)');
    else
      Af=zeros(Mm,1);
end
%
if mode == 1
%Compute Jaug*f by using forward mode
%Unpack the vector elements such that same kind of computation as
%in the function evaluation can be done.

%First hidden layer
TEMP=ones(1,q);
w1=reshape(f(1:Fst*Snd),Snd,Fst);
b1=f(Fst*Snd+1:Fst*Snd+Snd);
O_L2=b1*TEMP+w1*P;
S2=Deriv_Layer2 .* O_L2;

%Second hidden layer
w2temp=reshape(f(Om+1:Om+Snd*Thrd),Thrd,Snd);
b2=f(Om+Snd*Thrd+1:Om+Snd*Thrd+Thrd);
O_L3=b2*TEMP+w2*S2+w2temp*Output_Layer2;
S3=Deriv_Layer3 .* O_L3;

%Output layer
w3temp=reshape(f(Nm+1:Nm+Fth*Thrd),Fth,Thrd);
b3=f(Nm+Fth*Thrd+1:Mm);
O_L4=b3*TEMP+w3*S3+w3temp*Output_Layer3;
S4=Deriv_Layer4 .* O_L4;
Af=[S4(:); mu*f];

else % Now mode should be 2
   %Compute Jaug'*f
   e=reshape(f(1:Fth*q),Fth,q);
   d3=feval(df3,Output_Layer4,e);
   d2=feval(df2,Output_Layer3,d3,w3);
   d1=feval(df1,Output_Layer2,d2,w2);
   [dw1,db1]=learnbp(P,d1,1);
   [dw2,db2]=learnbp(Output_Layer2,d2,1);
   [dw3,db3]=learnbp(Output_Layer3,d3,1);
   Af=[dw1(:); db1(:); dw2(:); db2(:); dw3(:); db3(:)];
   Af=Af+mu*f(q*Fth+1:q*Fth+Mm);
end %of if mode == 1
end