Skip to main content Brad's PyNotes

PEP 465: Matrix Multiplication Operator (@) - Scientific Computing Revolution

TL;DR

PEP 465 introduced the @ operator for matrix multiplication, solving a major readability problem in scientific Python code by providing a dedicated operator that makes mathematical formulas translate directly into code.

Interesting!

The @ operator was so needed that it’s one of the few times Python added a completely new binary operator! It solved the conflict between element-wise (*) and matrix multiplication that plagued numerical computing for years.

The Problem Before @

Competing Uses for * Operator

python code snippet start

import numpy as np

# Before PEP 465 - confusing and inconsistent
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Element-wise multiplication (Hadamard product)
hadamard = A * B
print("Element-wise:", hadamard)
# [[5, 12], [21, 32]]

# Matrix multiplication - had to use function calls
matrix_mult = np.dot(A, B)  # Or A.dot(B)
print("Matrix multiplication:", matrix_mult)
# [[19, 22], [43, 50]]

# This made translating math formulas very difficult!
# Mathematical: C = A × B + D × E
# Python: C = np.dot(A, B) + np.dot(D, E)  # Ugly!

python code snippet end

Real-World Pain Points

python code snippet start

# Complex scientific formulas were unreadable
# Mathematical formula: y = X^T × X × β - r
# Before @ operator:
y = np.dot(np.dot(X.T, X), beta) - r

# Compare to what mathematicians write:
# y = X.T @ X @ beta - r  # Much cleaner!

python code snippet end

The @ Operator Solution

Clean Matrix Operations

python code snippet start

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = np.array([[1, 0], [0, 1]])

# Matrix multiplication with @
result = A @ B
print("A @ B =", result)
# [[19, 22], [43, 50]]

# Chain operations naturally
result = A @ B @ C
print("A @ B @ C =", result)
# [[19, 22], [43, 50]]

# Mixed operations are clear
result = A @ B + C * 2  # Matrix mult + element-wise
print("Mixed:", result)
# [[21, 22], [43, 52]]

python code snippet end

Mathematical Formula Translation

python code snippet start

# Linear regression: β = (X^T × X)^(-1) × X^T × y
# Before @:
beta_old = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)

# With @:
beta_new = np.linalg.inv(X.T @ X) @ X.T @ y

# Even cleaner with parentheses:
beta = (X.T @ X).I @ X.T @ y  # If using matrix objects

python code snippet end

Scientific Computing Applications

Machine Learning Examples

python code snippet start

import numpy as np

# Neural network forward pass
def forward_pass(X, W1, b1, W2, b2):
    """Two-layer neural network"""
    # Hidden layer
    z1 = X @ W1 + b1
    a1 = np.tanh(z1)
    
    # Output layer  
    z2 = a1 @ W2 + b2
    a2 = 1 / (1 + np.exp(-z2))  # Sigmoid
    
    return a2

# Before @ this would be:
# z1 = np.dot(X, W1) + b1
# z2 = np.dot(a1, W2) + b2

python code snippet end

Linear Algebra Operations

python code snippet start

# Gram-Schmidt orthogonalization
def gram_schmidt(vectors):
    """Orthogonalize a set of vectors"""
    orthogonal = []
    
    for v in vectors:
        # Remove projections onto previous vectors
        for u in orthogonal:
            # Projection: proj = (v @ u) / (u @ u) * u
            v = v - (v @ u) / (u @ u) * u
        
        # Normalize
        v = v / np.sqrt(v @ v)
        orthogonal.append(v)
    
    return np.array(orthogonal)

# Clear matrix operations with @
vectors = [np.array([1, 1, 0]), np.array([1, 0, 1]), np.array([0, 1, 1])]
orthogonal_vectors = gram_schmidt(vectors)

python code snippet end

Quantum Computing Simulation

python code snippet start

# Quantum gate operations
def apply_quantum_gates(state_vector, *gates):
    """Apply sequence of quantum gates to state vector"""
    result = state_vector
    
    for gate in gates:
        result = gate @ result
    
    return result

# Pauli gates
pauli_x = np.array([[0, 1], [1, 0]])
pauli_y = np.array([[0, -1j], [1j, 0]])
pauli_z = np.array([[1, 0], [0, -1]])

# Initial state |0⟩
state = np.array([1, 0])

# Apply gates: |ψ⟩ = Z @ Y @ X @ |0⟩
final_state = apply_quantum_gates(state, pauli_x, pauli_y, pauli_z)

python code snippet end

3D Graphics and Computer Vision

Transformation Matrices

python code snippet start

import numpy as np

def create_rotation_matrix(angle_x, angle_y, angle_z):
    """Create 3D rotation matrix"""
    # Rotation around X-axis
    Rx = np.array([
        [1, 0, 0],
        [0, np.cos(angle_x), -np.sin(angle_x)],
        [0, np.sin(angle_x), np.cos(angle_x)]
    ])
    
    # Rotation around Y-axis
    Ry = np.array([
        [np.cos(angle_y), 0, np.sin(angle_y)],
        [0, 1, 0],
        [-np.sin(angle_y), 0, np.cos(angle_y)]
    ])
    
    # Rotation around Z-axis
    Rz = np.array([
        [np.cos(angle_z), -np.sin(angle_z), 0],
        [np.sin(angle_z), np.cos(angle_z), 0],
        [0, 0, 1]
    ])
    
    # Combined rotation: R = Rz @ Ry @ Rx
    return Rz @ Ry @ Rx

# Transform 3D points
points = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]).T
rotation = create_rotation_matrix(0.1, 0.2, 0.3)
transformed_points = rotation @ points

python code snippet end

Camera Projection

python code snippet start

def project_3d_to_2d(points_3d, camera_matrix, world_to_camera):
    """Project 3D points to 2D image coordinates"""
    # Transform to camera coordinates
    camera_coords = world_to_camera @ points_3d
    
    # Project to image plane
    image_coords = camera_matrix @ camera_coords
    
    # Normalize by homogeneous coordinate
    image_coords = image_coords[:2] / image_coords[2]
    
    return image_coords

# Without @, this would be much less readable:
# np.dot(camera_matrix, np.dot(world_to_camera, points_3d))

python code snippet end

Library Adoption

NumPy Integration

python code snippet start

import numpy as np

# All NumPy arrays support @
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6], [7, 8]])
result = a @ b

# Works with different array types
sparse_a = np.array([[1, 0], [0, 1]])  # Identity
result = sparse_a @ b  # Still works

# Broadcasting rules apply
vector = np.array([1, 2])
matrix = np.array([[1, 2], [3, 4]])
result = matrix @ vector  # Matrix-vector multiplication

python code snippet end

Custom Class Implementation

python code snippet start

class Matrix:
    def __init__(self, data):
        self.data = np.array(data)
    
    def __matmul__(self, other):
        """Implement @ operator"""
        if isinstance(other, Matrix):
            return Matrix(self.data @ other.data)
        else:
            return Matrix(self.data @ other)
    
    def __rmatmul__(self, other):
        """Implement reverse @ operator"""
        return Matrix(other @ self.data)
    
    def __imatmul__(self, other):
        """Implement @= operator"""
        if isinstance(other, Matrix):
            self.data @= other.data
        else:
            self.data @= other
        return self

# Usage
A = Matrix([[1, 2], [3, 4]])
B = Matrix([[5, 6], [7, 8]])
C = A @ B  # Uses __matmul__

python code snippet end

Performance Benefits

Optimization Opportunities

python code snippet start

# The @ operator allows better optimization
# Chain of matrix multiplications
A, B, C, D = [np.random.rand(100, 100) for _ in range(4)]

# Before @: Manual optimization needed
temp1 = np.dot(A, B)
temp2 = np.dot(temp1, C) 
result = np.dot(temp2, D)

# With @: Automatic optimization possible
result = A @ B @ C @ D  # Optimizer can choose best order

python code snippet end

Migration from Old Code

Gradual Adoption

python code snippet start

# Old code
def old_linear_regression(X, y):
    return np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)

# New code
def new_linear_regression(X, y):
    return np.linalg.inv(X.T @ X) @ X.T @ y

# Both produce the same result!
assert np.allclose(old_linear_regression(X, y), new_linear_regression(X, y))

python code snippet end

The @ operator revolutionized scientific Python by making mathematical code readable and maintainable - it’s now the standard for matrix operations!

Reference: PEP 465 - Matrix Multiplication Operator