实验代码如下 lab3.m：%% Initialization clear ; close all...

实验代码如下
lab3.m：%% Initialization
clear ;
close all;
clc

%% Setup the parameters you will use for this lab
% 20x20 Input Images of Digits
input_layer_size = 400;
% 25 hidden units
hidden_layer_size = 25;
% 10 labels, from 1 to 10. "0" is mapped to label 10)
num_labels = 10;

%% =========== Part 1: Loading and Visualizing Data =============
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('lab3data.mat');
m = size(X, 1);

% Randomly select 100 data points to display
sel = randperm(size(X, 1));
sel = sel(1:100);

datademonstrate(X(sel, :));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 2: Loading Parameters ================
%

fprintf('\nLoading Saved Neural Network Parameters ...\n')

% Load the weights into variables Theta1 and Theta2
load('lab3weights.mat');

% Unroll parameters
nn_params = [Theta1(:) ; Theta2(:)];

%% ================ Part 3: Compute Cost (Feedforward) ================
% first start by implementing the feedforward part of the neural network
% that returns the cost only.Complete the code in nncfunction.m to return
% cost.
%
% first implementing the feedforward cost without regularization

%
fprintf('\nFeedforward Using Neural Network ...\n')

% Weight regularization parameter (we set this to 0 here).
lambda = 0;

J = nncfunction(nn_params, input_layer_size, hidden_layer_size, ...
num_labels, X, y, lambda);

fprintf(['Cost at parameters (loaded from lab3weights): %f '...
'\n(this value should be about 0.287629)\n'], J);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% =============== Part 4: Implement Regularization ===============
% implementing the regularization with the cost.
%

fprintf('\nChecking Cost Function (w/ Regularization) ... \n')

% Weight regularization parameter (we set this to 1 here).
lambda = 1;

J = nncfunction(nn_params, input_layer_size, hidden_layer_size, ...
num_labels, X, y, lambda);

fprintf(['Cost at parameters (loaded from lab3weights): %f '...
'\n(this value should be about 0.383770)\n'], J);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 5: Sigmoid Gradient ================
%

fprintf('\nEvaluating sigmoid gradient...\n')

g = sgmdgrad([1 -0.5 0 0.5 1]);
fprintf('Sigmoid gradient evaluated at [1 -0.5 0 0.5 1]:\n ');
fprintf('%f ', g);
fprintf('\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 6: Initializing Pameters ================
% starting to implment a two layer neural network that classifies digits.
%
fprintf('\nInitializing Neural Network Parameters ...\n')

initial_Theta1 = rdinitweights(input_layer_size, hidden_layer_size);
initial_Theta2 = rdinitweights(hidden_layer_size, num_labels);

% Unroll parameters
initial_nn_params = [initial_Theta1(:) ; initial_Theta2(:)];

%% =============== Part 7: Implement Backpropagation ===============
% investigate coresponding codes snippet in nncfunction.m
%
fprintf('\nChecking Backpropagation... \n');

% Check gradients
gradcheck;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% =============== Part 8: Implement Regularization ===============
% implementing the regularization with the cost and gradient.
%

fprintf('\nChecking Backpropagation (w/ Regularization) ... \n')

% Check gradients
lambda = 3;
gradcheck(lambda);

% output the costfunction debugging values
debug_J = nncfunction(nn_params, input_layer_size, ...
hidden_layer_size, num_labels, X, y, lambda);

fprintf(['\n\nCost at (fixed) debugging parameters (w/ lambda = 10): %f ' ...
'\n(this value should be about 0.576051)\n\n'], debug_J);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 8: Training NN ===================
% use "fmincg", which works similarly to "fminunc".
%
fprintf('\nTraining Neural Network... \n')

options = optimset('MaxIter', 50);

% You should also try different values of lambda
lambda = 1;

% Create "short hand" for the cost function to be minimized
costFunction = @(p) nncfunction(p, ...
input_layer_size, ...
hidden_layer_size, ...
num_labels, X, y, lambda);

[nn_params, cost] = 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));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================= Part 9: Visualize Weights =================
% visualize what the neural network is learning by displaying the
% hidden units information

fprintf('\nVisualizing Neural Network... \n')
figure;

datademonstrate(Theta1(:, 2:end));

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ================= Part 10: Implement Predict =================
% forecast the labels.

frcst = forecast(Theta1, Theta2, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(frcst == y)) * 100);

datademonstrate.m：function [h, display_array] = datademonstrate(X, example_width)
% datademonstrate demonstrates 2D data in a nice grid
%

% 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

text
	% 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

fmincg.m：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.
%
% [introtoai] 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 + z1s; % 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+z1RHOd1) || (d2 > -SIGd1)) && (M > 0)
limit = z1; % tighten the bracket
if f2 > f1
z2 = z3 - (0.5d3z3z3)/(d3z3+f2-f3); % quadratic fit
else
A = 6(f2-f3)/z3+3*(d2+d3); % cubic fit
B = 3*(f3-f2)-z3*(d3+2d2);
z2 = (sqrt(BB-Ad2z3z3)-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, INTz3),(1-INT)z3); % don't accept too close to limits
z1 = z1 + z2; % update the step
X = X + z2s;
[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+z1RHOd1 || d2 > -SIGd1
break; % this is a failure
elseif d2 > SIGd1
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+2d2);
z2 = -d2z3z3/(B+sqrt(BB-Ad2z3z3)); % 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 > z1EXT) % extrapolation beyond limit
z2 = z1*(EXT-1.0); % set to extrapolation limit
elseif z2 < -z3INT
z2 = -z3INT;
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');

sgmd.m：function g = sgmd(z)
% sgmd calculates sigmoid functoon

g = 1.0 ./ (1.0 + exp(-z));
end

numericalgrad.m：function numgrad = numericalgrad(J, theta)
% numericalgrad calculates the gradient using "finite differences"
%a nd gives us a numerical estimate of the gradient.
%

%

numgrad = zeros(size(theta));
perturb = zeros(size(theta));
e = 1e-4;
for p = 1:numel(theta)
% Set perturbation vector
perturb(p) = e;
loss1 = J(theta - perturb);
loss2 = J(theta + perturb);
% Compute Numerical Gradient
numgrad(p) = (loss2 - loss1) / (2*e);
perturb(p) = 0;
end

end

gradcheck.m：function gradcheck(lambda)
% gradcheck creates a small neural network to check the
% backpropagation gradients

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 = initialiseweights(hidden_layer_size, input_layer_size);
Theta2 = initialiseweights(num_labels, hidden_layer_size);
% Reusing initialiseweights to generate X
X = initialiseweights(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) nncfunction(p, input_layer_size, hidden_layer_size, ...
num_labels, X, y, lambda);

[cost, grad] = costFunc(nn_params);
numgrad = numericalgrad(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 numericalgrad.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
initialiseweights.m：function W = initialiseweights(fan_out, fan_in)
% initialiseweights Initialises the weights of a layer with fan_in
% incoming connections and fan_out outgoing connections using a fixed
% strategy, this will help you later in debugging
%

% Set W to zeros
W = zeros(fan_out, 1 + fan_in);

% Initialize W using "sin", this ensures that W is always of the same
% values and will be useful for debugging
W = reshape(sin(1:numel(W)), size(W)) / 10;

% =========================================================================

end
forecast.m：function p = forecast(Theta1, Theta2, X)
% forecast forecasts the label of an input given a trained neural network
%

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

h1 = sgmd([ones(m, 1) X] * Theta1');
h2 = sgmd([ones(m, 1) h1] * Theta2');
[dummy, p] = max(h2, [], 2);

% =========================================================================

end

sgmdgrad.m：function g = sgmdgrad(z)
% sgmdgrad returns the gradient of the sigmoid function evaluated at z
%

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% hints: Calculate the gradient of the sigmoid function evaluated at
% each value of z (z can be a matrix, vector or scalar).

g = sgmd(z).*(1-sgmd(z));

% =============================================================

end
rdinitweights.m：function W = rdinitweights(L_in, L_out)
% rdinitweights Randomly initialises the weights of a layer with L_in
% incoming connections and L_out outgoing connections
%

W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Hints: Initialize W randomly so that we break the symmetry while
% training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

epsilon_init = 0.12;
W = rand(L_out, 1+ L_in)2epsilon_init - epsilon_init;

% =========================================================================

end

nncfunction.m：function [J grad] = nncfunction(nn_params, ...
input_layer_size, ...
hidden_layer_size, ...
num_labels, ...
X, y, lambda)
% nncfunction implements the neural network cost function and gradients for a
% two layer neural network which performs classification
%

% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 2 layer neural network
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));

% 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));

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the code by working through the
% following parts.
%
% Part 1: Feedforward the neural network and return the cost in the
% variable J. After implementing Part 1, you can verify that your
% cost function computation is correct by verifying the cost
% computed in lab3.m
a_1 = [ones(m, 1) X];
Z_2 = a_1Theta1';
a_2 = sgmd(Z_2);
a_2 = [ones(m, 1) a_2];
Z_3 = a_2Theta2';
a_3 = sgmd(Z_3);
h_o = a_3;

Y = zeros(m, num_labels);
for num = 1: m
Y(num, y(num)) = 1;
end;

% row = (1:m)';
% label = ones(m,1);
% Y = sparse(row, y, label);

J = sum(sum(log(h_o).(-Y)-log(1-h_o).(1-Y),2))/m;
%
% Part 2: Implement the backpropagation algorithm to compute the gradients
% Theta1_grad and Theta2_grad. You should return the partial derivatives of
% the cost function with respect to Theta1 and Theta2 in Theta1_grad and
% Theta2_grad, respectively. After implementing Part 2, you can check
% that your implementation is correct by running checkNNGradients
%
% Note: The vector y passed into the function is a vector of labels
% containing values from 1..K. You need to map this vector into a
% binary vector of 1's and 0's to be used with the neural network
% cost function.
%
% Hint: We recommend implementing backpropagation using a for-loop
% over the training examples if you are implementing it for the
% first time.
% step 2: For each output unit k in layer 3, compute the error:
delta_3 = a_3 - Y;
% step 3: For the hidden layer l=2, compute the error:
delta_2 = (Theta2(:,2:end)'*delta_3')'.*sgmdgrad(Z_2);;

% delta_2 = delta_2(:,2:end);

% step 4: Accumulate the gradient from examples
Delta_2 = delta_3'*a_2;
Delta_1 = delta_2'*a_1;
% step 5: obtain the gradient for nn cost function by average
Theta1_grad = Delta_1/m;
Theta2_grad = Delta_2/m;

% compute the regularisation terms;
Theta_grad_reg_2 = lambda/mTheta2(:, 2:end);
Theta_grad_reg_1 = lambda/mTheta1(:, 2:end);
Theta_grad_reg_2 = [zeros(size(Theta2,1),1) Theta_grad_reg_2];
Theta_grad_reg_1 = [zeros(size(Theta1,1),1) Theta_grad_reg_1];

% reged :
Theta1_grad =Theta1_grad + Theta_grad_reg_1;
Theta2_grad =Theta2_grad + Theta_grad_reg_2;
% Part 3: Implement regularization with the cost function and gradients.
%
% Hint: You can implement this around the code for
% backpropagation. That is, you can compute the gradients for
% the regularization separately and then add them to Theta1_grad
% and Theta2_grad from Part 2.
%
reg_matrix_1 = Theta1.*Theta1;
reg_nobias_1 = reg_matrix_1(:, 2:size(reg_matrix_1,2));
reg_term_1 = sum(sum(reg_nobias_1));

reg_matrix_2 = Theta2.*Theta2;
reg_nobias_2 = reg_matrix_2(:, 2:size(reg_matrix_2,2));
reg_term_2 = sum(sum(reg_nobias_2));

J = J + lambda/(2m)(reg_term_1 + reg_term_2);

% -------------------------------------------------------------

% =========================================================================

% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];

end
帮我写实验报告，别照搬实验指导书，多分析，少分点