Python and Vectorization

deep-learning
vectorization
numpy
python
Replacing for loops with NumPy operations, fully vectorized logistic regression, broadcasting, rank 1 array pitfalls, and the maximum likelihood view of the cost.
Published

Jul 30, 2026

The previous section ended with a gradient descent implementation that needed two explicit for loops, and promised a way to get rid of them. This section delivers on that promise. You will see what vectorization is, why it can make your code hundreds of times faster, how to vectorize logistic regression end to end, how broadcasting works, and how to avoid a family of subtle NumPy bugs. It closes with an optional justification of where the logistic regression cost function comes from. You met the basic idea of vectorization in the machine learning course; here we push it much further.

Vectorization

Vectorization is basically the art of getting rid of explicit for loops in your code. In the deep learning era, you often find yourself training on relatively large datasets, because that is when deep learning algorithms tend to shine. So it is important that your code runs very quickly, because otherwise, when you are training on a big dataset, your code might take a long time to run, and you just find yourself waiting a very long time to get the result. The ability to perform vectorization has become a key skill.

Let us start with an example. In logistic regression you need to compute \(z = w^T x + b\), where \(w\) and \(x\) are both \(n_x\) dimensional column vectors, maybe very large vectors if you have a lot of features. A non-vectorized implementation would compute the sum one term at a time.

z = 0
for i in range(n_x):
    z += w[i] * x[i]
z += b

You would find that this is really slow. In contrast, a vectorized implementation just computes \(w^T x\) directly.

z = np.dot(w, x) + b

This is much faster. Let us actually illustrate this with a little demo. Create two arrays of a million random values each, then compute their dot product twice, once with np.dot and once with an explicit for loop, timing each version.

import numpy as np
import time

a = np.random.rand(1000000)
b = np.random.rand(1000000)

tic = time.time()
c = np.dot(a, b)
toc = time.time()
print("Vectorized version: " + str(1000 * (toc - tic)) + " ms")

c = 0
tic = time.time()
for i in range(1000000):
    c += a[i] * b[i]
toc = time.time()
print("Explicit for loop: " + str(1000 * (toc - tic)) + " ms")

Here is that exact demo executed while rendering this page, printing the value of c in both cases to make sure they compute the same thing.

Vectorized version: c = 249825.02, computed in 0.10 ms
Explicit for loop:  c = 249825.02, computed in 116.23 ms
The for loop took about 1183 times longer.

Both versions compute the same value, but the explicit for loop takes on the order of hundreds of times longer than the vectorized version (in the original demo, about 1.5 milliseconds against nearly 500 milliseconds, roughly 300 times slower). The exact numbers vary a little from run to run and machine to machine, but the gap stays enormous. A factor of 300 is the difference between your code taking one minute to run and taking five hours to run. When you are implementing deep learning algorithms, you get your result back much faster if you vectorize your code.

You might have heard that a lot of scalable deep learning implementations are done on a GPU (graphics processing unit), but the demo above ran on a CPU. It turns out that both GPUs and CPUs have parallelization instructions, sometimes called SIMD instructions, which stands for single instruction, multiple data. What this basically means is that if you use built-in functions such as np.dot, which do not require you to explicitly implement a for loop, NumPy can take much better advantage of parallelism to do the computations much faster. This is true both on CPUs and on GPUs. GPUs are remarkably good at these SIMD calculations, but CPUs are actually not too bad at them either, maybe just not as good as GPUs.

The rule of thumb to remember is, whenever possible, avoid using explicit for loops.

Review Questions

1. What is vectorization, and why has it become so important in deep learning?

Vectorization is the art of getting rid of explicit for loops in your code, replacing them with operations like np.dot that work on whole vectors at once. Deep learning shines on large datasets, and on large datasets non-vectorized code can take hundreds of times longer, the difference between a minute and hours of waiting.

1. Roughly how much slower was the explicit for loop than np.dot in the million-element dot product demo?

About 300 times slower (around 1.5 milliseconds for the vectorized version against nearly 500 milliseconds for the loop). Both computed the same value.

1. What are SIMD instructions, and what do they have to do with NumPy built-in functions?

SIMD stands for single instruction, multiple data, parallelization instructions available on both CPUs and GPUs. When you use built-in functions like np.dot instead of explicit for loops, NumPy can take much better advantage of this parallelism, which is what makes the vectorized version so much faster.

More Vectorization Examples

The rule of thumb to keep in mind is, when you are programming your neural networks, or even just a regression, whenever possible avoid explicit for loops. It is not always possible to never use a for loop, but when you can use a built-in function or find some other way to compute whatever you need, you will often go faster than with an explicit for loop.

Matrix Times Vector

Say you want to compute a vector \(u\) as the product of a matrix \(A\) and a vector \(v\). By the definition of matrix multiplication, \(u_i = \sum_j A_{ij} v_j\). A non-vectorized implementation needs two for loops, over both \(i\) and \(j\).

u = np.zeros((n, 1))
for i in range(n):
    for j in range(n):
        u[i] += A[i][j] * v[j]

The vectorized implementation eliminates both loops and is going to be way faster.

u = np.dot(A, v)

Element-Wise Functions

Say you already have a vector \(v\) in memory and you want to apply the exponential operation to every element, producing \(u = (e^{v_1}, e^{v_2}, \ldots, e^{v_n})\). The non-vectorized implementation initializes \(u\) to zeros and computes the elements one at a time.

u = np.zeros((n, 1))
for i in range(n):
    u[i] = math.exp(v[i])

But NumPy has many built-in functions that compute these vectors with a single function call.

u = np.exp(v)

With just one line of code, \(v\) as the input vector and \(u\) as the output vector, you have gotten rid of the explicit for loop. NumPy has many of these vector-valued functions.

NumPy element-wise operations
Call What it computes
np.exp(v) element-wise exponential
np.log(v) element-wise logarithm
np.abs(v) element-wise absolute value
np.maximum(v, 0) element-wise maximum of each element with 0
v**2 element-wise square
1/v element-wise inverse

So whenever you are tempted to write a for loop, take a look and see if there is a way to call a NumPy built-in function to do it without the for loop.

Getting Rid of One Loop in Logistic Regression

Let us take these learnings and apply them to the gradient descent implementation from the previous section, and see if we can at least get rid of one of the two for loops. The inner loop was over the features. With just \(n_x = 2\) features we wrote dw1 and dw2 by hand, but with more features you would need a for loop over dw1, dw2, dw3, and so on.

To eliminate this second loop, instead of explicitly initializing dw1, dw2, and so on to zeros, make dw a vector, and update it with one vector-valued operation per example.

dw = np.zeros((n_x, 1))

for i = 1 to m:
    # forward pass and dz(i) as before ...
    dw += x(i) dz(i)
    db += dz(i)

dw /= m

Now we have gone from two for loops to just one. We still have the loop over the individual training examples, but by getting rid of one for loop the code already runs faster. It turns out we can do even better, and process the entire training set without a single for loop.

Review Questions

1. How many for loops does the naive matrix-vector product \(u = Av\) need, and what replaces them?

Two, one over \(i\) and one over \(j\), since \(u_i = \sum_j A_{ij} v_j\). The single call u = np.dot(A, v) eliminates both.

1. Name some NumPy vector-valued functions that replace element-wise for loops.

np.exp(v), np.log(v), np.abs(v), np.maximum(v, 0) for the element-wise max with 0, v**2 for the element-wise square, and 1/v for the element-wise inverse.

1. How does making dw a vector remove one of the two for loops in the logistic regression gradient computation?

Instead of initializing and updating each dwj separately in a loop over the \(n_x\) features, you initialize dw = np.zeros((n_x, 1)) and replace the inner loop with the single vector operation dw += x(i) dz(i), followed by dw /= m after the loop. Only the loop over the \(m\) training examples remains.

1. Consider the array a = np.array([[2, 1], [1, 3]]). What is the result of a * a?

  1. \(\begin{bmatrix} 5 & 5 \\ 5 & 10 \end{bmatrix}\)
  2. \(\begin{bmatrix} 4 & 2 \\ 2 & 6 \end{bmatrix}\)
  3. \(\begin{bmatrix} 4 & 1 \\ 1 & 9 \end{bmatrix}\)
  4. The computation cannot happen because the sizes do not match. It is going to be an error.

c. The * operator is element-wise multiplication, so every entry is squared in place. Option a. is the matrix product \(a \cdot a\), which is what np.dot(a, a) would compute instead.

Vectorizing Logistic Regression

Now let us vectorize the implementation of logistic regression, so it can process an entire training set, that is, implement a single iteration of gradient descent with respect to the whole training set, without using even a single explicit for loop.

First examine the forward propagation steps of logistic regression. If you have \(m\) training examples, then to make a prediction on the first example you compute

\[ z^{(1)} = w^T x^{(1)} + b \qquad a^{(1)} = \sigma\big(z^{(1)}\big) \]

Then to make a prediction on the second example you compute \(z^{(2)}\) and \(a^{(2)}\), on the third \(z^{(3)}\) and \(a^{(3)}\), and so on. You might need to do this \(m\) times.

It turns out there is a way to compute all the predictions without an explicit for loop. Remember that we defined the matrix \(X\) to be the training inputs stacked together in columns, an \(n_x \times m\) matrix (as a Python NumPy shape, (n_x, m)). Now construct a \(1 \times m\) row vector holding \(z^{(1)}, z^{(2)}, \ldots, z^{(m)}\) all at the same time,

\[ Z = \begin{bmatrix} z^{(1)} & z^{(2)} & \cdots & z^{(m)} \end{bmatrix} = w^T X + \begin{bmatrix} b & b & \cdots & b \end{bmatrix} \]

Here is why this works. \(w^T\) is a row vector, so by the rules of matrix multiplication,

\[ w^T X = \begin{bmatrix} w^T x^{(1)} & w^T x^{(2)} & \cdots & w^T x^{(m)} \end{bmatrix} \]

and adding the row vector of \(b\)s adds \(b\) to each element. The first element is then exactly the definition of \(z^{(1)}\), the second element is exactly \(z^{(2)}\), and so on. Just as \(X\) was obtained by taking the training examples \(x^{(i)}\) and stacking them horizontally, \(Z\) is obtained by stacking the lowercase \(z^{(i)}\) values horizontally. In NumPy, the command is

Z = np.dot(w.T, X) + b

There is a subtlety in Python here. \(b\) is a real number (or if you want, a \(1 \times 1\) matrix), but when you add it to the row vector, Python automatically takes the real number \(b\) and expands it out to a \(1 \times m\) row vector. In case this operation seems a little bit mysterious, it is called broadcasting in Python, and you do not have to worry about it for now. The Broadcasting in Python section below covers it in more detail. The takeaway is that with this one line of code, you calculate \(Z\), a \(1 \times m\) matrix containing all of the lowercase \(z\) values.

How about the activations? Stacking the lowercase \(a^{(i)}\) values horizontally gives a new variable,

\[ A = \begin{bmatrix} a^{(1)} & a^{(2)} & \cdots & a^{(m)} \end{bmatrix} = \sigma(Z) \]

With a vector-valued sigmoid function, one that takes the matrix \(Z\) as input and efficiently outputs the matrix \(A\), you compute all the activations at the same time (you will implement such a sigmoid in a programming exercise).

So to recap, instead of needing to loop over the \(m\) training examples to compute \(z^{(i)}\) and \(a^{(i)}\) one at a time, you compute all the \(z\) values with one line of code and all the \(a\) values with another. That is a vectorized implementation of the forward propagation for all \(m\) training examples at the same time. It turns out vectorization can also compute the backward propagation, the gradients, just as efficiently. That comes next.

Review Questions

1. Write the vectorized formula that computes all \(m\) values of \(z\) at once, and its NumPy implementation.

\[ Z = w^T X + \begin{bmatrix} b & b & \cdots & b \end{bmatrix} \] where \(X\) is the \(n_x \times m\) matrix of inputs stacked in columns and \(Z\) is \(1 \times m\). In NumPy, Z = np.dot(w.T, X) + b.

1. In Z = np.dot(w.T, X) + b, the variable b is a single real number. Why does the addition still work?

Python automatically expands the real number \(b\) into a \(1 \times m\) row vector so it can be added element by element. This behavior is called broadcasting.

1. After computing \(Z\), how are all the activations computed without a loop?

\(A = \sigma(Z)\), applying a vector-valued sigmoid to the whole \(1 \times m\) matrix \(Z\) at once, giving the matrix \(A\) that stacks \(a^{(1)}\) through \(a^{(m)}\) horizontally.

Vectorizing the Gradient Output

You have seen how to use vectorization to compute the predictions, the lowercase \(a\) values, for an entire training set all at the same time. Now let us use vectorization to also perform the gradient computations for all \(m\) training examples at the same time, and then put it all together into a very efficient implementation of logistic regression.

Computing All the dz Values at Once

For the gradient computation, we computed \(dz^{(1)} = a^{(1)} - y^{(1)}\) for the first example, \(dz^{(2)} = a^{(2)} - y^{(2)}\) for the second, and so on for all \(m\) training examples. So define a new variable that stacks all the lowercase \(dz\) variables horizontally,

\[ dZ = \begin{bmatrix} dz^{(1)} & dz^{(2)} & \cdots & dz^{(m)} \end{bmatrix} \]

a \(1 \times m\) matrix, or alternatively an \(m\) dimensional row vector. We already know how to compute the matrix \(A\), which stacks \(a^{(1)}\) through \(a^{(m)}\), and we defined \(Y\) as \(y^{(1)}\) through \(y^{(m)}\), also stacked horizontally. Based on these definitions, maybe you can see for yourself that

\[ dZ = A - Y = \begin{bmatrix} a^{(1)} - y^{(1)} & a^{(2)} - y^{(2)} & \cdots & a^{(m)} - y^{(m)} \end{bmatrix} \]

The first element is exactly the definition of \(dz^{(1)}\), the second element is exactly \(dz^{(2)}\), and so on. So with just one line of code, dZ = A - Y, you compute all of the \(dz\) values at the same time.

Removing the Last Training Loop

In the previous implementation we had already gotten rid of one for loop, making dw a vector, but we still had the loop over the \(m\) training examples that accumulated dw += x(i) dz(i) and db += dz(i) before dividing both by \(m\). Let us vectorize those operations too.

For db, the computation is basically summing up all the \(dz\) values and dividing by \(m\),

\[ db = \frac{1}{m} \sum_{i=1}^{m} dz^{(i)} \]

and since all the \(dz\) values sit in the row vector \(dZ\), in Python this is

db = np.sum(dZ) / m

How about dw? It turns out to be

\[ dw = \frac{1}{m} X \, dZ^T \]

Here is why that is the right thing. \(X\) is the matrix with \(x^{(1)}\) through \(x^{(m)}\) stacked up in columns, and \(dZ^T\) is a column of the values \(dz^{(1)}\) down to \(dz^{(m)}\). If you work out what this matrix times this vector is, it turns out to be

\[ \frac{1}{m} X \, dZ^T = \frac{1}{m} \Big( x^{(1)} dz^{(1)} + \cdots + x^{(m)} dz^{(m)} \Big) \]

an \(n_x \times 1\) vector, which is exactly what dw was computing, taking the \(x^{(i)} dz^{(i)}\) terms and adding them up. So the matrix-vector multiplication does the accumulation for you, and with one line of code you compute dw.

Highly Efficient Logistic Regression

Now let us put it all together. A single iteration of gradient descent for logistic regression, processing the entire training set, becomes

Z = np.dot(w.T, X) + b
A = sigma(Z)

dZ = A - Y
dw = np.dot(X, dZ.T) / m
db = np.sum(dZ) / m

w = w - alpha * dw
b = b - alpha * db

The first two lines are the forward propagation, computing the predictions on all \(m\) examples, and the next three are the backpropagation, computing all the derivatives, without using a single for loop. Then the gradient descent update is \(w := w - \alpha \, dw\) and \(b := b - \alpha \, db\).

One caveat. We said to get rid of explicit for loops whenever you can, but if you want to implement multiple iterations of gradient descent, you still need a for loop over the number of iterations. If you want a thousand iterations of gradient descent, there is an outermost for loop over the iteration number, and there is no way to get rid of that one. But it is incredibly cool that you can implement at least one full iteration of gradient descent without needing a for loop.

There is just one more detail. The implementation above briefly leaned on the technique called broadcasting, which Python and NumPy let you use to make certain parts of your code much more efficient. Let us look at it properly.

Review Questions

1. How are all \(m\) values of \(dz\) computed with one line of code?

\(dZ = A - Y\), where \(A\) and \(Y\) stack the \(a^{(i)}\) and \(y^{(i)}\) values horizontally into \(1 \times m\) row vectors. Element \(i\) of the difference is \(a^{(i)} - y^{(i)}\), which is exactly the definition of \(dz^{(i)}\).

1. Give the vectorized formulas for dw and db, and explain why the dw formula is correct.

\[ dw = \frac{1}{m} X \, dZ^T \qquad db = \frac{1}{m} \sum_{i=1}^{m} dz^{(i)} \;=\; \texttt{np.sum(dZ) / m} \] \(X\) has the \(x^{(i)}\) stacked in columns and \(dZ^T\) is the column of \(dz^{(i)}\) values, so the product works out to \(\frac{1}{m}\big(x^{(1)} dz^{(1)} + \cdots + x^{(m)} dz^{(m)}\big)\), an \(n_x \times 1\) vector, exactly the accumulation the old for loop performed.

1. How do you compute the derivative of \(b\) in one line of code in Python NumPy?

  1. 1 - m(np.sum(dZ))
  2. 1 * m(np.sum(dZ))
  3. m(np.sum(dZ))
  4. 1 / m * (np.sum(dZ))

d. The derivative is the average of the \(dz^{(i)}\) values, \(db = \frac{1}{m} \sum_{i=1}^{m} dz^{(i)}\), and np.sum(dZ) computes the sum over the row vector \(dZ\), so multiplying by \(\frac{1}{m}\) gives the average. The other options subtract or multiply by \(m\) instead of dividing by it, so they do not compute an average at all.

1. After fully vectorizing one iteration of gradient descent, which for loop can you still not get rid of?

The outermost loop over the number of gradient descent iterations. If you want a thousand iterations, you still need a for loop over the iteration number. Everything within one iteration, forward propagation and backpropagation over all \(m\) examples, runs without any for loop.

Broadcasting in Python

Broadcasting is another technique that you can use to make your Python code run faster. Let us explore how broadcasting in Python actually works, with an example.

Calories Example

This matrix shows the number of calories from carbohydrates, proteins, and fats in 100 grams of four different foods.

\[ \begin{array}{ccl} & \begin{array}{cccc} \;\;\text{Apples} & \text{Beef} & \text{Eggs} & \text{Potatoes} \end{array} & \\ A = & \begin{bmatrix} \;\;56.0 & 0.0 & 4.4\phantom{0} & 68.0\;\; \\ \;\;1.2 & 104.0 & 52.0 & 8.0\;\; \\ \;\;1.8 & 135.0 & 99.0 & 0.9\;\; \end{bmatrix} & \begin{array}{l} \text{Carb} \\ \text{Protein} \\ \text{Fat} \end{array} \end{array} \]

For example, 100 grams of apples has 56 calories from carbs, and much less from proteins and fats, whereas 100 grams of beef has 104 calories from protein and 135 calories from fat. Say your goal is to calculate the percentage of calories from carbs, proteins, and fats for each of the four foods. If you add up the numbers in the apple column, you get \(56 + 1.2 + 1.8 = 59\) calories, so as a percentage, calories from carbohydrates in an apple are \(56 / 59 \approx 94.9\%\). Most of the calories in an apple come from carbs, whereas most of the calories in beef come from protein and fat.

The calculation you want is to sum each of the four columns to get the total calories for each food, and then divide throughout the matrix to get the percentages. The question is, can you do this without an explicit for loop? It takes two lines of Python.

cal = A.sum(axis=0)
percentage = 100 * A / cal.reshape(1, 4)

Here is the computation executed while rendering this page.

A =
[[ 56.    0.    4.4  68. ]
 [  1.2 104.   52.    8. ]
 [  1.8 135.   99.    0.9]]

cal = [ 59.  239.  155.4  76.9]

percentage =
[[94.9  0.   2.8 88.4]
 [ 2.  43.5 33.5 10.4]
 [ 3.1 56.5 63.7  1.2]]

The totals come out to 59 calories for the apple, 239 for the beef, and so on, and in the percentage matrix the first column confirms that 94.9% of the apple calories are from carbs.

A couple of details about those two lines.

  • The parameter axis=0 means to sum vertically, down the columns. The horizontal axis is axis 1, so summing across a row would use axis=1.
  • The division is an example of Python broadcasting. \(A\) is a \(3 \times 4\) matrix and it is divided by the \(1 \times 4\) matrix cal.reshape(1, 4). Technically, after the first line cal is already a \(1 \times 4\) matrix, so the reshape call is a little bit redundant. But when you are not entirely sure of the dimensions of a matrix, calling reshape to make sure it is the row or column vector you want is a constant-time, very cheap operation, so do not be shy about using it.

How Broadcasting Works

So how can you divide a \(3 \times 4\) matrix by a \(1 \times 4\) matrix? Let us go through a few more examples.

If you take a \(4 \times 1\) vector and add it to a number, Python auto-expands the number into a \(4 \times 1\) vector, adding 100 to every element,

\[ \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} + 100 = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} + \begin{bmatrix} 100 \\ 100 \\ 100 \\ 100 \end{bmatrix} = \begin{bmatrix} 101 \\ 102 \\ 103 \\ 104 \end{bmatrix} \]

This type of broadcasting works with both column vectors and row vectors. In fact, we used this form of broadcasting earlier, when the constant added to a vector was the parameter \(b\) in logistic regression.

Here is another example. If you have a \(2 \times 3\) matrix and add a \(1 \times 3\) matrix, Python copies the smaller matrix down twice to make it \(2 \times 3\), then adds element-wise, so 100 is added to the first column, 200 to the second, and 300 to the third,

\[ \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix} + \begin{bmatrix} 100 & 200 & 300 \end{bmatrix} = \begin{bmatrix} 101 & 202 & 303 \\ 104 & 205 & 306 \end{bmatrix} \]

This is basically what happened in the calories example, except with division instead of addition. One last example. If you have an \(m \times n\) matrix and add an \(m \times 1\) matrix, the column gets copied horizontally \(n\) times, here adding 100 to the first row and 200 to the second,

\[ \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix} + \begin{bmatrix} 100 \\ 200 \end{bmatrix} = \begin{bmatrix} 101 & 102 & 103 \\ 204 & 205 & 206 \end{bmatrix} \]

Here is the general principle of broadcasting in Python. If you have an \((m, n)\) matrix and you add, subtract, multiply, or divide with a \((1, n)\) matrix, the \((1, n)\) matrix is copied \(m\) times into an \((m, n)\) matrix, and the operation is applied element-wise. If instead you combine the \((m, n)\) matrix with an \((m, 1)\) matrix, that matrix is copied \(n\) times into an \((m, n)\) matrix, and again the operation applies element-wise. And if you take an \((m, 1)\) column vector and combine it with a real number (a \((1, 1)\) matrix), the number is copied \(m\) times so the operation applies to every element, with something similar for row vectors.

The fully general version of broadcasting can do even a little bit more than this. If you are interested, you can read the NumPy documentation on broadcasting for the slightly more general definition, but the forms above are the main ones you need when implementing a neural network.

For those used to programming in MATLAB or Octave, the function bsxfun does something similar to broadcasting, not quite the same, but it is often used for a similar purpose. This is really only for very advanced MATLAB and Octave users. If you have not heard of it, do not worry about it, you do not need it when coding neural networks in Python.

Broadcasting not only makes your code run faster, it also helps you get what you want done with fewer lines of code.

Review Questions

1. Which of the following NumPy lines of code would sum the values in a matrix A vertically?

  1. A.sum(axis)
  2. A.sum(axis=0)
  3. A.sum(axis=1)

b. axis=0 sums vertically, down the columns (giving the total calories of each food in the example). axis=1 sums horizontally, across each row. Option a. is not valid on its own, since axis is not defined.

1. In the calories example, the reshape call in A / cal.reshape(1, 4) was technically redundant. Why use it anyway?

cal was already a \(1 \times 4\) matrix, but when you are not entirely sure of a matrix’s dimensions, calling reshape guarantees the shape you want. It is a constant-time, very cheap operation, so there is no reason to be shy about using it.

1. State the general principle of broadcasting for an \((m, n)\) matrix.

When an \((m, n)\) matrix is added, subtracted, multiplied, or divided with a \((1, n)\) matrix, the \((1, n)\) matrix is copied \(m\) times to become \((m, n)\) and the operation applies element-wise. With an \((m, 1)\) matrix, it is copied \(n\) times horizontally and applied element-wise. A vector plus a real number copies the number into every position the same way.

1. Where had we already used broadcasting before this section named it?

In the vectorized forward pass Z = np.dot(w.T, X) + b, where the real number \(b\) is automatically expanded into a \(1 \times m\) row vector before the addition.

1. Consider the random arrays a = np.random.randn(3, 4) and b = np.random.randn(1, 4), and c = a + b. What will be the shape of c?

  1. c.shape = (3, 1)
  2. c.shape = (1, 4)
  3. c.shape = (3, 4)
  4. The computation cannot happen because it is not possible to broadcast more than one dimension.

c. Broadcasting is used, so the \((1, 4)\) row b is copied 3 times and summed with each row of a, giving a \((3, 4)\) result.

1. Consider the random arrays a = np.random.randn(4, 3) and b = np.random.randn(1, 3), and c = a * b. What will be the shape of c?

  1. c.shape = (4, 3)
  2. c.shape = (1, 3)
  3. The computation cannot happen because it is not possible to broadcast more than one dimension.
  4. The computation cannot happen because the sizes do not match.

a. Broadcasting is invoked, so the \((1, 3)\) row b is multiplied element-wise with each of the 4 rows of a, giving a \((4, 3)\) result. The general principle applies to addition, subtraction, multiplication, and division alike.

1. Suppose a.shape = (4, 3) and b.shape = (4, 1), and consider this code.

for i in range(3):
    for j in range(4):
        c[i][j] = a[j][i] + b[j]

How do you vectorize this?

  1. c = a.T + b.T
  2. c = a + b.T
  3. c = a.T + b
  4. c = a + b

a. The result c has shape \((3, 4)\). Reading a[j][i] where c[i][j] is written means rows and columns are swapped, so we need a.T, which is \((3, 4)\). And b[j] varies along the column index j of c, so the \((4, 1)\) column b must become the \((1, 4)\) row b.T, which then broadcasts down the 3 rows.

1. Consider a = np.random.randn(3, 3), b = np.random.randn(3, 1), and c = a * b. What will c be?

  1. It invokes broadcasting, so b is copied three times to become \((3, 3)\), and * invokes a matrix multiplication operation of two \(3 \times 3\) matrices, so c.shape will be \((3, 3)\).
  2. It invokes broadcasting, so b is copied three times to become \((3, 3)\), and * is an element-wise product, so c.shape will be \((3, 3)\).
  3. It leads to an error, since * cannot operate on these two matrices. You need to use np.dot(a, b) instead.
  4. It multiplies the \(3 \times 3\) matrix a with the \(3 \times 1\) vector b, resulting in a \(3 \times 1\) vector, so c.shape = (3, 1).

b. The \((3, 1)\) column b is broadcast (copied three times horizontally) into a \((3, 3)\) matrix, and * performs an element-wise product, never a matrix multiplication. Matrix multiplication would be np.dot(a, b), which is what option d. describes.

Note on Python and NumPy Vectors

The great flexibility of Python and NumPy, including broadcasting, is both a strength and a weakness of the language. It is a strength because the expressivity lets you get a lot done with a single line of code. But it is also a weakness, because with that flexibility it is possible to introduce very subtle, strange-looking bugs if you are not familiar with the intricacies. For example, if you add a column vector to a row vector, you might expect a dimension mismatch or a type error, but you actually get back a matrix as the sum. There is an internal logic to these effects, but they can produce very hard-to-find bugs. Here are some tips and tricks to eliminate the strange-looking bugs from your code.

Rank 1 Arrays

Consider this line.

a = np.random.randn(5)

This creates five random Gaussian variables stored in an array a whose shape is the funny-looking (5,). This is called a rank 1 array in Python, and it is neither a row vector nor a column vector, which leads to slightly non-intuitive effects. For example, a.T ends up looking exactly the same as a, and if you compute np.dot(a, a.T), you might expect an outer product, a matrix, but you instead get back a single number.

a         = [ 1.78862847  0.43650985  0.09649747 -1.8634927  -0.2773882 ]
a.shape   = (5,)
a.T       = [ 1.78862847  0.43650985  0.09649747 -1.8634927  -0.2773882 ]
np.dot(a, a.T) = 6.948593697290846

Instead, commit every array you create to being either a column vector or a row vector.

a = np.random.randn(5, 1)   # (5,1) column vector
a = np.random.randn(1, 5)   # (1,5) row vector

With a = np.random.randn(5, 1), a is a \((5, 1)\) column vector, and now a.T really is a row vector. Notice one subtle difference in the printout. The row vector has two square brackets, because it is really a \(1 \times 5\) matrix, whereas the rank 1 array printed with one square bracket. And np.dot(a, a.T) now gives the outer product of the vector, a matrix, as expected.

a =
[[ 1.78862847]
 [ 0.43650985]
 [ 0.09649747]
 [-1.8634927 ]
 [-0.2773882 ]]
a.shape = (5, 1)
a.T = [[ 1.78862847  0.43650985  0.09649747 -1.8634927  -0.2773882 ]]
np.dot(a, a.T) =
[[ 3.19919182  0.78075395  0.17259812 -3.33309611 -0.49614444]
 [ 0.78075395  0.19054085  0.0421221  -0.81343292 -0.12108268]
 [ 0.17259812  0.0421221   0.00931176 -0.17982233 -0.02676726]
 [-3.33309611 -0.81343292 -0.17982233  3.47260506  0.51691089]
 [-0.49614444 -0.12108268 -0.02676726  0.51691089  0.07694421]]

Assertions and Reshape

If you are not entirely sure what the dimension of one of your vectors is, throw in an assertion statement.

assert(a.shape == (5, 1))

These assertions are really inexpensive to execute, and they also help serve as documentation for your code, so do not hesitate to throw them in whenever you feel like it. And if for some reason you do end up with a rank 1 array, you can reshape it into an explicit column or row vector.

a = a.reshape((5, 1))

The takeaways are

  1. To simplify your code, do not use rank 1 arrays. Always use either \(n \times 1\) matrices (column vectors) or \(1 \times n\) matrices (row vectors).
  2. Feel free to toss in plenty of assertion statements to double-check the dimensions of your matrices and arrays.
  3. Do not be shy about calling reshape to make sure your vectors are the dimension you need.

Eliminating rank 1 arrays does not actually restrict what you can express in code, and it removes a whole cause of bugs.

Review Questions

1. What is a rank 1 array, and why is it best avoided?

An array with a shape like (5,), produced for example by np.random.randn(5). It is neither a row vector nor a column vector and does not behave consistently as either. a.T looks the same as a, and np.dot(a, a.T) gives a number instead of the outer product matrix, which leads to subtle, hard-to-find bugs.

1. How do you tell a \(1 \times 5\) row vector from a rank 1 array in a printout?

The row vector prints with two square brackets (it is really a \(1 \times 5\) matrix), while the rank 1 array prints with one square bracket.

1. What three habits keep NumPy vector code free of these shape bugs?

Commit every array to being an explicit column vector (\(n \times 1\)) or row vector (\(1 \times n\)) instead of a rank 1 array, add cheap assertion statements like assert(a.shape == (5, 1)) that also document the code, and use reshape whenever you need to force an array into the intended dimensions.

1. Consider the NumPy array x = np.array([[[1], [2]], [[3], [4]]]). What is the shape of x?

  1. (2, 2, 1)
  2. (4,)
  3. (2, 2)
  4. (1, 2, 2)

a. Count the nesting from the outside in. The outer list has 2 elements, each of those is a list of 2 elements, and each innermost list holds 1 number. So the array has two rows, each row has 2 arrays of size \(1 \times 1\), giving shape (2, 2, 1).

Justification of the Logistic Regression Cost Function

This closing section is optional. It gives a quick justification for why we use the particular cost function we do for logistic regression.

To recap, the prediction is \(\hat{y} = \sigma(w^T x + b)\), and we interpret \(\hat{y}\) as \(P(y = 1 \mid x)\), the chance that \(y = 1\) for a given set of input features \(x\). Another way to say this is

  • If \(y = 1\), then \(p(y \mid x) = \hat{y}\).
  • If \(y = 0\), then \(p(y \mid x) = 1 - \hat{y}\), because if \(\hat{y}\) is the chance that \(y = 1\), then \(1 - \hat{y}\) is the chance that \(y = 0\).

These two equations define \(p(y \mid x)\) for the two cases, and since this is binary classification, \(y = 0\) and \(y = 1\) are the only two possible cases. We can summarize the two equations into a single one,

\[ p(y \mid x) = \hat{y}^{\, y} \, (1 - \hat{y})^{(1 - y)} \]

Here is why this one line works. If \(y = 1\), the first term is \(\hat{y}\) to the power of 1, and the second term is \((1 - \hat{y})\) to the power of \(1 - 1 = 0\). Anything to the power of 0 equals 1, so that term goes away, leaving \(p(y \mid x) = \hat{y}\), exactly what we wanted. If \(y = 0\), the first term is \(\hat{y}\) to the power of 0, which is 1, and the second term is \((1 - \hat{y})\) to the power of \(1 - 0 = 1\), leaving \(p(y \mid x) = 1 - \hat{y}\), again exactly what we wanted.

Now, because the log function is a strictly monotonically increasing function, maximizing \(\log p(y \mid x)\) gives the same result as maximizing \(p(y \mid x)\). Taking the log,

\[ \log p(y \mid x) = y \log \hat{y} + (1 - y) \log(1 - \hat{y}) = -\mathcal{L}(\hat{y}, y) \]

This is exactly the negative of the loss function we defined previously. The negative sign is there because when training a learning algorithm we want to make the probability large, whereas in logistic regression we express this as minimizing the loss. Minimizing the loss corresponds to maximizing the log of the probability.

Cost on the Whole Training Set

How about the overall cost on the entire training set of \(m\) examples? If we assume the training examples were drawn i.i.d. (independently and identically distributed), then the probability of all the labels in the training set is the product of the probabilities,

\[ P(\text{labels in training set}) = \prod_{i=1}^{m} p\big(y^{(i)} \mid x^{(i)}\big) \]

Maximum likelihood estimation is the principle in statistics of choosing the parameters that maximize the chance of your observations. It is introduced from scratch, including the same log trick used here, on the point estimation and maximum likelihood page of the statistics course. Maximizing this product is the same as maximizing its log, and the log of a product is the sum of the logs,

\[ \log P(\text{labels in training set}) = \sum_{i=1}^{m} \log p\big(y^{(i)} \mid x^{(i)}\big) = -\sum_{i=1}^{m} \mathcal{L}\big(\hat{y}^{(i)}, y^{(i)}\big) \]

So carrying out maximum likelihood estimation means maximizing this quantity, which is the same as minimizing \(\sum_{i=1}^{m} \mathcal{L}\big(\hat{y}^{(i)}, y^{(i)}\big)\), since we got rid of the minus sign by flipping from maximizing to minimizing. Finally, for convenience, to make sure our quantities are better scaled, we add an extra \(\frac{1}{m}\) scaling factor. That gives exactly the cost we had,

\[ J(w, b) = \frac{1}{m} \sum_{i=1}^{m} \mathcal{L}\big(\hat{y}^{(i)}, y^{(i)}\big) \]

To summarize, by minimizing the cost \(J(w, b)\) we are really carrying out maximum likelihood estimation with the logistic regression model, under the assumption that the training examples were drawn i.i.d. The log loss page in the calculus notes and the cost function for logistic regression page in the machine learning course tell the same story from different angles, and the statistics course applies the same machinery to a different model in MLE and linear regression.

Review Questions

1. Why does the single formula \(p(y \mid x) = \hat{y}^{\, y} (1 - \hat{y})^{(1 - y)}\) capture both cases of a binary label?

When \(y = 1\), the second factor has exponent 0 and disappears (anything to the power of 0 is 1), leaving \(\hat{y}\). When \(y = 0\), the first factor has exponent 0 and disappears, leaving \(1 - \hat{y}\). Those match the two defining equations \(p(y=1 \mid x) = \hat{y}\) and \(p(y=0 \mid x) = 1 - \hat{y}\).

1. Why is it valid to maximize \(\log p(y \mid x)\) instead of \(p(y \mid x)\)?

Because the log is a strictly monotonically increasing function, the parameters that maximize \(\log p(y \mid x)\) are the same ones that maximize \(p(y \mid x)\).

1. True or False: minimizing the loss corresponds with maximizing \(\log p(y \mid x)\).

  1. False
  2. True

b. True. The loss is exactly the negative of the log probability, \(\log p(y \mid x) = -\mathcal{L}(\hat{y}, y)\). Making the loss as small as possible therefore makes \(\log p(y \mid x)\) as large as possible. The negative sign is there because training conventionally minimizes a loss, while the probabilistic view wants the probability of the observed label to be large.

1. How does the i.i.d. assumption lead from per-example probabilities to the cost function \(J(w, b)\)?

Under i.i.d. sampling, the probability of all the labels is the product \(\prod_{i=1}^{m} p\big(y^{(i)} \mid x^{(i)}\big)\). Taking logs turns the product into a sum of \(\log p\big(y^{(i)} \mid x^{(i)}\big) = -\mathcal{L}\big(\hat{y}^{(i)}, y^{(i)}\big)\). Maximizing that sum equals minimizing \(\sum_i \mathcal{L}\), and adding the \(\frac{1}{m}\) scaling factor for convenience gives \(J(w, b)\).

1. What is the relationship between minimizing \(J(w, b)\) and maximum likelihood estimation?

They are the same thing. Minimizing the logistic regression cost \(J(w, b)\) carries out maximum likelihood estimation of the parameters, under the assumption that the training examples were drawn i.i.d.

Back to top