mathematics-formal

Mathematics (Formal and Computational)

Overview

This skill bridges formal mathematics and scientific computing, covering symbolic computation, numerical methods, optimization theory, and common numerical pitfalls. Use it when mathematical rigor or numerical correctness is central to your research.

When to Use This Skill

• Deriving or verifying mathematical expressions symbolically
• Implementing numerically stable algorithms
• Choosing and applying optimization methods
• Working with probability distributions or information theory
• Checking gradient implementations or matrix computations

---

1. Symbolic Computation with SymPy

import sympy as sp

# Define symbolic variables
x, y, t, n = sp.symbols("x y t n", real=True)
alpha, beta = sp.symbols("alpha beta", positive=True)

# Algebra
expr = (x + y)**3
print(sp.expand(expr))       # x³ + 3x²y + 3xy² + y³
print(sp.factor(x**2 - 1))  # (x-1)(x+1)

# Calculus
f = sp.exp(-alpha * x**2)
df = sp.diff(f, x)           # derivative
integral = sp.integrate(f, (x, -sp.oo, sp.oo))  # definite integral
print(f"∫ exp(-αx²) dx = {integral}")  # √(π/α)

# Taylor series
series = sp.series(sp.sin(x), x, 0, n=7)
print(series)  # x - x³/6 + x⁵/120 - ...

# Solve equations
solutions = sp.solve(x**2 + 2*x - 3, x)  # [1, -3]

# ODEs
y_func = sp.Function("y")
ode = sp.Eq(y_func(t).diff(t) + alpha * y_func(t), 0)
sol = sp.dsolve(ode, y_func(t))
print(sol)  # y(t) = C1 * exp(-αt)

# Linear algebra
A = sp.Matrix([[1, 2], [3, 4]])
print(A.det())         # -2
print(A.eigenvals())   # {3 - √5: 1, 3 + √5: 1}
print(A.inv())

---

2. Numerical Linear Algebra

import numpy as np
from scipy import linalg

# ── Conditioning ──────────────────────────────────────────────
A = np.array([[1, 2], [2, 4.001]])  # nearly singular
cond = np.linalg.cond(A)
print(f"Condition number: {cond:.2e}")
# > 1e6 → ill-conditioned; results sensitive to input perturbations
# > 1/eps (≈ 4.5e15) → numerically singular

# ── Solving linear systems ─────────────────────────────────────
# NEVER: x = np.linalg.inv(A) @ b  (unstable, expensive)
# ALWAYS: x = np.linalg.solve(A, b)  (uses LU decomposition)
b = np.array([1.0, 2.0])
x = np.linalg.solve(A, b)

# Sparse systems (large n):
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
A_sparse = csr_matrix(A)
x_sparse = spsolve(A_sparse, b)

# ── SVD Decomposition ─────────────────────────────────────────
U, s, Vh = np.linalg.svd(A)
# Low-rank approximation (keep top k singular values)
k = 1
A_approx = (U[:, :k] * s[:k]) @ Vh[:k, :]

# Numerical rank (robust to noise)
rank = np.linalg.matrix_rank(A, tol=1e-10)

# Pseudo-inverse (for rank-deficient systems)
A_pinv = np.linalg.pinv(A)

# ── Eigendecomposition ────────────────────────────────────────
# For symmetric/Hermitian matrices (more stable):
eigenvalues, eigenvectors = np.linalg.eigh(A.T @ A)
# For general matrices:
eigenvalues_g, eigenvectors_g = np.linalg.eig(A)