From: <Saved by Microsoft Internet Explorer 5>
Date: Mon, 30 Dec 2002 11:01:09 -0800

function [w1, w2, w3] = backprog(I, J, K, L, ip, target)
%
% 4-layer (two hidden layers) backprogation training program
% the four layers are indexed by I (output), J, K, (hidden), L (input)
% ip: input feature vectors of dimension l x patterns
% target: output vectors of dimension i x patterns
% all other parameters are coded inside the backprog
%

% verify data correctness
[row, patterns] = size(ip);
if (row ~= L)
fprintf('input row dimension mut be: ', L);
return;
end
[row, patterns1] = size(target);
if (row ~= I)
fprintf('target row dimension must be: ', I);
return;
end
if (patterns1 ~= patterns)
fprintf('target col dimension must be: ', patterns);
return;
end

% initialize parameters
steep = 1; % steepness in tanh (beta x)
learning = 0.1; % learning rate
momentum = 0.9; % momentum
residue = 0.02; % maximum residue error
epochs = 1000; % maximum number of iterations
wKL = rand(K,L+1); % +1 is for bias
wJK = rand(J,K+1);
wIJ = rand(I,J+1);
dwKLold = zeros(K,L+1);% +1 is for bias
dwJKold = zeros(J,K+1);
dwIJold = zeros(I,J+1);
dwKLnew = zeros(K,L+1);% +1 is for bias
dwJKnew = zeros(J,K+1);
dwIJnew = zeros(I,J+1);

for (iter=1; iter<=epochs; iter++) % iterations or epochs
for (ptn=1; ptn<=patterns; ptn++) % patterns
aLp = ip(:,ptn); % ptn-th data
desired = target(:,ptn); % desired output
bKp = wKL*[aLp; -1];
aKp = tanh(bKp);
bJp = wJK*[aKp; -1];
aJp = tanh(bJp);
bIp = wIJ*[aJp; -1];
aIp = tanh(bIp);

% at the output layer
sigma = steep * (1-tanh(bIp).^2) .* (desired-aIp);
dwIJnew = learning * sigma * [aJp; -1]';

% push back one layer
sigma = steep * (1-tanh(bJp).^2) .* (wIJ' * sigma);
dwJKnew = learning * sigma * [aKp; -1]';

% push back one more layer
sigma = steep * (1-tanh(bKp).^2) .* (wJK' * sigma);
dwKLnew = learning * sigma * [aLp, -1]';

% update weight (incrementally) for each sample presented
wIJ = WIJ + dwIJnew + momentum * dwIJold;
wJK = WJK + dwJKnew + momentum * dwJKold;
wKL = WKL + dwKLnew + momentum * dwKLold;
dwIJold = dwIJnew; % record old delta for momentum
dwJKold = dwJKnew;
dwKLold = dwKLnew;
end
output = tanh(wIJ*tanh(wJK*tanh(wKL*ip)));
error = 0;
for (col=1; col<=patterns; col++)
error = error + norm(target(:,col)-output(:,col));
end
error = error / patterns;
fprintf ('after: ', iter, 'iteration, the error is: ', error);
if (error <= residue)
break;
end
end