%%Machine learning From Andrew
%%By youknowwho3_3 in 优快云 #GirlsHelpGirls#DOUBANEZU
%%Neural Networks Learning
%implement the backpropagation algorithm for neural networks and
%apply it to the task of hand-written digit recognition
%1.Neural Networks
%1.1 Visualizing the data
%1.2 Model representation
%1.3 Feedforward and cost function
%1.4 Regularized cost function
%2. Backpropagation
%2.1 Sigmoid gradient
%2.2 Random initialization
%2.3 Backpropagation
%2.4 Gradient checking
%2.5 Regularized neural networks
%2.6 Learning parameters using fmincg
%3. Visualizing the hidden layer
%3.1 Optional (ungraded) exercise
%%
%%%%%%%%%1.Neural Networks%%%%%%
load("ex4data1.mat");%size(X)=5000*400 size(y)=5000*1
m=size(X,1);%size(m)=5000
sel=randperm(size(X,1));%selis1*5000个随机数
sel=sel(1:100);%取sel的前100个
displayData(X(sel,:));%displayData provided
% Load the weights into variables Theta1 and Theta2
load('ex4weights.mat');%tTheta1 and Theta2
%size(Theta1)=25*401
%size(Theta2)=10*26
%%1.3Feedforward and costFunction
input_layer_size = 400; % 20x20 Input Images of Digits
hidden_layer_size = 25; % 25 hidden units
num_labels = 10; % 10 labels, from 1 to 10 (note that we have mapped "0" to label 10)
nn_params = [Theta1(:) ; Theta2(:)];%size(nn_params)= 100285*1 =25*401+10*26 %nn_params=[Theta1(:) ,Theta2(:)]; wrong
% Weight regularization parameter (we set this to 0 here).
lambda = 1;
J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
fprintf('Cost at parameters (loaded from ex4weights): %f', J);
% Call your sigmoidGradient function
sigmoidGradient(0);
%%%2.3Backpropagation
%详情见图片解说
%%%2.4Gradient Checking
checkNNGradients;%防止梯度爆炸
%%%2.5Regularized neural networks
% Check gradients by running checkNNGradients
lambda = 3;
checkNNGradients(lambda);
% Also output the costFunction debugging value
% This value should be about 0.576051
debug_J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
fprintf('Cost at (fixed) debugging parameters (w/lambda = 3): %f\n', debug_J);
%%%2.6Learning parameters using fmincg
options = optimset('MaxIter', 50);
lambda = 1;
% Create "short hand" for the cost function to be minimized
costFunction = @(p) nnCostFunction(p, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
% Now, costFunction is a function that takes in only one argument (the
% neural network parameters)
[nn_params, ~] = fmincg(costFunction, initial_nn_params, options);
% Obtain Theta1 and Theta2 back from nn_params
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), hidden_layer_size, (input_layer_size + 1));
Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), num_labels, (hidden_layer_size + 1));
pred = predict(Theta1, Theta2, X);
fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);
%%%3.Visualizing the hidden layer
% Visualize Weights
displayData(Theta1(:, 2:end));The cost function for neural networks with regularization is given by
Sigmoid gradient
%backpropagation反向传播图解
%%function nnCostFunction
%{
costFunction for Neural Networks with regularization
sigmoid gradient
formulat:
J=(1/m)*sum(i=1:m)sum(k=1:K)[-yk(i)*log(hkX(i)))-(1-yk(i)*log(1-hkX(i))]
+(lambda/2m)*{sum(j=1:25)sum(k=1:400)[thetajk(1)^2]+sum(j=1:10)sum(k=1:25)[thetajk(2)^2]}
g'(z)=d(g(z))/d(z)=g(z)*(1-g(z))
g(z)=1/(1+exp(-z))
%}
function [J grad] = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 2 layer neural network
%{
reshape:把A按照5*2reshape
% A = 1:10;
%B = reshape(A,[5,2])
%B = 5×2 = size(A)
% 1 6
% 2 7
% 3 8
% 4 9
% 5 10
%}
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), hidden_layer_size, (input_layer_size + 1));
%size(nn_params)=10285*1
%把nn_params(1:hidden_layer_size * (input_layer_size + 1))，按照hidden_layer_size*(input_layer_size + 1)reshape
%Theta1=nn_params(1:25*401),按照25*401reshape
%Theta1 is the original Theta1;不懂为啥要先变成nn_params
Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), num_labels, (hidden_layer_size + 1));
%Theta2=reshape(nn_params((1+25*401):end),10,25+1)
%size(Theta2)=10*26
% Setup some useful variables
m = size(X, 1);
% You need to return the following variables correctly
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));
%% Part 1 feedforword
a1=[ones(size(X,1),1),X];
%size(a1)=5000*401
z2=a1*Theta1'; %size(z2)=5000*25
a2=sigmoid(z2); %size(a2)=5000*25
a2=[ones(size(a2,1),1),a2]; %size(a2)=5000*26
z3=a2*Theta2';%size(z3)=5000*10
a3=sigmoid(z3);%size(a3)=5000*10
hThetaX=a3;%size(hThetaX)=5000*10
yVec = zeros(m,num_labels);%size(yVec)=5000*10;
for i = 1:m%m=5000
yVec(i,y(i)) = 1;%y变成了一个vec5000*10 y=几就在哪个位置上=1
end
%size(yVec)=5000*10
%J formulate: the costFunction for neural networks with regularizetion
J = 1/m * sum(sum(-1 * yVec .* log(hThetaX)-(1-yVec).* log(1-hThetaX)));
%size(yVec)=5000*10 size(hThetaX)=5000*10
%m=5000
%K=10
regularator = (sum(sum(Theta1(:,2:end).^2)) + sum(sum(Theta2(:,2:end).^2))) * (lambda/(2*m));
%size(Theta1)=25*401
%size(Theta2)=10*26
%刚好积分
J = J + regularator;
%% Part 2 Backpropagation
for t = 1:m
% For the input layer, where l=1:
a1 = [1; X(t,:)'];
%size(a1)=401*1
% For the hidden layers, where l=2:
z2 = Theta1 * a1; %size(z2)=25*1
a2 = [1; sigmoid(z2)]; %size(a2)=26*1
z3 = Theta2 * a2;%size(z3)=10*1
a3 = sigmoid(z3);%size(a3)=10*1
yy = ([1:num_labels]==y(t))';%size(yy)=10*1
%1到10里哪一个和y(t)相等，相应的位置等于1
% For the delta values:
delta_3 = a3 - yy;
%原理见图
delta_2 = (Theta2' * delta_3) .* [1; sigmoidGradient(z2)];
delta_2 = delta_2(2:end); % Taking of the bias row
% delta_1 is not calculated because we do not associate error with the input
% Big delta update
Theta1_grad = Theta1_grad + delta_2 * a1';%\Delta^{(l)} = \Delta^{(l)} + \delta^{(l+1)} (a^{(l)})^T
Theta2_grad = Theta2_grad + delta_3 * a2';
end
Theta1_grad = (1/m) * Theta1_grad + (lambda/m) * [zeros(size(Theta1, 1), 1) Theta1(:,2:end)];
Theta2_grad = (1/m) * Theta2_grad + (lambda/m) * [zeros(size(Theta2, 1), 1) Theta2(:,2:end)];
% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];
end
%%function sigmoidGradient
function g=sigmoidGradient(z)
%{
g'(z)=d(g(z))/d(z)=g(z)*(1-g(z))
g(z)=1/(1+exp(-z))
%}
g=zeros(size(z));
g=sigmoid(z).*(1-sigmoid(z));
end
function g=sigmoid(z)
g=zeros(size(z));
g=1./(1+exp(-z));
end
%%function displayData (provided)
function [h, display_array] = displayData(X, example_width)
%DISPLAYDATA Display 2D data in a nice grid
% [h, display_array] = DISPLAYDATA(X, example_width) displays 2D data
% stored in X in a nice grid. It returns the figure handle h and the
% displayed array if requested.
% Set example_width automatically if not passed in
if ~exist('example_width', 'var') || isempty(example_width)
example_width = round(sqrt(size(X, 2)));
end
% Gray Image
colormap(gray);
% Compute rows, cols
[m n] = size(X);
example_height = (n / example_width);
% Compute number of items to display
display_rows = floor(sqrt(m));
display_cols = ceil(m / display_rows);
% Between images padding
pad = 1;
% Setup blank display
display_array = - ones(pad + display_rows * (example_height + pad), ...
pad + display_cols * (example_width + pad));
% Copy each example into a patch on the display array
curr_ex = 1;
for j = 1:display_rows
for i = 1:display_cols
if curr_ex > m,
break;
end
% Copy the patch
% Get the max value of the patch
max_val = max(abs(X(curr_ex, :)));
display_array(pad + (j - 1) * (example_height + pad) + (1:example_height), ...
pad + (i - 1) * (example_width + pad) + (1:example_width)) = ...
reshape(X(curr_ex, :), example_height, example_width) / max_val;
curr_ex = curr_ex + 1;
end
if curr_ex > m,
break;
end
end
% Display Image
h = imagesc(display_array, [-1 1]);
% Do not show axis
axis image off
drawnow;
end
%function checkNNGradients
%provided
function checkNNGradients(lambda)
%CHECKNNGRADIENTS Creates a small neural network to check the
%backpropagation gradients
% CHECKNNGRADIENTS(lambda) Creates a small neural network to check the
% backpropagation gradients, it will output the analytical gradients
% produced by your backprop code and the numerical gradients (computed
% using computeNumericalGradient). These two gradient computations should
% result in very similar values.
%
if ~exist('lambda', 'var') || isempty(lambda)
lambda = 0;
end
input_layer_size = 3;
hidden_layer_size = 5;
num_labels = 3;
m = 5;
% We generate some 'random' test data
Theta1 = debugInitializeWeights(hidden_layer_size, input_layer_size);
Theta2 = debugInitializeWeights(num_labels, hidden_layer_size);
% Reusing debugInitializeWeights to generate X
X = debugInitializeWeights(m, input_layer_size - 1);
y = 1 + mod(1:m, num_labels)';
% Unroll parameters
nn_params = [Theta1(:) ; Theta2(:)];
% Short hand for cost function
costFunc = @(p)nnCostFunction(p, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
[cost, grad] = costFunc(nn_params);
numgrad = computeNumericalGradient(costFunc, nn_params);
% Visually examine the two gradient computations. The two columns
% you get should be very similar.
disp([numgrad grad]);
fprintf(['The above two columns you get should be very similar.\n' ...
'(Left-Your Numerical Gradient, Right-Analytical Gradient)\n\n']);
% Evaluate the norm of the difference between two solutions.
% If you have a correct implementation, and assuming you used EPSILON = 0.0001
% in computeNumericalGradient.m, then diff below should be less than 1e-9
diff = norm(numgrad-grad)/norm(numgrad+grad);
fprintf(['If your backpropagation implementation is correct, then \n' ...
'the relative difference will be small (less than 1e-9). \n' ...
'\nRelative Difference: %g\n'], diff);
end
%function fmincg
function [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5)
% Minimize a continuous differentialble multivariate function. Starting point
% is given by "X" (D by 1), and the function named in the string "f", must
% return a function value and a vector of partial derivatives. The Polack-
% Ribiere flavour of conjugate gradients is used to compute search directions,
% and a line search using quadratic and cubic polynomial approximations and the
% Wolfe-Powell stopping criteria is used together with the slope ratio method
% for guessing initial step sizes. Additionally a bunch of checks are made to
% make sure that exploration is taking place and that extrapolation will not
% be unboundedly large. The "length" gives the length of the run: if it is
% positive, it gives the maximum number of line searches, if negative its
% absolute gives the maximum allowed number of function evaluations. You can
% (optionally) give "length" a second component, which will indicate the
% reduction in function value to be expected in the first line-search (defaults
% to 1.0). The function returns when either its length is up, or if no further
% progress can be made (ie, we are at a minimum, or so close that due to
% numerical problems, we cannot get any closer). If the function terminates
% within a few iterations, it could be an indication that the function value
% and derivatives are not consistent (ie, there may be a bug in the
% implementation of your "f" function). The function returns the found
% solution "X", a vector of function values "fX" indicating the progress made
% and "i" the number of iterations (line searches or function evaluations,
% depending on the sign of "length") used.
%
% Usage: [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5)
%
% See also: checkgrad
%
% Copyright (C) 2001 and 2002 by Carl Edward Rasmussen. Date 2002-02-13
%
%
% (C) Copyright 1999, 2000 & 2001, Carl Edward Rasmussen
%
% Permission is granted for anyone to copy, use, or modify these
% programs and accompanying documents for purposes of research or
% education, provided this copyright notice is retained, and note is
% made of any changes that have been made.
%
% These programs and documents are distributed without any warranty,
% express or implied. As the programs were written for research
% purposes only, they have not been tested to the degree that would be
% advisable in any important application. All use of these programs is
% entirely at the user's own risk.
%
% [ml-class] Changes Made:
% 1) Function name and argument specifications
% 2) Output display
%
% Read options
if exist('options', 'var') && ~isempty(options) && isfield(options, 'MaxIter')
length = options.MaxIter;
else
length = 100;
end
RHO = 0.01; % a bunch of constants for line searches
SIG = 0.5; % RHO and SIG are the constants in the Wolfe-Powell conditions
INT = 0.1; % don't reevaluate within 0.1 of the limit of the current bracket
EXT = 3.0; % extrapolate maximum 3 times the current bracket
MAX = 20; % max 20 function evaluations per line search
RATIO = 100; % maximum allowed slope ratio
argstr = ['feval(f, X']; % compose string used to call function
for i = 1:(nargin - 3)
argstr = [argstr, ',P', int2str(i)];
end
argstr = [argstr, ')'];
if max(size(length)) == 2, red=length(2); length=length(1); else red=1; end
S=['Iteration '];
i = 0; % zero the run length counter
ls_failed = 0; % no previous line search has failed
fX = [];
[f1 df1] = eval(argstr); % get function value and gradient
i = i + (length<0); % count epochs?!
s = -df1; % search direction is steepest
d1 = -s'*s; % this is the slope
z1 = red/(1-d1); % initial step is red/(|s|+1)
while i < abs(length) % while not finished
i = i + (length>0); % count iterations?!
X0 = X; f0 = f1; df0 = df1; % make a copy of current values
X = X + z1*s; % begin line search
[f2 df2] = eval(argstr);
i = i + (length<0); % count epochs?!
d2 = df2'*s;
f3 = f1; d3 = d1; z3 = -z1; % initialize point 3 equal to point 1
if length>0, M = MAX; else M = min(MAX, -length-i); end
success = 0; limit = -1; % initialize quanteties
while 1
while ((f2 > f1+z1*RHO*d1) || (d2 > -SIG*d1)) && (M > 0)
limit = z1; % tighten the bracket
if f2 > f1
z2 = z3 - (0.5*d3*z3*z3)/(d3*z3+f2-f3); % quadratic fit
else
A = 6*(f2-f3)/z3+3*(d2+d3); % cubic fit
B = 3*(f3-f2)-z3*(d3+2*d2);
z2 = (sqrt(B*B-A*d2*z3*z3)-B)/A; % numerical error possible - ok!
end
if isnan(z2) || isinf(z2)
z2 = z3/2; % if we had a numerical problem then bisect
end
z2 = max(min(z2, INT*z3),(1-INT)*z3); % don't accept too close to limits
z1 = z1 + z2; % update the step
X = X + z2*s;
[f2 df2] = eval(argstr);
M = M - 1; i = i + (length<0); % count epochs?!
d2 = df2'*s;
z3 = z3-z2; % z3 is now relative to the location of z2
end
if f2 > f1+z1*RHO*d1 || d2 > -SIG*d1
break; % this is a failure
elseif d2 > SIG*d1
success = 1; break; % success
elseif M == 0
break; % failure
end
A = 6*(f2-f3)/z3+3*(d2+d3); % make cubic extrapolation
B = 3*(f3-f2)-z3*(d3+2*d2);
z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3)); % num. error possible - ok!
if ~isreal(z2) || isnan(z2) || isinf(z2) || z2 < 0 % num prob or wrong sign?
if limit < -0.5 % if we have no upper limit
z2 = z1 * (EXT-1); % the extrapolate the maximum amount
else
z2 = (limit-z1)/2; % otherwise bisect
end
elseif (limit > -0.5) && (z2+z1 > limit) % extraplation beyond max?
z2 = (limit-z1)/2; % bisect
elseif (limit < -0.5) && (z2+z1 > z1*EXT) % extrapolation beyond limit
z2 = z1*(EXT-1.0); % set to extrapolation limit
elseif z2 < -z3*INT
z2 = -z3*INT;
elseif (limit > -0.5) && (z2 < (limit-z1)*(1.0-INT)) % too close to limit?
z2 = (limit-z1)*(1.0-INT);
end
f3 = f2; d3 = d2; z3 = -z2; % set point 3 equal to point 2
z1 = z1 + z2; X = X + z2*s; % update current estimates
[f2 df2] = eval(argstr);
M = M - 1; i = i + (length<0); % count epochs?!
d2 = df2'*s;
end % end of line search
if success % if line search succeeded
f1 = f2; fX = [fX' f1]';
fprintf('%s %4i | Cost: %4.6e\r', S, i, f1);
s = (df2'*df2-df1'*df2)/(df1'*df1)*s - df2; % Polack-Ribiere direction
tmp = df1; df1 = df2; df2 = tmp; % swap derivatives
d2 = df1'*s;
if d2 > 0 % new slope must be negative
s = -df1; % otherwise use steepest direction
d2 = -s'*s;
end
z1 = z1 * min(RATIO, d1/(d2-realmin)); % slope ratio but max RATIO
d1 = d2;
ls_failed = 0; % this line search did not fail
else
X = X0; f1 = f0; df1 = df0; % restore point from before failed line search
if ls_failed || i > abs(length) % line search failed twice in a row
break; % or we ran out of time, so we give up
end
tmp = df1; df1 = df2; df2 = tmp; % swap derivatives
s = -df1; % try steepest
d1 = -s'*s;
z1 = 1/(1-d1);
ls_failed = 1; % this line search failed
end
if exist('OCTAVE_VERSION')
fflush(stdout);
end
end
fprintf('\n');
end
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
