Why is lists faster than numpy?

Why in this example are lists faster than numpy? numpy is usually considered faster than using lists? But here lists seem to be over 2x faster?

$ python main.py 
using_list: 1.8105032444000244 seconds
using_numpy: 4.414405345916748 seconds

# main.py
import numpy as np
import time

def using_list(N):
    x = [1] * N
    y = [0] * N
    for i in range(1, N):
        y[i] = x[i] + 0.5 * (y[i-1] - x[i])

def using_numpy(N):
    x = np.ones(N)
    y = np.zeros(N)
    for i in range(1, N):
        y[i] = x[i] + 0.5 * (y[i-1] - x[i])

def main():
    LOOPS = 100
    N = 100000
    for func in [using_list, using_numpy]:
        t = time.perf_counter()
        for i in range(LOOPS):
            func(N)
        print(f"{func.__name__}: {time.perf_counter() - t} seconds")

if __name__=="__main__":
    main()

Using Python 3.10.

For reference, I have found that using @njit from numba will dramatically improve my numpy code. But I was wondering why my plain numpy version is so much slower than lists? And is there a way to make my plain numpy version faster without using numba.

using_numba: 0.6930002999997669 seconds

from numba import njit

@njit
def using_numba(N):
    x = np.ones(N)
    y = np.zeros(N)
    for i in range(1, N):
        y[i] = x[i] + 0.5 * (y[i-1] - x[i])

Iterating or using indices to access numpy arrays element by element is inefficient (please refer to Why is for loop on native python list faster than for loop on numpy array). Consider using numpy's built-in functions or methods to accelerate the loop, which can significantly improve the performance on large arrays:

def mechanic_using_numpy(N):
    y = np.empty(N, float)
    y[:1] = 1.0
    y[1:] = 0.5
    y.cumprod(out=y)
    np.subtract(1, y, out=y)
    return y

>>> mechanic_using_numpy(10)
array([0.        , 0.5       , 0.75      , 0.875     , 0.9375    ,
       0.96875   , 0.984375  , 0.9921875 , 0.99609375, 0.99804688])
>>> np.array(using_list(10))
array([0.        , 0.5       , 0.75      , 0.875     , 0.9375    ,
       0.96875   , 0.984375  , 0.9921875 , 0.99609375, 0.99804688])

Benchmark with N = 100000:

In [_]: %timeit using_numpy(100000)
44.6 ms ± 227 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [_]: %timeit using_list(100000)
14.3 ms ± 38.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [_]: %timeit mechanic_using_numpy(100000)
286 µs ± 875 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In order to get the benefits from numpy, you need to approach the calculation differently. Think of it as a global operation on all values rather than an iterative calculation:

For example, this function gives the same result but runs up to 25 times faster than using a list:

def using_np(N):
    return 1 - 0.5**(np.arange(N))

If you need more speed, you can combine simpler operations (e.g. addition & multiplication instead of exponentiation):

def using_np2(N):
    y      = np.ones(N)
    y[0]   = 0
    y[1:] -= np.cumprod(y[1:]/2)
    return y

This one is more than 40 times faster than using a list

Numpy is fast when it can implicitly parallelize operations. If you didn't have the y[i-1] -> y[i] data dependency, you could replace a loop like

for i in range(N):
    y[i] = 0.5*x[i]

with

y = 0.5*x

Iterating over the elements of an array is no faster than iterating over the elements of a list, and is in fact slightly slower due to the need to convert back and forth between Numpy value types and Python native types.

For example, a list of float values is just a list of references. y[i] is a reference to a float value that exists outside the list. An array with dtype of float64, however, is just a single binary blob that represents a collection of floating-point values. An operation like y[i] has to extract the proper bits then use them to create a float value to return. A more complicated expression like

y[i] = x[i] + 0.5 * (y[i-1] - x[i])

needs to extract three such values from x and y, do the same arithmetic you would do with the list approach, then convert the result back to a byte sequence to insert back into y.

With y = 0.5*x, there are no conversions of byte sequences to float values, use of float.__mul__, or conversions back to bytes sequences. All arithmetic is done internally in Numpy's optimized C code.