Linear Algebra

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: ndarray, optional, output array
    
    Returns:
    - ndarray: Dot product of a and b
    """

def matmul(x1, x2, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True):
    """
    Matrix product of two arrays.
    
    Parameters:
    - x1: array_like, first input array
    - x2: array_like, second input array
    - out: ndarray, optional, output array
    - casting: casting rule
    - order: memory layout
    - dtype: data type, output type
    - subok: bool, whether to return subclass
    
    Returns:
    - ndarray: Matrix product of inputs
    """

def vdot(a, b):
    """
    Return dot product of two vectors (flattened arrays).
    
    Parameters:
    - a: array_like, first input array
    - b: array_like, second input array
    
    Returns:
    - scalar: Dot product of flattened arrays
    """

def inner(a, b):
    """
    Inner product of two arrays.
    
    Parameters:
    - a: array_like, first input array
    - b: array_like, second input array
    
    Returns:
    - ndarray: Inner product
    """

def outer(a, b, out=None):
    """
    Compute outer product of two vectors.
    
    Parameters:
    - a: array_like, first input vector
    - b: array_like, second input vector
    - out: ndarray, optional, output array
    
    Returns:
    - ndarray: Outer product matrix
    """

def tensordot(a, b, axes=2):
    """
    Compute tensor dot product along specified axes.
    
    Parameters:
    - a: array_like, first input array
    - b: array_like, second input array
    - axes: int or (2,) array_like, axes to sum over
    
    Returns:
    - ndarray: Tensor dot product
    """

def einsum(subscripts, *operands, **kwargs):
    """
    Evaluates Einstein summation convention on operands.
    
    Parameters:
    - subscripts: str, subscripts for summation
    - *operands: list of array_like, arrays for operation
    - optimize: {False, True, 'greedy', 'optimal'}, optimization strategy
    - out: ndarray, optional, output array
    
    Returns:
    - ndarray: Calculation based on Einstein summation
    """

def kron(a, b):
    """
    Kronecker product of two arrays.
    
    Parameters:
    - a: array_like, first input array
    - b: array_like, second input array
    
    Returns:
    - ndarray: Kronecker product
    """

def cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None):
    """
    Return cross product of two (arrays of) vectors.
    
    Parameters:
    - a: array_like, first input array
    - b: array_like, second input array
    - axisa: int, axis of a defining vectors
    - axisb: int, axis of b defining vectors  
    - axisc: int, axis of c containing cross product
    - axis: int, if defined, axis of inputs defining vectors
    
    Returns:
    - ndarray: Cross product of a and b
    """

def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
    """
    Singular Value Decomposition.
    
    Parameters:
    - a: array_like, input matrix (..., M, N)
    - full_matrices: bool, compute full or reduced SVD
    - compute_uv: bool, compute U and Vh matrices
    - hermitian: bool, assume input is Hermitian
    
    Returns:
    - U: ndarray, unitary array (..., M, M) or (..., M, K)
    - s: ndarray, singular values (..., K)
    - Vh: ndarray, unitary array (..., N, N) or (..., K, N)
    where K = min(M, N)
    """

def qr(a, mode='reduced'):
    """
    Compute QR decomposition of matrix.
    
    Parameters:
    - a: array_like, input matrix (..., M, N)
    - mode: {'reduced', 'complete', 'r', 'raw'}, decomposition mode
    
    Returns:
    - Q: ndarray, orthogonal matrix
    - R: ndarray, upper triangular matrix
    """

def cholesky(a):
    """
    Cholesky decomposition of positive-definite matrix.
    
    Parameters:
    - a: array_like, positive-definite matrix (..., N, N)
    
    Returns:
    - L: ndarray, lower triangular matrix such that a = L @ L.H
    """

def eigh(a, UPLO='L'):
    """
    Return eigenvalues and eigenvectors of Hermitian matrix.
    
    Parameters:
    - a: array_like, Hermitian matrix (..., N, N)
    - UPLO: {'L', 'U'}, whether to use lower or upper triangle
    
    Returns:
    - w: ndarray, eigenvalues in ascending order (..., N)
    - v: ndarray, normalized eigenvectors (..., N, N)
    """

def eigvalsh(a, UPLO='L'):
    """
    Return eigenvalues of Hermitian matrix.
    
    Parameters:
    - a: array_like, Hermitian matrix (..., N, N)
    - UPLO: {'L', 'U'}, whether to use lower or upper triangle
    
    Returns:
    - w: ndarray, eigenvalues in ascending order (..., N)
    """

def norm(x, ord=None, axis=None, keepdims=False):
    """
    Matrix or vector norm.
    
    Parameters:
    - x: array_like, input array
    - ord: {non-zero int, inf, -inf, 'fro', 'nuc'}, order of norm
    - axis: {None, int, 2-tuple of ints}, axis for norm calculation
    - keepdims: bool, keep reduced dimensions
    
    Returns:
    - ndarray: Norm of the matrix or vector
    """

def det(a):
    """
    Compute determinant of array.
    
    Parameters:
    - a: array_like, square matrix (..., N, N)
    
    Returns:
    - ndarray: Determinant of a
    """

def slogdet(a):
    """
    Compute sign and natural logarithm of determinant.
    
    Parameters:
    - a: array_like, square matrix (..., N, N)
    
    Returns:
    - sign: ndarray, sign of determinant
    - logabsdet: ndarray, natural log of absolute determinant
    """

def matrix_rank(M, tol=None, hermitian=False):
    """
    Return matrix rank using SVD method.
    
    Parameters:
    - M: array_like, input matrix
    - tol: float, optional, threshold for small singular values
    - hermitian: bool, assume M is Hermitian
    
    Returns:
    - rank: int, rank of matrix
    """

def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
    """
    Return sum along diagonals of array.
    
    Parameters:
    - a: array_like, input array
    - offset: int, offset from main diagonal
    - axis1: int, first axis for 2D sub-arrays
    - axis2: int, second axis for 2D sub-arrays
    - dtype: data type, output type
    - out: ndarray, optional, output array
    
    Returns:
    - ndarray: Sum along diagonal
    """

def solve(a, b):
    """
    Solve linear system ax = b.
    
    Parameters:
    - a: array_like, coefficient matrix (..., N, N)
    - b: array_like, ordinate values (..., N) or (..., N, K)
    
    Returns:
    - x: ndarray, solution to system (..., N) or (..., N, K)
    """

def lstsq(a, b, rcond=None):
    """
    Return least-squares solution to linear matrix equation.
    
    Parameters:
    - a: array_like, coefficient matrix (M, N)
    - b: array_like, ordinate values (M,) or (M, K)
    - rcond: float, optional, cutoff for small singular values
    
    Returns:
    - x: ndarray, least-squares solution (N,) or (N, K)
    - residuals: ndarray, residuals sum of squares
    - rank: int, rank of matrix a
    - s: ndarray, singular values of a
    """

def inv(a):
    """
    Compute multiplicative inverse of matrix.
    
    Parameters:
    - a: array_like, square matrix (..., N, N)
    
    Returns:
    - ainv: ndarray, multiplicative inverse (..., N, N)
    """

def pinv(a, rcond=1e-15, hermitian=False):
    """
    Compute Moore-Penrose pseudoinverse of matrix.
    
    Parameters:
    - a: array_like, matrix to be pseudo-inverted (..., M, N)
    - rcond: float, cutoff for small singular values
    - hermitian: bool, assume a is Hermitian
    
    Returns:
    - B: ndarray, pseudoinverse of a (..., N, M)  
    """

def tensorsolve(a, b, axes=None):
    """
    Solve tensor equation a x = b for x.
    
    Parameters:
    - a: array_like, coefficient tensor
    - b: array_like, right-hand side tensor
    - axes: list of ints, axes in a to reorder to right
    
    Returns:
    - x: ndarray, solution tensor
    """

def tensorinv(a, ind=2):
    """
    Compute 'inverse' of N-dimensional array.
    
    Parameters:
    - a: array_like, tensor to invert
    - ind: int, number of first indices to be involved in inverse
    
    Returns:
    - b: ndarray, inverse tensor
    """

def matrix_power(a, n):
    """
    Raise square matrix to integer power.
    
    Parameters:
    - a: array_like, square matrix (..., N, N)
    - n: int, exponent (can be negative)
    
    Returns:
    - a**n: ndarray, matrix raised to power n
    """

import cupy as cp

# Create matrices
A = cp.array([[1, 2], [3, 4]])
B = cp.array([[5, 6], [7, 8]])
v = cp.array([1, 2])

# Matrix multiplication
C = cp.dot(A, B)  # or A @ B
matrix_vector = cp.dot(A, v)  # or A @ v

# Matrix-matrix product with matmul
result = cp.matmul(A, B)

# Einstein summation
# Matrix multiplication: C_ij = A_ik * B_kj
C_einsum = cp.einsum('ik,kj->ij', A, B)

# Batch matrix multiplication
batch_A = cp.random.random((10, 3, 3))
batch_B = cp.random.random((10, 3, 3))
batch_result = cp.matmul(batch_A, batch_B)

import cupy as cp

# Create test matrix
A = cp.random.random((5, 4))
square_A = cp.dot(A.T, A)  # Make positive definite

# SVD decomposition
U, s, Vh = cp.linalg.svd(A, full_matrices=False)
print("SVD shapes:", U.shape, s.shape, Vh.shape)

# QR decomposition
Q, R = cp.linalg.qr(A)
print("QR reconstruction error:", cp.linalg.norm(A - cp.dot(Q, R)))

# Cholesky decomposition (for positive definite matrices)
L = cp.linalg.cholesky(square_A)
print("Cholesky reconstruction error:", cp.linalg.norm(square_A - cp.dot(L, L.T)))

import cupy as cp

# Create symmetric matrix
n = 100
A = cp.random.random((n, n))
symmetric_A = (A + A.T) / 2

# Eigenvalue decomposition
eigenvals, eigenvecs = cp.linalg.eigh(symmetric_A)
print("Eigenvalues range:", eigenvals.min(), "to", eigenvals.max())

# Verify eigenvector equation: A @ v = λ * v
idx = -1  # Largest eigenvalue
lambda_max = eigenvals[idx]
v_max = eigenvecs[:, idx]
Av = cp.dot(symmetric_A, v_max)
lambda_v = lambda_max * v_max
error = cp.linalg.norm(Av - lambda_v)
print("Eigenvector equation error:", error)

import cupy as cp

# Set up linear system Ax = b
n = 1000
A = cp.random.random((n, n)) + cp.eye(n) * 5  # Well-conditioned
x_true = cp.random.random(n)
b = cp.dot(A, x_true)

# Solve linear system
x_solved = cp.linalg.solve(A, b)
solve_error = cp.linalg.norm(x_solved - x_true)
print("Solution error:", solve_error)

# Matrix inversion
A_inv = cp.linalg.inv(A)
identity_error = cp.linalg.norm(cp.dot(A, A_inv) - cp.eye(n))
print("Inversion error:", identity_error)

# Least squares for overdetermined system
m, n = 1000, 500
A_over = cp.random.random((m, n))
x_true = cp.random.random(n)
b_over = cp.dot(A_over, x_true) + 0.01 * cp.random.random(m)

x_lstsq, residuals, rank, s = cp.linalg.lstsq(A_over, b_over, rcond=None)
print("Least squares residual:", residuals[0] if len(residuals) > 0 else 0)

import cupy as cp

# Batch operations
batch_size = 50
matrix_size = 10
batch_matrices = cp.random.random((batch_size, matrix_size, matrix_size))

# Batch determinants
batch_dets = cp.linalg.det(batch_matrices)
print("Batch determinants shape:", batch_dets.shape)

# Batch matrix norms
batch_norms = cp.linalg.norm(batch_matrices, axis=(1, 2))
print("Batch norms shape:", batch_norms.shape)

# Complex operations with broadcasting
A = cp.random.random((100, 3, 3))
B = cp.random.random((100, 3, 3))
C = cp.matmul(A, B)  # Batch matrix multiplication

# Tensor operations
tensor_A = cp.random.random((2, 3, 4, 5))
tensor_B = cp.random.random((2, 3, 5, 6))
tensor_result = cp.matmul(tensor_A, tensor_B)  # Shape: (2, 3, 4, 6)

class LinAlgError(Exception):
    """
    Generic linear algebra error.
    
    Raised when linear algebra operations encounter errors such as:
    - Singular matrices in inversion
    - Non-positive-definite matrices in Cholesky decomposition
    - Convergence failures in iterative algorithms
    """