Linear Algebra

from scipy import linalg
import numpy as np

# Or direct imports
from scipy.linalg import (
    solve, lstsq, inv, det, norm,
    eig, eigh, eigvals, eigvalsh,
    svd, svdvals,
    lu, qr, cholesky, ldl,
    expm, logm, sqrtm, funm,
    solve_triangular, solve_banded,
    block_diag, circulant, companion,
    kron, tril, triu
)

from scipy.linalg import solve

A = np.array([[3, 1], [1, 2]], dtype=float)
b = np.array([9, 8], dtype=float)

x = solve(A, b)           # solves A @ x = b
# x ≈ [2, 3]

# Multiple right-hand sides (b as 2-D)
B = np.random.rand(2, 5)
X = solve(A, B)           # shape (2, 5)

# Key options
x = solve(A, b,
          lower=False,        # A is upper triangular (skip factorization)
          overwrite_a=False,  # allow overwriting A for memory efficiency
          overwrite_b=False,
          check_finite=True,  # set False to skip input checking (faster)
          assume_a='gen')     # 'gen', 'sym', 'her', 'pos'

# Symmetric positive definite → faster
x = solve(A, b, assume_a='pos')

from scipy.linalg import solve_triangular

R = np.triu(np.random.rand(4, 4)) + np.eye(4)  # upper triangular
b = np.random.rand(4)

x = solve_triangular(R, b)                   # default: upper
x = solve_triangular(R, b, lower=True)       # treat as lower
x = solve_triangular(R, b, trans='T')        # solve R.T @ x = b
x = solve_triangular(R, b, unit_diagonal=True)

from scipy.linalg import solve_banded

# ab[0] = upper diagonal, ab[1] = main, ab[2] = lower (for l=1, u=1)
ab = np.array([[0, 1, 1, 1],    # upper diag (first element unused)
               [4, 4, 4, 4],    # main diagonal
               [1, 1, 1, 0]])   # lower diag (last element unused)
b = np.array([1, 2, 3, 4], dtype=float)

x = solve_banded((1, 1), ab, b)  # (l, u) = number of sub/superdiagonals

from scipy.linalg import lstsq

A = np.random.rand(10, 3)
b = np.random.rand(10)

x, residuals, rank, sv = lstsq(A, b)
# x         — solution minimizing ||A@x - b||
# residuals — sum of squared residuals (empty if rank < N)
# rank      — effective rank of A
# sv        — singular values

from scipy.linalg import lu, lu_factor, lu_solve

P, L, U = lu(A)             # A = P @ L @ U

# Efficient repeated solves with same A
lu_and_piv = lu_factor(A)   # factorize once
x = lu_solve(lu_and_piv, b) # solve many times
x2 = lu_solve(lu_and_piv, b2)

from scipy.linalg import qr, qr_update, qr_delete, qr_insert

Q, R = qr(A)                     # full QR: Q (m×m), R (m×n)
Q, R = qr(A, mode='economic')    # thin QR: Q (m×n), R (n×n)
R_only = qr(A, mode='r')         # R only (faster)
Q, R, P = qr(A, pivoting=True)   # column-pivoted: A[:,P] = Q @ R

from scipy.linalg import cholesky, cho_factor, cho_solve

# A must be symmetric positive definite
L = cholesky(A, lower=True)    # A = L @ L.T (lower triangular)
U = cholesky(A, lower=False)   # A = U.T @ U (upper triangular)

# Efficient solves
c, low = cho_factor(A, lower=True)
x = cho_solve((c, low), b)

from scipy.linalg import ldl

# For symmetric (not necessarily positive definite)
L, D, perm = ldl(A)   # A[np.ix_(perm,perm)] = L @ D @ L.T

from scipy.linalg import svd, svdvals

U, s, Vh = svd(A)              # full SVD: A = U @ diag(s) @ Vh
U, s, Vh = svd(A, full_matrices=False)  # thin/economy SVD

s_only = svdvals(A)            # just singular values (faster)

# Reconstruct
A_reconstructed = U @ np.diag(s) @ Vh

# Truncated SVD (use scipy.sparse.linalg for large matrices)
from scipy.sparse.linalg import svds
U_k, s_k, Vh_k = svds(A, k=5)  # top-5 singular values/vectors

from scipy.linalg import schur, rsf2csf

T, Z = schur(A)               # A = Z @ T @ Z.H (T quasi-upper triangular)
T, Z = schur(A, output='complex')  # force complex Schur form

from scipy.linalg import eig, eigh, eigvals, eigvalsh

# General (non-symmetric)
eigenvalues, eigenvectors = eig(A)
eigenvalues = eigvals(A)          # eigenvalues only (faster)

# Symmetric/Hermitian (real eigenvalues, orthonormal eigenvectors)
eigenvalues, eigenvectors = eigh(A)    # sorted ascending
eigenvalues = eigvalsh(A)

# Generalized eigenvalue problem: A @ v = lambda * B @ v
eigenvalues, eigenvectors = eig(A, b=B)
eigenvalues, eigenvectors = eigh(A, b=B)  # symmetric A, B

# Subset (faster for large matrices)
eigenvalues, eigenvectors = eigh(A, subset_by_index=[0, 4])  # first 5
eigenvalues, eigenvectors = eigh(A, subset_by_value=[-1, 1]) # in [-1,1]

# Left and right eigenvectors
eigenvalues, vl, vr = eig(A, left=True, right=True)

from scipy.linalg import inv, det, norm, matrix_rank, pinv

inv_A = inv(A)                         # inverse (avoid if possible — use solve)
d = det(A)                             # determinant
r = matrix_rank(A)                     # numerical rank
A_pinv = pinv(A)                       # Moore-Penrose pseudoinverse

# Norms
n1 = norm(A, ord=1)     # max column sum
n2 = norm(A, ord=2)     # spectral norm (= largest singular value)
nf = norm(A, 'fro')     # Frobenius
ni = norm(A, np.inf)    # max row sum

# Condition number
cond = norm(A) * norm(inv(A))
# or faster:
cond = np.linalg.cond(A)

from scipy.linalg import expm, logm, sqrtm, signm, funm, cosm, sinm

eA = expm(A)          # matrix exponential e^A
lA = logm(A)          # matrix logarithm
sA = sqrtm(A)         # matrix square root
sgA = signm(A)        # matrix sign

# General matrix function
fA = funm(A, np.cos)  # cos(A) via Schur decomposition

# Trigonometric
cA = cosm(A)
sA = sinm(A)

from scipy.linalg import block_diag, circulant, companion, toeplitz, hankel, hadamard, hilbert, leslie

# Block diagonal
BD = block_diag(A1, A2, A3)

# Circulant (from first column)
C = circulant([1, 2, 3, 4])

# Companion matrix of polynomial
# For x^3 - 6x^2 + 11x - 6 = 0:
M = companion([1, -6, 11, -6])

# Toeplitz
T = toeplitz([1, 2, 3], [1, 4, 5])   # (col, row) vectors

# Hilbert (ill-conditioned test matrix)
H = hilbert(5)

# Hadamard
W = hadamard(8)   # n must be power of 2

# Kronecker product
from scipy.linalg import kron
K = kron(A, B)    # Kronecker product

# Solve with preconditioning (large systems — see also scipy.sparse.linalg)
from scipy.sparse.linalg import spsolve, cg, gmres
from scipy.sparse import csr_array

# For large, sparse A:
A_sparse = csr_array(A)
x = spsolve(A_sparse, b)

# Iterative solvers
x, info = cg(A_sparse, b, tol=1e-8, maxiter=1000)    # Conjugate Gradient (SPD)
x, info = gmres(A_sparse, b, tol=1e-8)               # GMRES (general)