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?

  • $\begingroup$ I'm not sure what this would even mean. A for loop at the lowest level is just preforming a task over and over, like accumulating a collection of numbers into a sum. At some point in training you have to add together some numbers, which is part of the mathematics of what you are doing. Adding together a collection of numbers is fundamentally a for loop. You can write that for loop in assembly or machine code, but that just changes its form, not its essence. $\endgroup$ – Matthew Drury Jul 28 '17 at 2:12
  • $\begingroup$ On the topic of vectorization, it's important to realize that when you do something like sum(x * (y == 1)), what people call vectorized code, in a language like python or R, this is still a for loop. It's just not a for loop in python or R, it's a pre-compiled for loop in C or FORTRAN. So, without reference to a language, the concept of "is a for loop" is either under-defined, or pointless. $\endgroup$ – Matthew Drury Jul 28 '17 at 2:15
  • $\begingroup$ @MatthewDrury I've added the language - Octave. $\endgroup$ – Tom Hale Jul 28 '17 at 2:36
  • 1
    $\begingroup$ @MatthewDrury. I agree with your general point, with a small caveat: Although performing a bunch of identical operations and accumulating the results is often implemented as a loop--even at very low levels, as in your nice machine code example--I wouldn't say it must be a loop. For example, parallelizing the operations in hardware can eliminate loops (even entirely in some cases). This is a current strategy used to accelerate neural nets (e.g. GPU computation and custom hardware at the extreme end). $\endgroup$ – user20160 Jul 28 '17 at 10:57

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:

enter image description here

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.


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


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

    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)


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

        # 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 = 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.

  • $\begingroup$ FYI, I worked out how to do it. $\endgroup$ – Tom Hale Sep 15 '17 at 10:54
  • $\begingroup$ Great!! This is exactly what I mean (coded in python in my case). Simply note that you can join the 2 loops in a single loop and avoid to store the delta term d{layer} for all the layers as the Theta_grad only depends on the last delta. Congratulations! $\endgroup$ – asterio gonzalez Sep 19 '17 at 7:14
  • $\begingroup$ Good optimisation! I've moved on for now, but perhaps you'd like to post your code which does that? $\endgroup$ – Tom Hale Sep 19 '17 at 7:38

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.