NumPy & SciPy-compatible array library for GPU-accelerated computing with Python
—
GPU-accelerated linear algebra operations including matrix multiplication, decompositions, eigenvalue computation, and solving linear systems. CuPy provides both basic linear algebra functions in the main namespace and advanced operations in the cupy.linalg module.
Core matrix multiplication and related operations.
def dot(a, b, out=None):
"""
Dot product of two arrays.
Parameters:
- a: array-like, first input array
- b: array-like, second input array
- out: cupy.ndarray, output array
Returns:
cupy.ndarray: Dot product
- If both a and b are 1-D: inner product
- If both a and b are 2-D: matrix multiplication
- If a is N-D and b is 1-D: sum product over last axis of a and b
- If a is N-D and b is M-D (M>=2): sum product over last axis of a and second-to-last axis of b
"""
def matmul(x1, x2, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True):
"""
Matrix product of two arrays.
Parameters:
- x1, x2: array-like, input arrays
- out: cupy.ndarray, output array
- casting: {'no', 'equiv', 'safe', 'same_kind', 'unsafe'}, casting rule
- order: {'K', 'A', 'C', 'F'}, memory layout
- dtype: data type, type of output
- subok: bool, whether to allow subclasses
Returns:
cupy.ndarray: Matrix product
- No scalars allowed (unlike dot)
- Handles batch matrix multiplication for N-D arrays
"""
def inner(a, b):
"""
Inner product of two arrays.
Parameters:
- a, b: array-like, input arrays
Returns:
cupy.ndarray: Inner product over last axes
"""
def outer(a, b, out=None):
"""
Outer product of two vectors.
Parameters:
- a: array-like, first vector (flattened)
- b: array-like, second vector (flattened)
- out: cupy.ndarray, output array
Returns:
cupy.ndarray: Outer product matrix
"""
def vdot(a, b):
"""
Dot product of two vectors with complex conjugation.
Parameters:
- a, b: array-like, input arrays (flattened)
Returns:
scalar: sum(conjugate(a) * b)
"""
def tensordot(a, b, axes=2):
"""
Compute tensor dot product along specified axes.
Parameters:
- a, b: array-like, input arrays
- axes: int or (2,) array-like, axes to sum over
Returns:
cupy.ndarray: Tensor dot product
"""
def kron(a, b):
"""
Kronecker product of two arrays.
Parameters:
- a, b: array-like, input arrays
Returns:
cupy.ndarray: Kronecker product
"""
def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
"""
Cross product of two arrays.
Parameters:
- a, b: array-like, input arrays
- axisa, axisb: int, axis of a and b to take cross product over
- axisc: int, axis of output to store cross product
- axis: int, axis to take cross product over (alternative to axisa/axisb)
Returns:
cupy.ndarray: Cross product
"""Efficient tensor operations using Einstein summation convention.
def einsum(subscripts, *operands, out=None, dtype=None, order='K', casting='safe', optimize=False):
"""
Einstein summation convention on operands.
Parameters:
- subscripts: str, subscripts for summation
- operands: arrays, input arrays
- out: cupy.ndarray, output array
- dtype: data type, type of output
- order: {'K', 'A', 'C', 'F'}, memory layout
- casting: {'no', 'equiv', 'safe', 'same_kind', 'unsafe'}, casting rule
- optimize: bool or str, whether to optimize contraction order
Returns:
cupy.ndarray: Result of Einstein summation
Examples:
- 'ij,jk->ik': matrix multiplication
- 'ii->i': diagonal extraction
- 'ii': trace
- 'ij->ji': transpose
"""Advanced matrix decompositions from cupy.linalg.
def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
"""
Singular Value Decomposition.
Parameters:
- a: array-like, input matrix
- full_matrices: bool, compute full or reduced SVD
- compute_uv: bool, compute U and Vh or just singular values
- hermitian: bool, assume input is Hermitian
Returns:
tuple: (U, s, Vh) where a = U @ diag(s) @ Vh
- U: left singular vectors
- s: singular values
- Vh: right singular vectors (conjugate transposed)
"""
def qr(a, mode='reduced'):
"""
QR decomposition.
Parameters:
- a: array-like, input matrix
- mode: {'reduced', 'complete', 'r', 'raw'}, decomposition mode
Returns:
tuple: (Q, R) where a = Q @ R
- Q: orthogonal matrix
- R: upper triangular matrix
"""
def cholesky(a):
"""
Cholesky decomposition of positive definite matrix.
Parameters:
- a: array-like, Hermitian positive-definite matrix
Returns:
cupy.ndarray: Lower triangular Cholesky factor L where a = L @ L.H
"""Eigenvalue decompositions for Hermitian matrices.
def eigh(a, UPLO='L'):
"""
Eigenvalues and eigenvectors of Hermitian matrix.
Parameters:
- a: array-like, Hermitian matrix
- UPLO: {'L', 'U'}, use lower or upper triangle
Returns:
tuple: (eigenvalues, eigenvectors)
- eigenvalues: 1-D array of eigenvalues in ascending order
- eigenvectors: 2-D array with eigenvectors as columns
"""
def eigvalsh(a, UPLO='L'):
"""
Eigenvalues of Hermitian matrix.
Parameters:
- a: array-like, Hermitian matrix
- UPLO: {'L', 'U'}, use lower or upper triangle
Returns:
cupy.ndarray: 1-D array of eigenvalues in ascending order
"""Functions to compute matrix properties and norms.
def norm(x, ord=None, axis=None, keepdims=False):
"""
Matrix or vector norm.
Parameters:
- x: array-like, input array
- ord: norm type (None, 'fro', 'nuc', inf, -inf, int, float)
- axis: int or 2-tuple of ints, axes to compute norm over
- keepdims: bool, keep reduced dimensions as size 1
Returns:
cupy.ndarray: Norm of x
Matrix norms (2-D arrays):
- ord=None/'fro': Frobenius norm
- ord='nuc': nuclear norm (sum of singular values)
- ord=1: max(sum(abs(x), axis=0))
- ord=-1: min(sum(abs(x), axis=0))
- ord=2: largest singular value
- ord=-2: smallest singular value
- ord=inf: max(sum(abs(x), axis=1))
- ord=-inf: min(sum(abs(x), axis=1))
Vector norms (1-D arrays):
- ord=None/2: 2-norm (Euclidean)
- ord=inf: max(abs(x))
- ord=-inf: min(abs(x))
- ord=0: number of non-zero elements
- ord>0: sum(abs(x)**ord)**(1/ord)
"""
def det(a):
"""
Determinant of array.
Parameters:
- a: array-like, square matrix
Returns:
cupy.ndarray: Determinant of a
"""
def slogdet(a):
"""
Sign and logarithm of determinant.
Parameters:
- a: array-like, square matrix
Returns:
tuple: (sign, logdet) where det = sign * exp(logdet)
- More numerically stable than det for large matrices
"""
def matrix_rank(M, tol=None, hermitian=False):
"""
Matrix rank using SVD.
Parameters:
- M: array-like, input matrix
- tol: float, threshold for singular values
- hermitian: bool, assume input is Hermitian
Returns:
int: Rank of matrix
"""
def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
"""
Sum along diagonals of array.
Parameters:
- a: array-like, input array
- offset: int, diagonal offset
- axis1, axis2: int, axes to trace over
- dtype: data type, type of output
- out: cupy.ndarray, output array
Returns:
cupy.ndarray: Sum along diagonal
"""Functions for solving linear equations and matrix inversion.
def solve(a, b):
"""
Solve linear system ax = b.
Parameters:
- a: array-like, coefficient matrix
- b: array-like, dependent variable values
Returns:
cupy.ndarray: Solution x to ax = b
"""
def lstsq(a, b, rcond=None):
"""
Least-squares solution to linear system.
Parameters:
- a: array-like, coefficient matrix
- b: array-like, dependent variable values
- rcond: float, cutoff for small singular values
Returns:
tuple: (x, residuals, rank, s)
- x: least-squares solution
- residuals: sum of residuals
- rank: rank of matrix a
- s: singular values of a
"""
def inv(a):
"""
Matrix inverse.
Parameters:
- a: array-like, square matrix to invert
Returns:
cupy.ndarray: Inverse of a
"""
def pinv(a, rcond=1e-15, hermitian=False):
"""
Moore-Penrose pseudoinverse.
Parameters:
- a: array-like, matrix to pseudoinvert
- rcond: float, cutoff for small singular values
- hermitian: bool, assume input is Hermitian
Returns:
cupy.ndarray: Pseudoinverse of a
"""
def tensorsolve(a, b, axes=None):
"""
Solve tensor equation ax = b for x.
Parameters:
- a: array-like, coefficient tensor
- b: array-like, dependent variable tensor
- axes: list of ints, axes in a to reorder to front
Returns:
cupy.ndarray: Solution x
"""
def tensorinv(a, ind=2):
"""
Inverse of tensor using SVD.
Parameters:
- a: array-like, tensor to invert
- ind: int, number of first indices forming square matrix
Returns:
cupy.ndarray: Inverse of tensor
"""Specialized matrix operations.
def matrix_power(a, n):
"""
Raise square matrix to integer power.
Parameters:
- a: array-like, square matrix
- n: int, exponent
Returns:
cupy.ndarray: a**n
"""import cupy as cp
# Matrix multiplication
A = cp.random.random((1000, 500))
B = cp.random.random((500, 800))
# Different multiplication methods
C1 = cp.dot(A, B) # Basic dot product
C2 = cp.matmul(A, B) # Matrix multiplication
C3 = A @ B # Operator syntax
# Batch matrix multiplication
batch_A = cp.random.random((10, 100, 50))
batch_B = cp.random.random((10, 50, 75))
batch_C = batch_A @ batch_B # Shape: (10, 100, 75)# Solve system Ax = b
A = cp.random.random((100, 100))
b = cp.random.random((100, 10))
# Direct solution for square systems
x = cp.linalg.solve(A, b)
# Least squares for overdetermined systems
A_over = cp.random.random((200, 100))
b_over = cp.random.random((200,))
x_ls, residuals, rank, s = cp.linalg.lstsq(A_over, b_over)
# Matrix inversion (avoid when possible, use solve instead)
A_inv = cp.linalg.inv(A)# Singular Value Decomposition
matrix = cp.random.random((500, 300))
U, s, Vh = cp.linalg.svd(matrix)
# QR decomposition for least squares
A = cp.random.random((1000, 100))
Q, R = cp.linalg.qr(A)
# Eigendecomposition for symmetric matrices
symmetric = cp.random.random((200, 200))
symmetric = (symmetric + symmetric.T) / 2 # Make symmetric
eigenvals, eigenvecs = cp.linalg.eigh(symmetric)# Einstein summation for complex tensor operations
A = cp.random.random((10, 20, 30))
B = cp.random.random((30, 40))
# Contract last axis of A with first axis of B
result = cp.einsum('ijk,kl->ijl', A, B)
# Batch trace operation
batch_matrices = cp.random.random((50, 100, 100))
traces = cp.einsum('ijj->i', batch_matrices) # Trace of each matrix
# Matrix norms and properties
frobenius_norm = cp.linalg.norm(matrix, 'fro')
spectral_norm = cp.linalg.norm(matrix, 2)
matrix_rank = cp.linalg.matrix_rank(matrix)
determinant = cp.linalg.det(matrix[:300, :300]) # Square submatrix# Use out parameter to avoid allocations
result = cp.empty((1000, 800))
cp.dot(A, B, out=result)
# In-place operations when possible
matrix += cp.eye(matrix.shape[0]) # Add identity to diagonal
# Reuse decompositions
U, s, Vh = cp.linalg.svd(data_matrix)
# Use U, s, Vh for multiple operations without recomputingInstall with Tessl CLI
npx tessl i tessl/pypi-cupy