Linear Algebra

import numpy as np

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

# Matrix multiplication (dot product for 2-D)
A @ B                    # preferred syntax (PEP 465)
np.matmul(A, B)          # equivalent — no scalar broadcasting
np.dot(A, B)             # same for 2-D; for N-D behaves differently

# Dot product for 1-D arrays (inner product)
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.dot(a, b)             # 32  (scalar)
a @ b                    # 32

# Outer product
np.outer(a, b)           # (3, 3) matrix
np.outer(a, b) == a[:, None] * b[None, :]  # True

# Inner product (generalized)
np.inner(A, B)           # sum over last axis of A and last axis of B

# Batch matrix multiply (3-D: batch × M × K  @  batch × K × N)
batch_A = np.random.rand(10, 3, 4)
batch_B = np.random.rand(10, 4, 5)
batch_A @ batch_B        # (10, 3, 5)

# einsum — expressive contraction notation
np.einsum("ij,jk->ik", A, B)        # matrix multiply
np.einsum("i,i->", a, b)            # dot product
np.einsum("ij->ji", A)              # transpose
np.einsum("ij,ij->ij", A, B)        # element-wise multiply
np.einsum("ijk,ikl->ijl", batch_A, batch_B)  # batch matmul
np.einsum("ii->i", A)               # diagonal
np.einsum("ii", A)                  # trace
np.einsum("ij->i", A)               # row sums

from numpy import linalg as la
# or: import numpy.linalg as la

vals, vecs = la.eig(A)      # eigenvalues, eigenvectors (cols of vecs)
vals                         # array of eigenvalues (may be complex)
vecs[:, 0]                   # first eigenvector

vals = la.eigvals(A)         # eigenvalues only (faster)

# Symmetric / Hermitian matrices (real eigenvalues guaranteed)
vals, vecs = la.eigh(A)     # eigenvalues sorted ascending
vals = la.eigvalsh(A)        # eigenvalues only

U, S, Vh = la.svd(A)               # A = U @ np.diag(S) @ Vh
U, S, Vh = la.svd(A, full_matrices=False)  # economy / thin SVD
S = la.svdvals(A)                    # singular values only (NumPy ≥ 2.0)

# Reconstruct: A ≈ U @ np.diag(S) @ Vh  (full) or  (U * S) @ Vh (thin)
k = 2
A_approx = (U[:, :k] * S[:k]) @ Vh[:k, :]   # rank-k approximation

Q, R = la.qr(A)              # A = Q @ R
Q, R = la.qr(A, mode="reduced")   # economy form (default for non-square)
Q, R = la.qr(A, mode="complete")  # full Q
Q, R = la.qr(A, mode="r")         # R only

# A must be symmetric positive-definite
L = la.cholesky(A)           # A = L @ L.T.conj()
# NumPy returns lower-triangular L

# Solve Ax = b
A = np.array([[2.0, 1.0], [1.0, 3.0]])
b = np.array([5.0, 10.0])

x = la.solve(A, b)                   # exact solution (A must be square)
la.solve(A, np.stack([b, b], axis=1))  # multiple right-hand sides

x, res, rank, sv = la.lstsq(A, b, rcond=None)  # least-squares (A can be non-square)
# res = residuals (sum of squared residuals), rank = effective rank, sv = singular values

# Triangular solve
np.linalg.solve(np.tril(A), b)   # use la.solve; no dedicated triangular variant in NumPy

la.inv(A)          # matrix inverse (A must be square and non-singular)
la.pinv(A)         # Moore-Penrose pseudo-inverse (works for any shape)
la.det(A)          # determinant (scalar)
la.slogdet(A)      # (sign, log|det|) — numerically stable for large matrices
la.matrix_rank(A)  # effective rank (using SVD + tolerance)
la.matrix_rank(A, tol=1e-10)   # custom tolerance

la.norm(a)                     # L2 norm of vector
la.norm(a, ord=1)              # L1 norm
la.norm(a, ord=np.inf)        # max absolute value
la.norm(a, ord=-np.inf)       # min absolute value

la.norm(A)                     # Frobenius norm (default for matrices)
la.norm(A, ord="fro")          # same
la.norm(A, ord=2)              # spectral norm (largest singular value)
la.norm(A, ord=1)              # max column sum
la.norm(A, ord=np.inf)        # max row sum
la.norm(A, ord=-1)             # min column sum
la.norm(A, ord="nuc")          # nuclear norm (sum of singular values)

la.norm(A, axis=0)             # per-column norm → shape (n,)
la.norm(A, axis=1, keepdims=True)   # per-row norm → shape (m, 1)

la.cond(A)          # condition number (ratio of largest to smallest singular value)
la.cond(A, p=1)     # condition number in 1-norm

np.trace(A)               # sum of diagonal elements
np.trace(A, offset=1)     # sum of k-th diagonal
np.diagonal(A)            # extract main diagonal as view
np.diagonal(A, offset=1)  # k-th diagonal
np.fill_diagonal(A, 5)    # fill main diagonal in-place (no return value)

# Triangular extraction
np.triu(A)           # upper triangle (zeros below diagonal)
np.tril(A)           # lower triangle (zeros above diagonal)
np.triu(A, k=1)      # above main diagonal only
np.tril(A, k=-1)     # below main diagonal only

# Kronecker product
np.kron(A, B)        # block matrix

# Cross product (3-D vectors)
np.cross(a, b)       # a × b
np.cross(a, b, axisa=0, axisb=0)  # specify which axis holds the vectors

Q, _ = la.qr(A)               # orthonormal basis
proj = Q @ (Q.T @ b)          # project vector b onto column space of A

A = np.random.rand(4, 4)
B = np.random.rand(4, 100)    # 100 right-hand sides
X = la.solve(A, B)            # X.shape == (4, 100)

def is_pos_def(M):
    try:
        la.cholesky(M)
        return True
    except la.LinAlgError:
        return False

def null_space(A, tol=1e-10):
    _, S, Vh = la.svd(A, full_matrices=True)
    return Vh[np.sum(S > tol):].T

try:
    la.inv(singular_matrix)
except la.LinAlgError as e:
    print(e)   # "Singular matrix"

try:
    la.cholesky(non_posdef)
except la.LinAlgError as e:
    print(e)   # "Matrix is not positive definite"