NumPy - np.add, np.subtract, np.multiply, np.divide

Key Insights

NumPy’s arithmetic functions (np.add, np.subtract, np.multiply, np.divide) provide vectorized operations that are 10-100x faster than Python loops for numerical computations
These functions support broadcasting, allowing operations between arrays of different shapes following specific rules, eliminating the need for manual array reshaping
Understanding when to use ufuncs versus operators (+, -, *, /) matters for performance optimization and handling special cases like integer division and error handling

Basic Arithmetic Operations

NumPy’s core arithmetic functions operate element-wise on arrays. While Python operators work identically for most cases, the explicit functions offer additional parameters for advanced control.

import numpy as np

a = np.array([10, 20, 30, 40])
b = np.array([1, 2, 3, 4])

# Addition
result_add = np.add(a, b)
print(result_add)  # [11 22 33 44]

# Subtraction
result_sub = np.subtract(a, b)
print(result_sub)  # [ 9 18 27 36]

# Multiplication
result_mul = np.multiply(a, b)
print(result_mul)  # [ 10  40  90 160]

# Division
result_div = np.divide(a, b)
print(result_div)  # [10. 20. 10. 10.]

The operator equivalents produce identical results:

print(a + b)  # Same as np.add(a, b)
print(a - b)  # Same as np.subtract(a, b)
print(a * b)  # Same as np.multiply(a, b)
print(a / b)  # Same as np.divide(a, b)

Broadcasting Mechanics

Broadcasting allows arithmetic operations between arrays of different shapes without explicit replication. NumPy automatically expands dimensions according to specific rules.

# Scalar broadcasting
arr = np.array([[1, 2, 3],
                [4, 5, 6]])

result = np.add(arr, 10)
print(result)
# [[11 12 13]
#  [14 15 16]]

# 1D array broadcasting across rows
row_vector = np.array([10, 20, 30])
result = np.multiply(arr, row_vector)
print(result)
# [[10 40 90]
#  [40 100 180]]

# Column vector broadcasting
col_vector = np.array([[10], [20]])
result = np.add(arr, col_vector)
print(result)
# [[11 12 13]
#  [24 25 26]]

Broadcasting follows these rules:

Compare dimensions from right to left
Dimensions are compatible if they’re equal or one of them is 1
Missing dimensions are treated as 1

# Shape compatibility examples
a = np.ones((3, 4, 5))
b = np.ones((4, 5))     # Compatible: (4,5) broadcasts to (3,4,5)
c = np.ones((5,))       # Compatible: (5,) broadcasts to (3,4,5)
d = np.ones((3, 1, 5))  # Compatible: middle dimension broadcasts

result1 = np.add(a, b)  # Works
result2 = np.add(a, c)  # Works
result3 = np.add(a, d)  # Works

# Incompatible shapes
e = np.ones((3, 4))
# np.add(a, e)  # ValueError: operands could not be broadcast together

Output Arrays and In-Place Operations

The out parameter allows you to specify a pre-allocated array for results, avoiding memory allocation overhead in tight loops.

a = np.array([1.0, 2.0, 3.0, 4.0])
b = np.array([0.5, 1.5, 2.5, 3.5])

# Pre-allocate output array
output = np.empty(4)
np.add(a, b, out=output)
print(output)  # [1.5 3.5 5.5 7.5]

# In-place operation (modifies first array)
np.add(a, b, out=a)
print(a)  # [1.5 3.5 5.5 7.5]

# Performance comparison for large arrays
large_a = np.random.rand(1000000)
large_b = np.random.rand(1000000)
result_array = np.empty(1000000)

# Reusing output array in loops
for i in range(100):
    np.multiply(large_a, large_b, out=result_array)
    # Process result_array without creating new arrays

Division Variants and Integer Handling

NumPy provides multiple division functions for different use cases, particularly important when working with integer arrays.

a = np.array([10, 15, 20, 25])
b = np.array([3, 4, 6, 7])

# True division (always returns float)
true_div = np.divide(a, b)
print(true_div)  # [3.33333333 3.75 3.33333333 3.57142857]

# Floor division (rounds toward negative infinity)
floor_div = np.floor_divide(a, b)
print(floor_div)  # [3 3 3 3]

# Modulo operation
remainder = np.mod(a, b)
print(remainder)  # [1 3 2 4]

# Combined quotient and remainder
quotient, remainder = np.divmod(a, b)
print(quotient)   # [3 3 3 3]
print(remainder)  # [1 3 2 4]

Integer division behavior differs from Python 2 to Python 3. NumPy’s divide always performs true division:

int_a = np.array([10, 15], dtype=np.int32)
int_b = np.array([3, 4], dtype=np.int32)

# np.divide always returns float
result = np.divide(int_a, int_b)
print(result.dtype)  # float64
print(result)        # [3.33333333 3.75]

# Operator / also performs true division
result2 = int_a / int_b
print(result2)  # [3.33333333 3.75]

# Use // for floor division
result3 = int_a // int_b
print(result3)  # [3 3]

Error Handling and Special Values

Control how NumPy handles division by zero, overflow, and invalid operations using np.errstate or function parameters.

a = np.array([1.0, 0.0, -1.0])
b = np.array([0.0, 0.0, 0.0])

# Default behavior: warning
result = np.divide(a, b)
print(result)  # [inf nan -inf]

# Suppress warnings
with np.errstate(divide='ignore', invalid='ignore'):
    result = np.divide(a, b)
    print(result)  # [inf nan -inf]

# Raise exception
try:
    with np.errstate(divide='raise'):
        result = np.divide(a, b)
except FloatingPointError as e:
    print(f"Error caught: {e}")

# Using where parameter to avoid division by zero
a = np.array([10.0, 20.0, 30.0, 40.0])
b = np.array([2.0, 0.0, 5.0, 0.0])
result = np.divide(a, b, where=b!=0)
print(result)  # Undefined values where b==0

Performance Optimization

NumPy’s ufuncs are implemented in C and vectorized, providing significant performance advantages over Python loops.

import time

# Python loop approach
def python_loop_add(a, b):
    result = []
    for i in range(len(a)):
        result.append(a[i] + b[i])
    return result

# NumPy approach
size = 1000000
a_list = list(range(size))
b_list = list(range(size))
a_np = np.arange(size)
b_np = np.arange(size)

# Benchmark Python loop
start = time.time()
result_py = python_loop_add(a_list, b_list)
py_time = time.time() - start

# Benchmark NumPy
start = time.time()
result_np = np.add(a_np, b_np)
np_time = time.time() - start

print(f"Python loop: {py_time:.4f}s")
print(f"NumPy: {np_time:.4f}s")
print(f"Speedup: {py_time/np_time:.1f}x")

For complex expressions, consider using np.einsum or numexpr for additional optimization:

# Multiple operations
a = np.random.rand(1000, 1000)
b = np.random.rand(1000, 1000)
c = np.random.rand(1000, 1000)

# Standard approach
result1 = np.add(np.multiply(a, b), c)

# Equivalent with operators
result2 = a * b + c

# Pre-allocate for memory efficiency
temp = np.empty_like(a)
output = np.empty_like(a)
np.multiply(a, b, out=temp)
np.add(temp, c, out=output)

Multi-Dimensional Array Operations

Arithmetic functions work seamlessly with multi-dimensional arrays, respecting shape compatibility.

# 3D array operations
tensor_a = np.random.rand(10, 20, 30)
tensor_b = np.random.rand(10, 20, 30)

result = np.add(tensor_a, tensor_b)
print(result.shape)  # (10, 20, 30)

# Broadcasting with 3D arrays
matrix = np.random.rand(20, 30)
result = np.multiply(tensor_a, matrix)  # Broadcasts across first dimension
print(result.shape)  # (10, 20, 30)

# Axis-specific operations with broadcasting
weights = np.random.rand(10, 1, 1)
weighted = np.multiply(tensor_a, weights)
print(weighted.shape)  # (10, 20, 30)

These arithmetic functions form the foundation of NumPy’s computational capabilities. Understanding their broadcasting rules, performance characteristics, and special parameters enables efficient numerical computing for data science, machine learning, and scientific applications.