Vectorization

machine-learning
supervised-learning
regression
In multiple linear regression, you learned to write the model compactly as:
Published

July 6, 2026

What Is Vectorization?

In multiple linear regression, you learned to write the model compactly as:

\[ f_{\vec{w},b}(\vec{x}) = \vec{w} \cdot \vec{x} + b \]

Vectorization is the practice of writing code that takes advantage of this vector form. When you implement a learning algorithm with vectorization, your code becomes shorter and runs much faster.

Vectorized code also lets you use modern numerical linear algebra libraries (like NumPy in Python) and hardware such as GPUs (graphics processing units). GPUs were originally designed to speed up computer graphics, but they are also excellent at the kind of parallel math that vectorized machine learning code needs.

Review Questions

1. What are two benefits of writing vectorized code?

Vectorized code is shorter (easier to write and read) and runs much faster than equivalent code written with many separate operations or for loops.

1. What library does this course use for numerical linear algebra in Python?

NumPy. It is the most widely used numerical linear algebra library in Python and machine learning.

Vectors in Math and Python

Note

Vectors in math notation use an arrow: \(\vec{w}\), \(\vec{x}\). Scalars use plain symbols: \(w_j\), \(x_j\), \(b\). See Multiple Linear Regression for the full notation guide.

Here is a concrete example. Suppose \(n = 3\), so we have three features. The parameter vector \(\vec{w}\) and feature vector \(\vec{x}\) each contain three numbers:

\[ \vec{w} = \begin{bmatrix} w_1 & w_2 & w_3 \end{bmatrix}, \qquad \vec{x} = \begin{bmatrix} x_1 & x_2 & x_3 \end{bmatrix} \]

In linear algebra, indexing usually starts at 1: \(w_1\) is the first weight and \(x_1\) is the first feature.

In Python with NumPy, array indexing starts at 0:

import numpy as np

w = np.array([1.0, 2.0, 3.0])   # w_1, w_2, w_3 in math notation
x = np.array([4.0, 5.0, 6.0])   # x_1, x_2, x_3 in math notation
b = 0.5

# Math index 1  -> Python index 0
w[0]   # first weight (w_1)
w[1]   # second weight (w_2)
w[2]   # third weight (w_3)
Math notation Python (NumPy)
\(w_1\) w[0]
\(w_2\) w[1]
\(w_3\) w[2]
\(x_1\) x[0]
\(x_2\) x[1]
\(x_3\) x[2]

Many programming languages, including Python, count from 0 rather than 1. When you read math formulas and write Python code, keep this offset in mind.

Review Questions

1. In math notation, what is the first element of \(\vec{w}\) called?

\[ w_1 \]

Math indexing starts at 1.

1. In NumPy, how do you access the first element of w?

w[0]. Python array indexing starts at 0.

1. If there are three features, what is the last valid Python index for accessing elements of w?

\[ n = 3 \]

\[ 2 \]

The valid indices are:

\[ 0, \; 1, \; 2 \]

Without Vectorization

Manual Term-by-Term Calculation

One way to compute the prediction is to multiply each weight by its matching feature and add the results by hand:

f = w[0] * x[0] + w[1] * x[1] + w[2] * x[2] + b

This works when \(n = 3\), but what if \(n = 100\) or \(n = 100{,}000\)? You would have to type an enormous expression. It is inefficient for you to write and inefficient for the computer to run.

Using a For Loop

In math, you can use summation notation:

\[ f_{\vec{w},b}(\vec{x}) = \sum_{j=1}^{n} w_j x_j + b \]

For \(n = 3\), the index \(j\) runs through 1, 2, and 3.

In code, you can initialize f to zero and loop over all features:

f = 0
for j in range(n):
    f = f + w[j] * x[j]
f = f + b

Here range(n) makes j go from 0 to n - 1 (it does not include n itself). So for \(n = 3\), j takes values 0, 1, and 2, which correctly maps to \(w_1, w_2, w_3\) in math notation.

This loop version is a bit better than typing every term manually, but it still does not use vectorization and is not as efficient as it could be.

Manual: 32.5
Loop:   32.5

Review Questions

1. What does this summation compute?

\[ \sum_{j=1}^{n} w_j x_j + b \]

The model prediction:

\[ f_{\vec{w},b}(\vec{x}) = \sum_{j=1}^{n} w_j x_j + b \]

1. In Python, what values does j take when you write for j in range(n)?

\[ n = 3 \]

\[ j \in \{0, 1, 2\} \]

range(n) starts at 0 and stops before n, so it never includes 3.

1. Why is the manual term-by-term approach impractical when \(n\) is large?

You would need to write an extremely long expression with one term per feature. It is hard to maintain and slow to compute.

Vectorized Implementation

The math expression is a dot product plus a bias:

\[ f_{\vec{w},b}(\vec{x}) = \vec{w} \cdot \vec{x} + b \]

In NumPy, you can compute this in a single line:

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

np.dot(w, x) computes the dot product between the two vectors. The first np refers to the NumPy library; dot is the function that performs the dot product.

This one line replaces the long manual expression or the for loop. When \(n\) is large, np.dot runs much faster than either non-vectorized approach.

Vectorized: 32.5
Matches manual and loop? True

Review Questions

1. Which of the following is a vectorized implementation for computing a linear regression model’s prediction?

f = np.dot(w, x) + b
f = w[0] * x[0] + w[1] * x[1] + w[2] * x[2] + b
f = 0
for j in range(n):
    f = f + w[j] * x[j]
f = f + b

a. np.dot(w, x) + b is the vectorized implementation. NumPy’s dot function uses parallel hardware to efficiently compute:

\[ \vec{w} \cdot \vec{x} \]

in one call. Option (b) is correct math for three features but expands every term by hand (not vectorized). Option (c) uses a for loop and computes the same result sequentially.

1. What does f = np.dot(w, x) + b compute?

\[ f_{\vec{w},b}(\vec{x}) = \vec{w} \cdot \vec{x} + b \]

It computes the dot product of w and x, then adds the scalar bias b.

1. How many lines of core computation does the vectorized version need, compared to a for loop?

One line (np.dot(w, x) + b) instead of initializing a variable, looping over all indices, and adding b afterward.

Why Vectorization Is Faster

Vectorization gives you two distinct benefits:

  1. Shorter code. One line is easier to write, read, and share with others.
  2. Faster execution. np.dot runs much faster than manual term-by-term code or a for loop when \(n\) is large.

The speed comes from what happens behind the scenes. NumPy is designed to use parallel hardware on your computer. Whether you run code on a normal CPU or on a GPU, vectorized operations can perform many multiplications and additions at once instead of one at a time.

That is why the for loop (which updates f one step at a time) is slower: it processes each product sequentially. np.dot delegates the work to optimized linear algebra routines that take full advantage of your hardware.

Note

The non-vectorized manual version saves typing compared to writing hundreds of terms, but it is still not computationally efficient because it does not use vectorization. Only np.dot (or similar vectorized functions) gets the full speed benefit.

Review Questions

1. What are the two main benefits of vectorization?

Shorter, more readable code and much faster execution, especially when the number of features \(n\) is large.

1. Why does np.dot(w, x) run faster than a for loop when the number of features is large?

np.dot(w, x) uses optimized, parallel linear algebra routines that can perform many operations at once on CPU or GPU hardware. A for loop processes one multiplication at a time sequentially. The speed gap grows as \(n\) increases.

1. Does vectorization change the mathematical result of the computation?

No. Vectorized and non-vectorized code compute the same value. Vectorization changes how efficiently the computer performs the calculation, not the formula itself.

What Happens Behind the Scenes

Vectorization can feel like a magic trick the first time you see it: the same math, but dramatically faster. Here is what your computer is actually doing.

Sequential For Loop

Consider a for loop that updates one index at a time:

for j in range(16):
    f = f + w[j] * x[j]

If j runs from 0 to 15, the computer performs the operations one after another:

  • At time \(t_0\): compute using index 0
  • At time \(t_1\): compute using index 1
  • \(\ldots\)
  • At time \(t_{15}\): compute using index 15

Each step waits for the previous one to finish. There is no parallelism.

Parallel Vectorized Operations

np.dot(w, x) is implemented to take advantage of vectorized hardware. Instead of 16 separate steps, the computer can:

  1. Fetch all values in w and x
  2. Multiply every matching pair \((w_j, x_j)\) at the same time
  3. Add all 16 products together using specialized hardware, rather than 16 separate additions one after another

The result is the same number. The path to get there is much faster.

This matters most when you train on large datasets or models with many features, which is typical in modern machine learning. Vectorization is a key reason learning algorithms can scale to the data sizes used in practice today.

Review Questions

1. How does a for loop process these indices?

\[ j = 0, 1, \ldots, 15 \]

One index at a time, in sequence. Each step completes before the next begins:

\[ t_0, \; t_1, \; \ldots, \; t_{15} \]

1. What two operations does np.dot(w, x) perform in parallel behind the scenes?

First, multiply all corresponding pairs at once:

\[ w_1 x_1, \; w_2 x_2, \; \ldots, \; w_n x_n \]

Then add all of those products together efficiently.

1. When does vectorization make the biggest practical difference?

On large datasets and models with many features (large \(n\)), where sequential code would be too slow to be practical.

Vectorized Weight Updates

Vectorization helps not only when computing predictions, but also when updating parameters during gradient descent.

Suppose you have \(n = 16\) features and 16 weights \(w_1, \ldots, w_{16}\) (ignoring \(b\) for this example). After computing derivatives, you store the weight values in a NumPy array w and the derivative values in an array d.

Each weight update is:

\[ w_j \leftarrow w_j - \alpha \, d_j \]

Say the learning rate is \(\alpha = 0.1\). For all 16 weights, that means:

\[ w_1 \leftarrow w_1 - 0.1 \, d_1, \quad w_2 \leftarrow w_2 - 0.1 \, d_2, \quad \ldots, \quad w_{16} \leftarrow w_{16} - 0.1 \, d_{16} \]

Without Vectorization

You could write 16 separate lines, or use a for loop:

for j in range(16):
    w[j] = w[j] - 0.1 * d[j]

This updates one weight at a time, sequentially.

With Vectorization

The same update for all 16 weights fits in one line:

w = w - 0.1 * d

Behind the scenes, the computer subtracts \(0.1 \times d_j\) from every \(w_j\) in parallel and writes all 16 results back to w at once.

For loop: [0.95 1.9  2.98 3.92]
Vectorized: [0.95 1.9  2.98 3.92]
Equal? True

With only 16 features, the speed difference may be modest. With thousands of features and large training sets, vectorized updates can mean the difference between code that finishes in minutes versus code that runs for hours.

You can time these two implementations yourself in Python and see how much vectorization helps as the feature count grows.

In the next section, you will combine multiple linear regression with vectorization to implement gradient descent for multiple linear regression.

Review Questions

1. What update rule is applied to each weight \(w_j\)?

\[ w_j \leftarrow w_j - \alpha \, d_j \]

Each weight is reduced by the learning rate times the corresponding derivative. For example, when:

\[ \alpha = 0.1 \]

each update is:

\[ w_j \leftarrow w_j - 0.1 \, d_j \]

1. Which line of code is the vectorized weight update for all 16 parameters at once?

for j in range(16):
    w[j] = w[j] - 0.1 * d[j]
w = w - 0.1 * d

b. w = w - 0.1 * d updates the entire weight vector in one vectorized operation. Option (a) is correct math but runs sequentially in a for loop.

1. Does vectorization change the values of the weights after an update step?

No. Vectorized and for-loop code produce the same updated weights:

\[ w_j \leftarrow w_j - \alpha \, d_j \]

Vectorization only changes how quickly the computer performs the updates.

Lab: Python, NumPy and Vectorization

In this lab, you will practice the NumPy skills used throughout Course 1. You will create vectors and matrices, index and slice them, implement a dot product with a for loop, compare it to np.dot, and time the difference on large arrays.

Goals

In this lab, you will:

  • create NumPy 1-D arrays (vectors) and 2-D arrays (matrices)
  • practice indexing, slicing, and element-wise operations
  • implement my_dot with a for loop and compare it to np.dot
  • time vectorized vs loop dot products on large arrays
  • preview how training data shapes appear in later labs

Setup

NumPy is the standard numerical library for Python in machine learning. It is conventional to import it as np. We also import time to compare how long different implementations take.

import numpy as np
import time

Creating Vectors

NumPy arrays have a shape (dimensions) and a dtype (data type). A 1-D array with shape (n,) is how we represent a vector of \(n\) numbers in this course.

Data-creation routines often take a shape as the first argument. Below are several ways to create a length-4 vector:

a = np.zeros(4)
print(f"np.zeros(4) :   a = {a}, a shape = {a.shape}, a data type = {a.dtype}")

a = np.zeros((4,))
print(f"np.zeros(4,) :  a = {a}, a shape = {a.shape}, a data type = {a.dtype}")

a = np.random.random_sample(4)
print(f"np.random.random_sample(4): a = {a}, a shape = {a.shape}, a data type = {a.dtype}")
np.zeros(4) :   a = [0. 0. 0. 0.], a shape = (4,), a data type = float64
np.zeros(4,) :  a = [0. 0. 0. 0.], a shape = (4,), a data type = float64
np.random.random_sample(4): a = [0.55475831 0.85430005 0.83537495 0.06412351], a shape = (4,), a data type = float64

Some routines take a size argument instead of a shape tuple:

a = np.arange(4.)
print(f"np.arange(4.):     a = {a}, a shape = {a.shape}, a data type = {a.dtype}")

a = np.random.rand(4)
print(f"np.random.rand(4): a = {a}, a shape = {a.shape}, a data type = {a.dtype}")
np.arange(4.):     a = [0. 1. 2. 3.], a shape = (4,), a data type = float64
np.random.rand(4): a = [0.00632104 0.85053139 0.86341432 0.31146075], a shape = (4,), a data type = float64

You can also pass explicit values. Note how 5. makes the array use floating-point dtype:

a = np.array([5, 4, 3, 2])
print(f"np.array([5,4,3,2]):  a = {a},     a shape = {a.shape}, a data type = {a.dtype}")

a = np.array([5., 4, 3, 2])
print(f"np.array([5.,4,3,2]): a = {a}, a shape = {a.shape}, a data type = {a.dtype}")
np.array([5,4,3,2]):  a = [5 4 3 2],     a shape = (4,), a data type = int64
np.array([5.,4,3,2]): a = [5. 4. 3. 2.], a shape = (4,), a data type = float64

All of these create a 1-D vector with four elements. The shape (4,) means one dimension with length 4.

Indexing Vectors

NumPy indexing starts at 0. The element at math position \(x_2\) is a[2] in code. Accessing a single element returns a scalar (no shape).

a = np.arange(10)
print(a)

print(f"a[2].shape: {a[2].shape} a[2]  = {a[2]}, Accessing an element returns a scalar")

print(f"a[-1] = {a[-1]}")

try:
    c = a[10]
except Exception as e:
    print("The error message you'll see is:")
    print(e)
[0 1 2 3 4 5 6 7 8 9]
a[2].shape: () a[2]  = 2, Accessing an element returns a scalar
a[-1] = 9
The error message you'll see is:
index 10 is out of bounds for axis 0 with size 10

Slicing Vectors

Slicing uses start:stop:step to grab a subset of elements:

a = np.arange(10)
print(f"a         = {a}")

c = a[2:7:1];  print("a[2:7:1] = ", c)
c = a[2:7:2];  print("a[2:7:2] = ", c)
c = a[3:];     print("a[3:]    = ", c)
c = a[:3];     print("a[:3]    = ", c)
c = a[:];      print("a[:]     = ", c)
a         = [0 1 2 3 4 5 6 7 8 9]
a[2:7:1] =  [2 3 4 5 6]
a[2:7:2] =  [2 4 6]
a[3:]    =  [3 4 5 6 7 8 9]
a[:3]    =  [0 1 2]
a[:]     =  [0 1 2 3 4 5 6 7 8 9]

Operations on One Vector

NumPy can negate elements, sum them, take the mean, or square them element-wise:

a = np.array([1, 2, 3, 4])
print(f"a             : {a}")

b = -a
print(f"b = -a        : {b}")

b = np.sum(a)
print(f"b = np.sum(a) : {b}")

b = np.mean(a)
print(f"b = np.mean(a): {b}")

b = a**2
print(f"b = a**2      : {b}")
a             : [1 2 3 4]
b = -a        : [-1 -2 -3 -4]
b = np.sum(a) : 10
b = np.mean(a): 2.5
b = a**2      : [ 1  4  9 16]

Element-Wise Operations on Two Vectors

Arithmetic operators apply element by element. For vectors \(\vec{a}\) and \(\vec{b}\) of the same length:

\[ c_i = a_i + b_i \]

a = np.array([1, 2, 3, 4])
b = np.array([-1, -2, 3, 4])
print(f"Binary operators work element wise: {a + b}")
Binary operators work element wise: [0 0 6 8]

The vectors must be the same size:

c = np.array([1, 2])
try:
    d = a + c
except Exception as e:
    print("The error message you'll see is:")
    print(e)
The error message you'll see is:
operands could not be broadcast together with shapes (4,) (2,)

Scalar Times Vector

A scalar multiplies every element:

a = np.array([1, 2, 3, 4])
b = 5 * a
print(f"b = 5 * a : {b}")
b = 5 * a : [ 5 10 15 20]

Dot Product: For Loop vs np.dot

The dot product multiplies corresponding entries and sums the results:

\[ \vec{a} \cdot \vec{b} = \sum_{i=0}^{n-1} a_i b_i \]

First, implement it with a for loop:

def my_dot(a, b):
    """
    Compute the dot product of two vectors

    Args:
      a (ndarray (n,)):  input vector
      b (ndarray (n,)):  input vector with same dimension as a

    Returns:
      x (scalar):
    """
    x = 0
    for i in range(a.shape[0]):
        x = x + a[i] * b[i]
    return x

Test on small vectors:

a = np.array([1, 2, 3, 4])
b = np.array([-1, 4, 3, 2])
print(f"my_dot(a, b) = {my_dot(a, b)}")
my_dot(a, b) = 24

The result should be a scalar. Compare with NumPy:

a = np.array([1, 2, 3, 4])
b = np.array([-1, 4, 3, 2])
c = np.dot(a, b)
print(f"NumPy 1-D np.dot(a, b) = {c}, np.dot(a, b).shape = {c.shape}")
c = np.dot(b, a)
print(f"NumPy 1-D np.dot(b, a) = {c}, np.dot(b, a).shape = {c.shape}")
NumPy 1-D np.dot(a, b) = 24, np.dot(a, b).shape = ()
NumPy 1-D np.dot(b, a) = 24, np.dot(b, a).shape = ()

Timing: Vectorized vs Loop

On very large vectors, np.dot is much faster because NumPy uses optimized, parallel hardware (SIMD on CPUs, and GPUs when available).

np.random.seed(1)
a = np.random.rand(10000000)
b = np.random.rand(10000000)

tic = time.time()
c = np.dot(a, b)
toc = time.time()
print(f"np.dot(a, b) =  {c:.4f}")
print(f"Vectorized version duration: {1000*(toc-tic):.4f} ms")

tic = time.time()
c = my_dot(a, b)
toc = time.time()
print(f"my_dot(a, b) =  {c:.4f}")
print(f"loop version duration: {1000*(toc-tic):.4f} ms")

del a; del b
np.dot(a, b) =  2501072.5817
Vectorized version duration: 3.1238 ms
my_dot(a, b) =  2501072.5817
loop version duration: 879.0026 ms

You should see the vectorized version finish in a fraction of the time. The exact numbers depend on your machine.

Shapes in Course 1

Later labs store training data in a 2-D array X_train with shape (m, n) (\(m\) examples, \(n\) features). Each row X[i] is a 1-D vector of shape (n,). Here is a small preview:

X = np.array([[1], [2], [3], [4]])
w = np.array([2])
c = np.dot(X[1], w)

print(f"X[1] has shape {X[1].shape}")
print(f"w has shape {w.shape}")
print(f"c has shape {c.shape}")
X[1] has shape (1,)
w has shape (1,)
c has shape ()

Paying attention to .shape prevents silent bugs when you combine vectors and matrices.

Matrices in NumPy

A matrix is a 2-D array with shape (m, n) (rows, columns). In Course 1, X_train will have shape (m, n).

Create 2-D arrays with a shape tuple:

a = np.zeros((1, 5))
print(f"a shape = {a.shape}, a = {a}")

a = np.zeros((2, 1))
print(f"a shape = {a.shape}, a = {a}")

a = np.random.random_sample((1, 1))
print(f"a shape = {a.shape}, a = {a}")
a shape = (1, 5), a = [[0. 0. 0. 0. 0.]]
a shape = (2, 1), a = [[0.]
 [0.]]
a shape = (1, 1), a = [[0.44236513]]

Or specify values directly (each inner list is a row):

a = np.array([[5], [4], [3]])
print(f" a shape = {a.shape}, np.array: a = {a}")

a = np.array([[5],
              [4],
              [3]])
print(f" a shape = {a.shape}, np.array: a = {a}")
 a shape = (3, 1), np.array: a = [[5]
 [4]
 [3]]
 a shape = (3, 1), np.array: a = [[5]
 [4]
 [3]]

Matrix Indexing

Use [row, column]. Accessing one index like a[2] returns the entire row as a 1-D vector:

a = np.arange(6).reshape(-1, 2)
print(f"a.shape: {a.shape}, \na= {a}")

print(f"\na[2,0].shape:   {a[2, 0].shape}, a[2,0] = {a[2, 0]},     type(a[2,0]) = {type(a[2, 0])} Accessing an element returns a scalar\n")

print(f"a[2].shape:   {a[2].shape}, a[2]   = {a[2]}, type(a[2])   = {type(a[2])}")
a.shape: (3, 2), 
a= [[0 1]
 [2 3]
 [4 5]]

a[2,0].shape:   (), a[2,0] = 4,     type(a[2,0]) = <class 'numpy.int64'> Accessing an element returns a scalar

a[2].shape:   (2,), a[2]   = [4 5], type(a[2])   = <class 'numpy.ndarray'>

reshape(-1, 2) turns a length-6 vector into 3 rows and 2 columns. The -1 tells NumPy to infer the number of rows.

Matrix Slicing

a = np.arange(20).reshape(-1, 10)
print(f"a = \n{a}")

print("a[0, 2:7:1] = ", a[0, 2:7:1], ",  a[0, 2:7:1].shape =", a[0, 2:7:1].shape, "a 1-D array")

print("a[:, 2:7:1] = \n", a[:, 2:7:1], ",  a[:, 2:7:1].shape =", a[:, 2:7:1].shape, "a 2-D array")

print("a[:,:] = \n", a[:,:], ",  a[:,:].shape =", a[:,:].shape)

print("a[1,:] = ", a[1,:], ",  a[1,:].shape =", a[1,:].shape, "a 1-D array")
print("a[1]   = ", a[1],   ",  a[1].shape   =", a[1].shape, "a 1-D array")
a = 
[[ 0  1  2  3  4  5  6  7  8  9]
 [10 11 12 13 14 15 16 17 18 19]]
a[0, 2:7:1] =  [2 3 4 5 6] ,  a[0, 2:7:1].shape = (5,) a 1-D array
a[:, 2:7:1] = 
 [[ 2  3  4  5  6]
 [12 13 14 15 16]] ,  a[:, 2:7:1].shape = (2, 5) a 2-D array
a[:,:] = 
 [[ 0  1  2  3  4  5  6  7  8  9]
 [10 11 12 13 14 15 16 17 18 19]] ,  a[:,:].shape = (2, 10)
a[1,:] =  [10 11 12 13 14 15 16 17 18 19] ,  a[1,:].shape = (10,) a 1-D array
a[1]   =  [10 11 12 13 14 15 16 17 18 19] ,  a[1].shape   = (10,) a 1-D array

a[1,:] and a[1] both return row 1 as a 1-D array. This pattern appears often when you pull one training example from X_train.

Summary

In this lab you:

  • created NumPy vectors and matrices and inspected .shape and .dtype
  • practiced indexing, slicing, and element-wise operations
  • implemented my_dot and verified it matches np.dot on small inputs
  • timed np.dot vs a for loop on 10 million elements
  • previewed how row vectors from a matrix feed into dot products in later labs

Keep these examples handy as a NumPy reference while you work through the rest of the course.

Review Questions

1. What does the shape (4,) mean for a NumPy array?

It is a 1-D array (vector) with 4 elements, indexed from 0 to 3.

1. What is the difference between a[2, 0] and a[2] when a is a 2-D matrix?

a[2, 0] returns a single scalar (row 2, column 0). a[2] returns the entire row 2 as a 1-D vector.

1. What does my_dot(a, b) compute?

\[ \vec{a} \cdot \vec{b} = \sum_{i=0}^{n-1} a_i b_i \]

The dot product: multiply each pair \(a_i b_i\) and sum all terms. It returns a scalar.

1. Why is np.dot(a, b) faster than my_dot(a, b) on very large vectors?

NumPy uses optimized, parallel linear algebra routines (SIMD on the CPU) that perform many multiplications and additions at once. my_dot runs one multiplication per loop iteration.