function prob2
figure(1)
traj1 = prob_2a( [-1; -1], .03, .0001 );
traj2 = prob_2a( [-1; -1], .02, .0001 );
traj3 = prob_2a( [-1; -1], .005, .0001 );
line( traj1(:,1), traj1(:,2), 'color', 'r', 'linewidth', 3 );
line( traj2(:,1), traj2(:,2), 'color', 'y', 'linewidth', 3 );
line( traj3(:,1), traj3(:,2), 'color', 'w', 'linewidth', 3 );
xlabel('w1');
ylabel('w2');
title('Problem 2A');
figure(2)
momentum = prob_2b( [-1; -1], .005, .0001, 0.8 );
line( traj3(:,1), traj3(:,2), 'color', 'r', 'linewidth', 3 );
line( momentum(:,1), momentum(:,2), 'color', 'g', 'linewidth', 3 );
xlabel('w1');
ylabel('w2');
title('Problem 2B');
wMin = prob_2c( [-1; -1] );
disp( 'Newton''s Method' )
disp(sprintf('w1:  %e', wMin(1)))
disp(sprintf('w2:  %e', wMin(2)))
return





function trajectory = prob_2a(           initialWeight,                    ...
                                         learningRate,                     ...
                                         stopMagnitude)

%------------------------------------------------------------------------------
% This script will perform a steepest-descent iteration of the cost 
% function given below and return the weight vectors that were encountered 
% during iteration, ending with the weight vector that gave the cost function
% minimum.
%
% If the simulation does not converge, the weight vectors that were computed 
% will be returned.
%
%------------------------------------------------------------------------------
%   Inputs
%------------------------------------------------------------------------------
% initialWeight     the initial weight vector.  The script will begin iterating 
%                   at the cost function corresponding to this weight.
%
% learningRate      the amount of the gradient that will be applied to each 
%                   iteration.  A large value may cause oscillations around 
%                   the minimum, while a small value may result in slow 
%                   computation time.
%
% stopMagnitude     iteration will halt when the magnitude of the correction 
%                   becomes smaller than this value.
%------------------------------------------------------------------------------
%   Outputs
%------------------------------------------------------------------------------
% trajectory        a two-column matrix in which successive rows are successive 
%                   values of the weight vector as the iteration progresses.




% Just in case I screwed something up, I'll make a maximum iteration count
maxIterations = 800;



% I'm going to do most of this symbolically since I can.  Define an arbitrary
% weight vector
syms w1 w2;
w = [   w1;                                                                 ...
        w2];



% Define the cost function in terms of the weight vector
costFunction = w(1)^2  +  30*w(2)^2  +  5*w(1)*w(2);



% Find the gradient of the cost function
grad = [    diff(costFunction, w(1));                                       ...
            diff(costFunction, w(2)) ];

    
    
% Do the iteration
wMin = initialWeight;
trajectory = [ wMin(1) wMin(2) ];
iteration = 0;
while ( iteration < maxIterations );
    
    % Find the update to the weight vector
    w1 = wMin(1);
    w2 = wMin(2);
    correction = -learningRate * eval( grad );
    wMin = wMin + correction;
    trajectory = [ trajectory; wMin(1) wMin(2) ];
    
    % If the magnitude of the correction is smaller than stopMagnitude, 
    % we're done
    if sqrt(dot(correction, correction)) <= stopMagnitude;
        break;
    end;
    
    iteration = iteration + 1;
end;
if iteration >= maxIterations;
    disp 'The iteration did not converge.  Try lowering the learning rate.';
end;



% Plot the cost function - the bounds will be based upon the minimum and 
% maximum values that were found while computing the trajectory.  The bounds 
% will be 'spread' so that the trajectory does not lie right on the edge of 
% the plot.
spread = 1.6;
w1_min = min(trajectory(:,1));
w1_max = max(trajectory(:,1));
w1_spread = spread * (w1_max - w1_min);
w1_min = w1_min - ( w1_spread - (w1_max - w1_min) ) /2;
w2_min = min(trajectory(:,2));
w2_max = max(trajectory(:,2));
w2_spread = spread * (w2_max - w2_min);
w2_min = w2_min - ( w2_spread - (w2_max - w2_min) ) /2;
discretizationCount = 100;
weight1 = w1_min: w1_spread/discretizationCount: w1_min + w1_spread;
weight2 = w2_min: w2_spread/discretizationCount: w2_min + w2_spread;
for i=1:discretizationCount+1;       for j=1:discretizationCount+1;
    w1 = weight1(i);
    w2 = weight2(j);
    evaluatedCost(i,j) = eval(costFunction);
end; end;
contourf(weight1, weight2, evaluatedCost, 30);
line( trajectory(:,1), trajectory(:,2), 'color', 'r' );
disp(sprintf('Iterations for a learning rate of %1.4f: %d',                ...
                learningRate, iteration + 1))
return






function trajectory = prob_2b(           initialWeight,                    ...
                                         learningRate,                     ...
                                         stopMagnitude,                    ...
                                         inertia)

%------------------------------------------------------------------------------
% This script will perform a steepest-descent iteration of the cost 
% function given below and return the weight vectors that were encountered 
% during iteration, ending with the weight vector that gave the cost function
% minimum.
%
% If the simulation does not converge, the weight vectors that were computed 
% will be returned.
%
%------------------------------------------------------------------------------
%   Inputs
%------------------------------------------------------------------------------
% initialWeight     the initial weight vector.  The script will begin iterating 
%                   at the cost function corresponding to this weight.
%
% learningRate      the amount of the gradient that will be applied to each 
%                   iteration.  A large value may cause oscillations around 
%                   the minimum, while a small value may result in slow 
%                   computation time.
%
% stopMagnitude     iteration will halt when the magnitude of the correction 
%                   becomes smaller than this value.
%
% inertia           an inertial weighting factor used to speed up calculation
%------------------------------------------------------------------------------
%   Outputs
%------------------------------------------------------------------------------
% trajectory        a two-column matrix in which successive rows are successive 
%                   values of the weight vector as the iteration progresses.




% Just in case I screwed something up, I'll make a maximum iteration count
maxIterations = 800;



% I'm going to do most of this symbolically since I can.  Define an arbitrary
% weight vector
syms w1 w2;
w = [   w1;                                                                 ...
        w2];



% Define the cost function in terms of the weight vector
costFunction = w(1)^2  +  30*w(2)^2  +  5*w(1)*w(2);



% Find the gradient of the cost function
grad = [    diff(costFunction, w(1));                                       ...
            diff(costFunction, w(2)) ];

    
    
% Do the iteration
wMin = initialWeight;
trajectory = [ wMin(1) wMin(2) ];
iteration = 0;
lastCorrection = 0;
while ( iteration < maxIterations );
    
    % Find the update to the weight vector
    w1 = wMin(1);
    w2 = wMin(2);
    correction = -learningRate * eval( grad )  +  inertia * lastCorrection;
    lastCorrection = correction;
    wMin = wMin + correction;
    trajectory = [ trajectory; wMin(1) wMin(2) ];
    
    % If the magnitude of the correction is smaller than stopMagnitude, 
    % we're done
    if sqrt(dot(correction, correction)) <= stopMagnitude;
        break;
    end;
    
    iteration = iteration + 1;
end;
if iteration >= maxIterations;
    disp 'The iteration did not converge.  Try lowering the learning rate.';
end;



% Plot the cost function - the bounds will be based upon the minimum and 
% maximum values that were found while computing the trajectory.  The bounds 
% will be 'spread' so that the trajectory does not lie right on the edge of 
% the plot.
spread = 1.6;
w1_min = min(trajectory(:,1));
w1_max = max(trajectory(:,1));
w1_spread = spread * (w1_max - w1_min);
w1_min = w1_min - ( w1_spread - (w1_max - w1_min) ) /2;
w2_min = min(trajectory(:,2));
w2_max = max(trajectory(:,2));
w2_spread = spread * (w2_max - w2_min);
w2_min = w2_min - ( w2_spread - (w2_max - w2_min) ) /2;
discretizationCount = 100;
weight1 = w1_min: w1_spread/discretizationCount: w1_min + w1_spread;
weight2 = w2_min: w2_spread/discretizationCount: w2_min + w2_spread;
for i=1:discretizationCount+1;       for j=1:discretizationCount+1;
    w1 = weight1(i);
    w2 = weight2(j);
    evaluatedCost(i,j) = eval(costFunction);
end; end;
contourf(weight1, weight2, evaluatedCost, 30);
%line( trajectory(:,1), trajectory(:,2), 'color', 'r' );
disp(sprintf('Iterations for a learning rate of %1.4f and momentum of %1.4f: %d',                ...
                learningRate, inertia, iteration + 1))
return





function wMin = prob_2c(     initialWeight      )

%------------------------------------------------------------------------------
% This script will perform Newton's method to find the weight vector that minimizes
% the cost function given below.  The script returns the weight vector that 
% minimized the cost function.
%
% If the simulation does not converge, the weight vectors that were computed 
% will be returned.
%
%------------------------------------------------------------------------------
%   Inputs
%------------------------------------------------------------------------------
% initialWeight     the initial weight vector.  The script will begin iterating 
%                   at the cost function corresponding to this weight.
%
%------------------------------------------------------------------------------
%   Outputs
%------------------------------------------------------------------------------
% wMin              the weight vector that minimizes the cost function





% I'm going to do most of this symbolically since I can.  Define an arbitrary
% weight vector
syms w1 w2;
w = [   w1;                                                                 ...
        w2];



% Define the cost function in terms of the weight vector
costFunction = w(1)^2  +  30*w(2)^2  +  5*w(1)*w(2);



% Find the Hessian and the gradient of the cost function
d1 = diff( costFunction, w(1) );
d2 = diff( costFunction, w(2) );
hessian = [     diff( d1,w(1) )             diff( d2,w(1) );                ...
                diff( d1,w(2) )             diff( d2,w(2) )];
invHessian = inv(hessian);
grad = [        d1;                                                         ...
                d2 ];


    
   
% Find the weight change - since our cost function is quadratic, no iteration is 
% necessary
w1 = initialWeight(1);
w2 = initialWeight(2);
correction = - eval(invHessian)  *  eval(grad);
wMin = initialWeight + correction;
return