Hessian-free NN algo summary

I'm working on an implementation for the algorithm James Martens published over the summer, and the following is a summary of what I get out of the paper. If anyone would like to critique it I would sincerely appreciate it.

(...)

EDIT:

I finally had some extra time to jump back into this. Here's a mix between Matlab syntax and pseudo-code sketch of the algorithm. Please pick it apart, or use at your own risk.

layer_dimensions = [ 784 500 2000 2000 2000 500 784 ]
ld_size          = size( layer_dimensions, 1 )
num_neurons      = sum( layer_dimensions )

% 'deal', as used here, just copies parameters; so x, rx, etc are
% duplicates of the cell
% 
% rx = R{x}
% dx = derivative of the error w/respect to x
% dphi = y * (1 - y)
[x rx dx rdx y ry dy rdy dphi] = deal( cell( ld_size,   1 ) )
[W rW dW rdW]                  = deal( cell( ld_size-1, 1 ) )

for ii = 1:ld_size
  x{ii} = zeros( layer_dimensions( ii, 1 ), num_samples )
  y{ii} = x{ii}

  if ( ii < ld_size )
    W{ii} = randn( layer_dimensions( ii:ii+1, 1 ) )
  end
end

% Initialize first layer
y{1}    = inputs
dphi{1} = zeros( size( y{1} ) )  % phi depends on x; since x{1}=0, phi{1} is
                                 % constant, so dphi{1} = 0
ry{1}   = dphi{1}                % y{1} receives no input from earlier layers

function hv = HV( V )
  % f1 pass begin
  for kk = 1:ld_size-1
    % Using material discussed near eqn (4.2) from Schraudolph's paper,
    % do a full f1 pass through M ...
    rx{kk + 1} = ( ry{kk}' * W{kk} + y{kk}' * V{kk} )'
    ry{kk + 1} = rx{kk + 1} .* dphi{kk + 1}
  end
  % f1 pass end

  % The matrix 'A' in eqn (4.2) is the identity matrix (section 2)
  %   ry{.} = A * JM * JN * V
  %
  % Now, set u = ry{ld_size} and backprop this; then we're done

  % r1 pass begin
  dy{ld_size}     = ry{ld_size}
  dx{ld_size}     = dphi{ld_size}  .* dy{ld_size}
  dW{ld_size - 1} = y{ld_size - 1} .* dx{ld_size}

  for kk = 1:ld_size-1
    dy{ld_size - kk}     = W{ld_size - kk} * dx{ld_size - kk + 1}
    dx{ld_size - kk}     = dphi{ld_size - kk} .* dy{ld_size - kk}
    dW{ld_size - 1 - kk} = y{ld_size - kk}    .* dx{ld_size - kk + 1}
  end
  % r1 pass end

  hv = dW
end

for ii = 1:numepochs
  for jj = 1:numbatches
    % f0 pass begin
    for kk = 1:ld_size-1
      x{kk + 1}    = ( y{kk}' * W{kk} )'
      y{kk + 1}    = s( x{kk + 1} )
      dphi{kk + 1} = y{kk + 1} .* (1 - y{kk + 1})
    end
    % f0 pass end

    % r1 pass begin
    dy{ld_size}     = y{ld_size} - y{1}
    dx{ld_size}     = dphi{ld_size}  .* dy{ld_size}
    dW{ld_size - 1} = y{ld_size - 1} .* dx{ld_size}

    for kk = 1:ld_size-1
      dy{ld_size - kk}     = W{ld_size - kk} * dx{ld_size - kk + 1}
      dx{ld_size - kk}     = dphi{ld_size - kk} .* dy{ld_size - kk}
      dW{ld_size - 1 - kk} = y{ld_size - kk}    .* dx{ld_size - kk + 1}
    end
    % r1 pass end

    % first candidate solution for conjugate-gradient:
    %   dW = J_LoMoN
    V = dW

    % used inside the CG loop
    normr_last = 0

    while ( not converged )
      hv = HV( V )

      % Initialize some stuff for conjugate-gradient
      if( this is the first iteration of the CG loop )
        r = V
        for kk = 1:ld_size-1
          r{kk}       = r{kk} - hv{kk}
          normr_last  = normr_last + r{kk}' * r{kk}
        end

        normr_last = sqrt( normr_last )

        p = r
      end

      % The rest from here on is a single pass of conjugate-gradient
      Hp = HV( p )
      normr = 0

      for kk = 1:ld_size-1        
        alpha   = sum(r{kk}.^2) / ( p{kk}' Hp{kk} )
        V{kk}   = V{kk} + alpha .* p{kk}
        r{kk}   = r{kk} - alpha .* Hp{kk}
        normr   = normr + r{kk}' * r{kk}
      end

      normr = sqrt( normr )

      if( normr is small enough, or max # of iterations exceeded )
        set condition to exit CG loop
      else

        beta       = normr / normr_last
        normr_last = normr

        for kk = 1:ld_size-1
          p{kk} = r{kk} + beta .* p{kk}
        end
      end % if statement

      % At this point, V contains the next best guess for the next iteration
    end % CG-loop

    % At this point, hopefully the CG algorithm converged.
    % V should contain a good weight update
  end % mini-batches
end % epochs

-Brian