Multivariable Calculus

Multivariable Calculus

---

Multivariable Functions

Example: A neural network layer with 256 inputs and 64 outputs is a vector function $\vec{f}: \mathbb{R}^{256} \to \mathbb{R}^{64}$.

---

Partial Derivatives Review

Key Rules:

• Sum Rule: $\frac{\partial}{\partial x_i}(f + g) = \frac{\partial f}{\partial x_i} + \frac{\partial g}{\partial x_i}$
• Product Rule: $\frac{\partial}{\partial x_i}(f \cdot g) = \frac{\partial f}{\partial x_i} g + f \frac{\partial g}{\partial x_i}$
• Chain Rule: $\frac{\partial}{\partial x_i} f(g_1(\vec{x}), \ldots) = \sum_j \frac{\partial f}{\partial g_j} \frac{\partial g_j}{\partial x_i}$
• Clairaut's Theorem: $\frac{\partial^2 f}{\partial x_i \partial x_j} = \frac{\partial^2 f}{\partial x_j \partial x_i}$ (for continuous second derivatives)

---

The Gradient Vector

Example: For $f(x, y) = x^2 y + \sin(y)$:

$\nabla f = \left( 2xy, \; x^2 + \cos(y) \right)$

At point $(1, \pi)$: $\nabla f(1, \pi) = (2\pi, \; 1 - 1) = (2\pi, 0)$

---

Directional Derivative

Example: For $f(x,y) = x^2 + y^2$ at $(1,1)$, direction $\vec{u} = (1/\sqrt{2}, 1/\sqrt{2})$:

$\nabla f(1,1) = (2, 2)$, so $D_{\vec{u}} f = 2 \cdot \frac{1}{\sqrt{2}} + 2 \cdot \frac{1}{\sqrt{2}} = 2\sqrt{2}$

---

Hessian Matrix

Example: For $f(x, y) = x^4 + y^4 - 4xy + 1$:

$\nabla f = (4x^3 - 4y, \; 4y^3 - 4x)$ $H = \begin{bmatrix} 12x^2 & -4 \\ -4 & 12y^2 \end{bmatrix}$

At critical point $(1,1)$: $H = \begin{bmatrix} 12 & -4 \\ -4 & 12 \end{bmatrix}$

---

Multivariable Taylor Series

Second-order approximation (scalar case):

$f(x_0 + h) \approx f(x_0) + f'(x_0)h + \frac{1}{2}f''(x_0)h^2$

Vector form: $\vec{h}^T H \vec{h} = \sum_{i,j} H_{ij} h_i h_j$ is a quadratic form.

---

Critical Points Classification

Classification rules for $f(x,y)$:

Example: $f(x,y) = x^2 + y^2 - 2x - 4y + 5$

$\nabla f = (2x - 2, \; 2y - 4) = (0,0) \implies (x,y) = (1,2)$

$H = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}$, eigenvalues $\{2, 2\} > 0 \implies$ local minimum

---

Double Integrals

---

Change of Variables

Common transformations:

---

Python Implementation

import numpy as np
import sympy as sp

# ============================================
# 1. Partial Derivatives with SymPy
# ============================================
x, y = sp.symbols('x y')
f = x**2 * y + sp.sin(y)

# Partial derivatives
df_dx = sp.diff(f, x)
df_dy = sp.diff(f, y)
print(f"f = {f}")
print(f"∂f/∂x = {df_dx}")      # 2*x*y
print(f"∂f/∂y = {df_dy}")      # x**2 + cos(y)

# Gradient
grad_f = [df_dx, df_dy]
print(f"∇f = {grad_f}")

# ============================================
# 2. Hessian Matrix with SymPy
# ============================================
f_hess = x**4 + y**4 - 4*x*y + 1
H = sp.hessian(f_hess, (x, y))
print(f"\nHessian of {f_hess}:")
sp.pprint(H)

eigenvals = H.eigenvals()
print(f"Eigenvalues at (1,1): {H.subs({x:1, y:1}).eigenvals()}")

# ============================================
# 3. Gradient Descent Implementation
# ============================================
def gradient_descent(grad_f, x0, lr=0.01, max_iter=1000, tol=1e-8):
    """Simple gradient descent optimization."""
    x = np.array(x0, dtype=float)
    history = [x.copy()]
    
    for i in range(max_iter):
        g = np.array(grad_f(x), dtype=float)
        if np.linalg.norm(g) < tol:
            print(f"Converged in {i} iterations")
            break
        x = x - lr * g
        history.append(x.copy())
    
    return x, np.array(history)

# Minimize f(x,y) = x^2 + y^2 - 2x - 4y + 5
f = lambda v: v[0]**2 + v[1]**2 - 2*v[0] - 4*v[1] + 5
grad = lambda v: np.array([2*v[0]-2, 2*v[1]-4])

x_min, hist = gradient_descent(grad, [5.0, 5.0], lr=0.1)
print(f"\nMinimum at: {x_min}")  # Should be [1, 2]

# ============================================
# 4. Numerical Jacobian
# ============================================
def numerical_jacobian(f, x, h=1e-7):
    """Compute Jacobian matrix numerically using central differences."""
    n = len(x)
    m = len(f(x))
    J = np.zeros((m, n))
    for j in range(n):
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[j] += h
        x_minus[j] -= h
        J[:, j] = (f(x_plus) - f(x_minus)) / (2 * h)
    return J

# Example: f(x,y) = [x^2 + y, x*y^2]
f_vec = lambda v: np.array([v[0]**2 + v[1], v[0]*v[1]**2])
x_point = np.array([1.0, 2.0])
J = numerical_jacobian(f_vec, x_point)
print(f"\nJacobian at (1,2):\n{J}")
# Expected: [[2, 1], [4, 4]]

# ============================================
# 5. Numerical Hessian
# ============================================
def numerical_hessian(f, x, h=1e-5):
    """Compute Hessian matrix numerically."""
    n = len(x)
    H = np.zeros((n, n))
    f0 = f(x)
    
    for i in range(n):
        for j in range(n):
            x_pp = x.copy()
            x_pm = x.copy()
            x_mp = x.copy()
            x_mm = x.copy()
            
            x_pp[i] += h; x_pp[j] += h
            x_pm[i] += h; x_pm[j] -= h
            x_mp[i] -= h; x_mp[j] += h
            x_mm[i] -= h; x_mm[j] -= h
            
            H[i, j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (4 * h**2)
    
    return H

# Test on f(x,y) = x^4 + y^4 - 4xy + 1
f_test = lambda v: v[0]**4 + v[1]**4 - 4*v[0]*v[1] + 1
H_num = numerical_hessian(f_test, np.array([1.0, 1.0]))
print(f"\nNumerical Hessian at (1,1):\n{H_num}")

# ============================================
# 6. Double Integration with SciPy
# ============================================
from scipy import integrate

def integrand(y, x):
    return np.exp(-(x**2 + y**2))

# Integrate over circle of radius R
result, error = integrate.dblquad(
    integrand,
    -2, 2,           # x bounds
    lambda x: -np.sqrt(4-x**2),  # y lower
    lambda x: np.sqrt(4-x**2)    # y upper
)
print(f"\nDouble integral over circle: {result:.6f}")
print(f"Expected (≈ π(1-e^(-4))): {np.pi*(1-np.exp(-4)):.6f}")

---

Applications in AI/ML

1. Gradient Descent: The most fundamental optimization algorithm uses $\theta_{t+1} = \theta_t - \eta \nabla L(\theta_t)$. The gradient $\nabla L$ provides the direction of steepest descent.

2. Second-Order Methods: Newton's method uses the Hessian to achieve quadratic convergence:

$\theta_{t+1} = \theta_t - H^{-1}(\theta_t) \nabla L(\theta_t)$

The Hessian's eigenvalues determine convergence rate: eigenvalues close to zero slow convergence, while large condition numbers indicate ill-conditioning.

3. Natural Gradient: The Fisher information matrix (an expected Hessian) accounts for the geometry of the parameter space, providing better updates than standard gradient descent.

4. Saddle Points in High Dimensions: In high-dimensional loss landscapes, saddle points are exponentially more common than local minima. The Hessian's indefinite eigenvalues reveal saddle points that trap first-order methods.

5. Normalizing Flows: Use the change of variables formula with the Jacobian determinant to transform simple distributions into complex ones:

$\log p_X(x) = \log p_Z(f^{-1}(x)) + \log |\det(J_{f^{-1}}(x))|$

6. Backpropagation: Each layer's Jacobian matrix propagates gradients through the chain rule. Vanishing/exploding gradients correspond to Jacobian eigenvalues less than or greater than 1.

---

Common Mistakes

---

Interview Questions

---

Practice Problems

---

Quick Reference

---

Cross-References

• Linear Algebra: Eigenvalues and Eigenvectors — Essential for Hessian analysis and classification of critical points
• Optimization: Calculus Optimization — Direct application of multivariable gradients
• Probability: Probability Distributions — Joint distributions and expectations via double integrals
• Calculus I: Calculus Derivatives — Foundation for partial derivatives and gradients
• Calculus II: Calculus Integrals — Foundation for multiple integrals
• Matrix Calculus: Matrix Calculus — Efficient computation of gradients in vectorized form