This article takes an in-depth look at the performance bottlenecks that occur when subtracting NumPy arrays and lists. The core reasons are the overhead of NumPy's internal iterators handling small-sized broadcast arrays, and the type promotion caused by the implicit conversion of Python floating point numbers to `np.float64`. By analyzing the performance differences of different implementation methods, the article reveals the key impact of data type, broadcast mechanism and memory layout on NumPy operating efficiency, and provides optimization solutions.

NumPy's performance is critical when working with large multi-dimensional arrays such as image data. However, sometimes seemingly simple operations can lead to unexpected performance differences. For example, when a NumPy array of shape 4000x4000x3 needs to subtract a set of channel values, two different implementations may result in a performance gap of tens of times.

Consider the following two subtraction implementations:

Implementation 1: Subtract the list directly

 import time
import numpy as np

image = np.random.rand(4000, 4000, 3).astype("float32")
values ??= [0.43, 0.44, 0.45]

st = time.time()
image -= values
et = time.time()
print("Realize 1 (directly subtract the list)", et - st)

Implementation 2: Loop to subtract list elements channel by channel

 import time
import numpy as np

image = np.random.rand(4000, 4000, 3).astype("float32")
values ??= [0.43, 0.44, 0.45]

st = time.time()
for i in range(3):
    image[..., i] -= values[i]
et = time.time()
print("Implementation 2 (loop channel by channel subtraction)", et - st)

In actual testing, implementation 2 ran about 20 times faster than implementation 1. This significant performance difference is no accident. It involves NumPy's internal mechanisms, including broadcasting, data type conversion, and memory access patterns.

Performance bottleneck analysis

The main reasons why the performance of implementation 1 is much lower than that of implementation 2 are as follows:

1. NumPy internal iterators and broadcast overhead

In order to achieve versatility and support advanced features such as broadcasting, NumPy introduces an internal iterator mechanism. When doing image -= values ??, NumPy will try to convert values ??(a Python list) into a NumPy array and broadcast it to the shape of image . For small arrays like values, NumPy's internal iterators introduce significant overhead when handling broadcasts. This is because the iterator needs to iterate over this small array many times to match the dimensions of the large array.

In addition, since the size of the values ??array is very small (for example, the shape is (3,)), it cannot even completely fill the SIMD (Single Instruction Multiple Data) register of mainstream CPUs, resulting in the inability to fully utilize the parallel computing advantages brought by SIMD instructions.

To verify this, we can observe performance changes by changing the size of the broadcast array. When the size of the broadcast array gradually increases, the performance will first improve because the iterator overhead is relatively reduced, but when the array is too large to fit completely into the CPU cache, the performance will decrease due to memory access delays.

 import numpy as np

# Assume that image has been defined as np.random.rand(4000, 4000, 3).astype("float32")
# For testing, we create a copy to avoid in-place modifications affecting subsequent tests test_image = np.random.rand(4000, 4000, 3).astype("float32")
values ??= [0.43, 0.44, 0.45]

# Flatten the image, and then test broadcast arrays of different sizes view = test_image.reshape(-1, 3)
print("Test np.tile(values, 1)")
%time view -= np.tile(values, 1)

view = test_image.reshape(-1, 6)
print("Test np.tile(values, 2)")
%time view -= np.tile(values, 2)

view = test_image.reshape(-1, 384)
print("Test np.tile(values, 128)")
%time view -= np.tile(values, 128)

view = test_image.reshape(-1, 3*4000)
print("Test np.tile(values, 4000)")
%time view -= np.tile(values, 4000)

It can be seen from the above experimental results that when the broadcast array size (generated through np.tile) gradually increases, the operation speed will speed up until a certain critical point (usually the array exceeds the CPU cache), and the performance will decrease again. This shows that, within certain limits, reducing the relative cost of NumPy iterators can improve performance.

2. Data types and implicit conversion

Another key issue is data type. In image -= values, values ??is a list of Python floating point numbers. When NumPy performs operations, it will be implicitly converted to a 1D array of type np.float64. The image array is of np.float32 type. According to NumPy's Type Promotion rules, in order to avoid precision loss, the entire subtraction operation will be performed with the precision of np.float64.

The np.float64 operation is generally slower than the np.float32 operation because it needs to process twice the amount of data and may not take full advantage of the hardware's float32 optimizations. This unnecessary type conversion and high-precision operations significantly reduce the performance of Implementation 1.

We can avoid this problem by explicitly specifying the data type of values:

 import numpy as np

# Assume that image has been defined as np.random.rand(4000, 4000, 3).astype("float32")
test_image = np.random.rand(4000, 4000, 3).astype("float32")
values_np_float32 = np.array([0.43, 0.44, 0.45], dtype=np.float32)

view = test_image.reshape(-1, 3)
print("Use np.float32 array for broadcasting")
%time view -= np.tile(values_np_float32, 1) # Performance will be significantly improved

After explicitly converting values ??to the np.float32 type, the performance of the broadcast operation will be greatly improved because the overhead of type conversion from float64 to float32 is avoided.

Why is implementation 2 faster?

Implementation 2 (loop-by-channel subtraction) is faster because it avoids the two main problems mentioned above:

1. No broadcast overhead: In image[..., i] -= values[i] , values[i] is a Python floating point scalar. When NumPy handles scalar and array operations, its internal mechanism is different from broadcast arrays. It converts the scalar directly to the same np.float32 type as the array element (for performance optimization), and then performs element-wise subtraction, avoiding the iterator overhead of a small broadcast array.
2. Data type consistency: Since image is np.float32, values[i] will also be efficiently converted to np.float32. The entire operation is performed in the np.float32 domain, avoiding the performance loss caused by np.float64.

Implementation 2 is not perfect, however. It requires traversing the entire image array 3 times (once for each channel), which means the entire array needs to be read from DRAM (main memory) and written 3 times, which is not very efficient in memory-intensive operations.

Optimization plan

Combining the above analysis, we can build a more optimized solution that avoids broadcast overhead and type conversion issues while reducing the number of memory accesses:

 import time
import numpy as np

image = np.random.rand(4000, 4000, 3).astype("float32")
values ??= [0.43, 0.44, 0.45]

st = time.time()
# Convert the values ??to an np.float32 array and use np.tile to repeat so that its shape matches the last dimension of the image # Then reshape to the form of (N, 3) and perform a broadcast operation with image.reshape(-1, 3) # Note: In order to broadcast with the shape of the original image here, more refined shape processing is required.
# A more direct way is to construct a shape (1, 1, 3) that can be broadcast directly
optimized_values ??= np.array(values, dtype=np.float32).reshape(1, 1, 3)
image -= optimized_values
et = time.time()
print("Optimized implementation (broadcast np.float32 array)", et - st)

In this optimized implementation:

• np.array(values, dtype=np.float32) ensures that the operation is performed with np.float32 precision.
• .reshape(1, 1, 3) changes the shape of the values ??array to (1, 1, 3), allowing it to be broadcast directly with the image array (4000, 4000, 3), and NumPy is generally more efficient in handling such broadcasts because it avoids frequent iterations of small arrays.

Memory layout considerations

In addition to the above factors, the memory layout of the array also has a significant impact on performance. For multi-channel image data, a common layout is (height, width, channels). However, this layout may not be optimal for NumPy and SIMD operations.

Usually, placing the channel dimension first, that is, using a (channels, height, width) layout, can bring better performance. This is because this layout is contiguous in memory, allowing SIMD instructions to process data more efficiently and aiding CPU cache utilization. For some operations, a (height, channels, width) layout might also be a good compromise.

For example, if you frequently need to operate on all channels, placing the channel dimension first can make data access more sequential, thereby improving cache hit ratio and SIMD parallelism.

 # Original layout (H, W, C)
image_hwc = np.random.rand(4000, 4000, 3).astype("float32")

# Convert to (C, H, W) layout image_chw = image_hwc.transpose(2, 0, 1)

# It may be more efficient to operate under the (C, H, W) layout values_chw = np.array(values, dtype=np.float32).reshape(3, 1, 1)
# image_chw -= values_chw # Example operation

Summary and best practices

To optimize the performance of NumPy array operations, especially when broadcasting and data type conversions are involved, the following points should be kept in mind:

1. Explicitly manage data types: Always explicitly specify the data type of a NumPy array (such as np.float32) to avoid the performance overhead caused by implicit conversion of Python lists or default floating point numbers to np.float64.
2. Use broadcasts with caution: avoid broadcasting on very small arrays, which can cause significant overhead on NumPy's internal iterators. If broadcasting is required, try to construct a NumPy array that can be broadcast efficiently (e.g., use reshape to adjust the shape to match the dimensions).
3. Optimize memory access patterns: For large multi-dimensional arrays, consider their memory layout. Put the most frequently iterated dimension at the end, or adjust it to (channels, height, width) and other layouts that are more suitable for SIMD and caching according to specific operations.
4. Avoid unnecessary loops: Take advantage of NumPy's vectorized operations whenever possible and avoid explicit loops at the Python level. Although the loop implementation is faster in this example, that's because it avoids the pitfalls of broadcasting and type conversion. Pure NumPy's vectorized operations are often the fastest when both data types and broadcasting are optimized.

By deeply understanding the internal mechanisms of NumPy, we can write more efficient and robust code, thereby fully utilizing its powerful numerical computing capabilities.

The above is the detailed content of Optimizing NumPy array subtraction operation: in-depth understanding of broadcast mechanism and data type impact. For more information, please follow other related articles on the PHP Chinese website!