Backward propagation without a loop over training examples At this point in a Coursera Machine Learning lecture on implementing neural networks, Andrew Ng says:

To implement back-prop, usually we will do that with a for-loop over 
  the training examples.
Some of you may have heard of
  frankly very advanced
  vectorization methods where you don't
  have a for-loop over the $m$ training
  examples 

The language used in the course is Octave / MATLAB.
How does one implement vectorised back-prop without a loop over the training examples?
 A: Reading the quote, it's not clear to me what exactly he means by 'loop'. But, we can consider a couple possibilities. As an example, say we want to compute the dot product between two vectors, $x$ and $y$. A naive way to do this would be to write a loop in some high level programming language to manually add up the products of the elements. In pseudo-code:
result = 0
for i = 1 to n:
    result = result + x[i] * y[i]

A more efficient approach would be to use numerical linear algebra library that implements dot products. The code might look something like:
result = dot(x, y)

In the context of high level, interpreted languages like Python, Matlab, etc. this would be called 'vectorized code'. The code we've written no longer contains a loop, although looping might still be happening at a lower level (I'll get to this in a second). In an interpreted language, the naive code would cause the interpreter to repeatedly execute the instructions in the loop. To do its job, the interpreter must perform additional computations beyond the mathematical operations we're interested in. Doing this repeatedly over the loop incurs a lot of overhead (but a smarter interpreter might use a JIT compiler to avoid this). In the vectorized code, the interpreter would pass the entire operation to the linear algebra library, which can execute it much more efficiently. In compiled languages like C, the high level code would be compiled to machine code, avoiding the overhead of an interpreter. But, using a numerical linear algebra library (e.g. BLAS) will still be more efficient than the naive code because it's highly optimized and can accelerate the computation using special features of the hardware.
The vectorized code doesn't explicitly contain loops. But, as Matthew Drury pointed out in the comments, looping might still be happening at a lower level (e.g. as the computer steps through the operations in the linear algebra library's machine code). The loop(s) won't look exactly like the naive code because much of the efficiency of numerical computing libraries comes from executing multiple operations simultaneously (i.e. parallelism). This can be achieved by taking advantage of special CPU instructions, using multiple CPU cores, or even multiple networked machines or specialized hardware.
At a fundamental level, we can consider a loop to be the repeated, sequential execution of identical operations in time. So, is it possible to perform the computation without any loops whatsoever? In principle, the answer is yes (for some problems), by using parallelism and executing all operations simultaneously. Let's ignore the specifics of existing hardware. In the dot product example, imagine we made a custom piece of hardware (e.g. an analog or digital circuit) with the following logical structure. It would compute dot products without any loops taking place:

Not all computations can be parallelized. For example, sometimes a later part of the computation depends on an earlier result. Some problems can be broken into smaller pieces that can be parallelized. Getting back to the original question about backprop, the gradient is a sum of terms computed for each training point. These terms don't depend on each other, so can be computed in parallel. If the number of terms exceeds the number of processing units, some amount of looping would be required, but fewer iterations of the loop would be needed.
In practice, high performance neural nets are often implemented using libraries like Theano and TensorFlow. These libraries allow one to describe the structure and operation of the network at a somewhat abstract level. They compile these descriptions into efficient code in a way that takes care of details like computing gradients (see automatic differentiation). They can achieve parallelism in multiple ways, including using specialized CPU instructions, multiple CPU cores, multiple networked machines, and GPUs. A current research topic and emerging trend in industry is to accelerate neural nets using custom hardware.
A: The same question came into my mind when I did the course.
From my understanding, It's mean that we can user a tensor with an extra dimension to deal with all m examples at the same time using the same vectorization approach, instead of looping over the m training examples in the backpropagation method.
Let's say we have 5000 training examples and a digit classification problems, so the output of the FNN for the looping approach would be (10,) (a vector).
Now with the full vectorized approach, the output will be (5000 x 10) matrix, each row containing the i-th training result.
A piece of python code from my Neural Network class looks like:

def forward(self, X):
    """Make the forward step, storing the value of the hypothesis
    H in each layer for backpropagation.

    I decided to store the H values (H = sigmoid(Z)) instead the Z
    as many people do for saving computation time as coping and
    assign 1 in a column (the bias term) is faster than evaluate the
    sigmoid(Z) again to compute the derivate for 2nd time.
    """
    # initial H is X + bias term
    H = np.insert(X, 0, values=1, axis=len(X.shape) > 1)  # axis 0 or 1
    self.Hs = [H]
    for Theta in self.Theta:
        # make the logistic regression for this layer
        # and store H values for backpropagation
        H = np.copy(H)
        H[:, 0] = 1

        Z = H.dot(Theta)
        Z = np.insert(Z, 0, values=1, axis=len(Z.shape) > 1)  # axis 0 or 1

        H = sigmoid(Z)     # Z --> H
        self.Hs.append(H)  # store H with all sigmoid values done

    self.H = H[:, 1:]  # remove the bias term in last step
    return self.H      # return the hypothesis

and
def grad_backprop(self, X, Y, lamb=0.0):
    "Compute the gradients using Back Propagation algorithm"
    self.forward(X)

    outputs = list()
    delta = self.H - Y
    for i in reversed(range(len(self.Theta))):
        H = self.Hs[i]
        # we need to add the full bias term
        # so make a copy and overwrite the 1st column
        H = np.copy(H)
        H[:, 0] = 1
        DTheta = np.einsum('...i,...j->...ij', H, delta)
        if len(DTheta.shape) > 2:
            # average of all gradients per sample
            DTheta = np.average(DTheta, axis=0)

        outputs.append(DTheta)

        # we can avoid the rest of the computation
        # in the last step
        if i == 0:
            break

        # sigmoid gradient
        H = self.Hs[i]
        sigrad = H * (1 - H)

        Theta = self.Theta[i]
        delta = delta.dot(Theta.T) * sigrad
        delta = delta[:, 1:]  # remove the bias

    outputs.reverse()
    outputs = np.array(outputs)

    # regularization
    if lamb > 0.0:
        m = Y.shape[0]
        Theta = np.copy(self.Theta)
        for th in Theta:
            th[0] = 0

        outputs += Theta * lamb / m

    return np.array(outputs)


I hope this helps Tom
Don't hesitate to contact me if you want to see the full code or whatever you need.
A: I'm quite stoked that I eventually worked it out:
%
% Compute gradients via backpropagation
%

% Get error of output layer
layers = 1 + length(Theta);
d{layers} = hX - Y;

% Propagate errors backwards through hidden layers
for layer = layers-1 : -1 : 2
  d{layer} = d{layer+1} * Theta{layer};
  d{layer} = d{layer}(:, 2:end); % Remove "error" for constant bias term
  d{layer} .*= sigmoidGradient(z{layer});
end

% Calculate Theta gradients
for l = 1:layers-1
  Theta_grad{l} = zeros(size(Theta{l}));

  % Sum of outer products
  Theta_grad{l} += d{l+1}' * [ones(m,1) activation{l}];

  % Add regularisation term
  Theta_grad{l}(:, 2:end) += lambda * Theta{l}(:, 2:end);
  Theta_grad{l} /= m;
end

A: i don't know if you're still looking for answer but here's a totally vectorized code for back prop in octave without the for loop and it works just as good.
a1 = X; % X already have a bias Line 44 (5000*401)
a1 = a1'; % (401*5000)
z2 = Theta1 * a1; % (25*401)*(401*5000)
a2 = sigmoid(z2); % (25*5000)

a2 = [ones(1,m) ; a2]; % adding a bias (26*5000)
z3 = Theta2 * a2; % (10*26)*(26*5000)
a3 = sigmoid(z3); % final activation layer a3 == h(theta) (10*5000)

% Step 2
delta_3 = a3 - y_new; % (10*5000)

z2=[ones(1,m); z2]; % bias (26*5000)
% Step 3
delta_2 = (Theta2' * delta_3) .* sigmoidGradient(z2); % ((26*10)*(10*5000))=(26*5000)

% Step 4
delta_2 = delta_2(2:end,:); % skipping sigma2(0) (25*5000)

Theta2_grad = Theta2_grad + delta_3 * a2'; % (10*5000)*(5000*26)
Theta1_grad = Theta1_grad + delta_2 * a1'; % (25*5000)*(5000*401)



% Step 5
Theta2_grad = (1/m) * Theta2_grad; % (10*26)
Theta1_grad = (1/m) * Theta1_grad; % (25*401)

% Regularization

% Theta1_grad(:, 1) = Theta1_grad(:, 1) ./ m; % for j = 0
% 
Theta1_grad(:, 2:end) = Theta1_grad(:, 2:end) + ((lambda/m) * Theta1(:, 2:end)); % for j >= 1 
% 
% Theta2_grad(:, 1) = Theta2_grad(:, 1) ./ m; % for j = 0
% 
Theta2_grad(:, 2:end) = Theta2_grad(:, 2:end) + ((lambda/m) * Theta2(:, 2:end)); % for j >= 1

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

