6a:Tc01,import numpy as np from scipy.linalg import cho_solve, cho_factor, solve from typing import Callable, Tuple, Optional def newtons_method( f: Callable[[np.ndarray], float], grad_f: Callable[[np.ndarray], np.ndarray], hess_f: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, tol: float = 1e-8, max_iter: int = 100, verbose: bool = False ) -> Tuple[np.ndarray, list]: """ Pure Newton's method for unconstrained optimization. Parameters: ----------- f : callable Objective function f(x) -> scalar grad_f : callable Gradient function ∇f(x) -> n-vector hess_f : callable Hessian function ∇²f(x) -> n×n matrix x0 : np.ndarray Starting point (n,) tol : float Gradient norm tolerance for convergence max_iter : int Maximum number of iterations verbose : bool Print iteration information Returns: -------- x : np.ndarray Approximate minimizer history : list List of (x, f(x), ||grad||) at each iteration """ x = x0.copy().astype(float) history = [] for k in range(max_iter): # Compute gradient and Hessian grad = grad_f(x) grad_norm = np.linalg.norm(grad) f_val = f(x) history.append((x.copy(), f_val, grad_norm)) if verbose: print(f"Iter {k:3d}: f = {f_val:.6e}, ||∇f|| = {grad_norm:.6e}") # Check convergence if grad_norm < tol: if verbose: print(f"Converged in {k} iterations") break # Compute Hessian H = hess_f(x) # Solve Newton system: H @ d = -grad # Using Cholesky factorization (assumes H is positive definite) try: L, lower = cho_factor(H) d = cho_solve((L, lower), -grad) except np.linalg.LinAlgError: # Hessian not positive definite, fall back to general solver d = solve(H, -grad) # Newton update (full step) x = x + d return x, history # Example: Minimize Rosenbrock function def rosenbrock(x): """Rosenbrock function: f(x,y) = (a-x)² + b(y-x²)²""" a, b = 1.0, 100.0 return (a - x[0])**2 + b * (x[1] - x[0]**2)**2 def rosenbrock_grad(x): """Gradient of Rosenbrock function""" a, b = 1.0, 100.0 dx = -2*(a - x[0]) - 4*b*x[0]*(x[1] - x[0]**2) dy = 2*b*(x[1] - x[0]**2) return np.array([dx, dy]) def rosenbrock_hess(x): """Hessian of Rosenbrock function""" a, b = 1.0, 100.0 h11 = 2 - 4*b*x[1] + 12*b*x[0]**2 h12 = -4*b*x[0] h22 = 2*b return np.array([[h11, h12], [h12, h22]]) if __name__ == "__main__": # Starting point far from minimum at (1, 1) x0 = np.array([-1.0, 1.0]) x_opt, history = newtons_method( rosenbrock, rosenbrock_grad, rosenbrock_hess, x0, verbose=True ) print(f" Optimum found at: {x_opt}") print(f"Function value: {rosenbrock(x_opt):.2e}")6b:T8bf,def modified_newton_direction( grad: np.ndarray, H: np.ndarray, min_eigval: float = 1e-8 ) -> np.ndarray: """ Compute modified Newton direction with Hessian regularization. If H is not positive definite, add τI to make it so. Parameters: ----------- grad : np.ndarray Gradient vector (n,) H : np.ndarray Hessian matrix (n, n) min_eigval : float Minimum eigenvalue threshold for positive definiteness Returns: -------- d : np.ndarray Modified Newton direction (n,) """ n = len(grad) # Compute eigendecomposition eigenvalues, eigenvectors = np.linalg.eigh(H) # Find minimum eigenvalue lambda_min = eigenvalues.min() if lambda_min > min_eigval: # H is already positive definite tau = 0.0 else: # Regularize to ensure positive definiteness tau = -lambda_min + min_eigval # Regularized Hessian H_reg = H + tau * np.eye(n) # Solve modified Newton system d = np.linalg.solve(H_reg, -grad) return d def trust_region_cauchy_point( grad: np.ndarray, H: np.ndarray, delta: float ) -> np.ndarray: """ Compute Cauchy point for trust region subproblem. The Cauchy point minimizes the quadratic model along the steepest descent direction within the trust region. Parameters: ----------- grad : np.ndarray Gradient vector (n,) H : np.ndarray Hessian matrix (n, n) delta : float Trust region radius Returns: -------- p_c : np.ndarray Cauchy point (n,) """ grad_norm = np.linalg.norm(grad) if grad_norm < 1e-12: return np.zeros_like(grad) # Normalized steepest descent direction g_hat = grad / grad_norm # Curvature along steepest descent gHg = np.dot(grad, H @ grad) if gHg <= 0: # Negative or zero curvature: go to boundary tau_star = delta / grad_norm else: # Positive curvature: unconstrained minimizer or boundary tau_unc = grad_norm**2 / gHg tau_star = min(tau_unc, delta / grad_norm) return -tau_star * grad6c:Ta9a,def backtracking_line_search( f: Callable[[np.ndarray], float], x: np.ndarray, d: np.ndarray, grad: np.ndarray, alpha_init: float = 1.0, c1: float = 1e-4, rho: float = 0.5, max_iter: int = 50 ) -> float: """ Backtracking line search satisfying Armijo condition. Finds step size α such that: f(x + α·d) ≤ f(x) + c1·α·(∇f)ᵀd Parameters: ----------- f : callable Objective function x : np.ndarray Current point d : np.ndarray Search direction (should be descent direction) grad : np.ndarray Gradient at x alpha_init : float Initial step size (1.0 for Newton methods) c1 : float Armijo parameter (typically 1e-4) rho : float Backtracking factor (typically 0.5) max_iter : int Maximum backtracking iterations Returns: -------- alpha : float Step size satisfying Armijo condition """ alpha = alpha_init f_x = f(x) grad_d = np.dot(grad, d) # Directional derivative if grad_d >= 0: raise ValueError("d is not a descent direction (∇fᵀd ≥ 0)") for _ in range(max_iter): if f(x + alpha * d) <= f_x + c1 * alpha * grad_d: return alpha alpha *= rho return alpha def damped_newton_method( f: Callable[[np.ndarray], float], grad_f: Callable[[np.ndarray], np.ndarray], hess_f: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, tol: float = 1e-8, max_iter: int = 100, verbose: bool = False ) -> Tuple[np.ndarray, list]: """ Damped Newton's method with backtracking line search. Combines global convergence of line search with local quadratic convergence of Newton's method. """ x = x0.copy().astype(float) history = [] for k in range(max_iter): grad = grad_f(x) grad_norm = np.linalg.norm(grad) f_val = f(x) history.append((x.copy(), f_val, grad_norm)) if verbose: print(f"Iter {k:3d}: f = {f_val:.6e}, ||∇f|| = {grad_norm:.6e}") if grad_norm < tol: if verbose: print(f"Converged in {k} iterations") break # Compute Newton direction with regularization H = hess_f(x) d = modified_newton_direction(grad, H) # Line search for step size alpha = backtracking_line_search(f, x, d, grad) if verbose and alpha < 1.0: print(f" Step size reduced to α = {alpha:.4f}") # Damped Newton update x = x + alpha * d return x, history6d:T454,**Comparison with gradient descent:** For an optimization problem with n parameters requiring K iterations to converge: - **Gradient Descent**: O(K_GD × n) total cost, where K_GD can be O(1/ε) or worse - **Newton's Method**: O(K_N × n³) total cost, where K_N is often O(log log(1/ε)) The crossover point depends on: - **Problem dimension n**: For n > 10,000, the O(n³) cost dominates - **Condition number κ**: Gradient descent iterations scale with κ, Newton's do not - **Available Hessian structure**: Sparse or structured Hessians can reduce O(n³) to O(n) or O(n²) **When Newton's method wins:** 1. **Small to moderate dimensions (n < 1000)**: The per-iteration cost is manageable, and quadratic convergence delivers huge iteration savings 2. **Ill-conditioned problems**: Gradient descent struggles with high condition numbers, making many more iterations necessary 3. **High accuracy requirements**: When ε = 10⁻¹² is needed, Newton's O(log log(1/ε)) scaling is unbeatable 4. **Cheap Hessians**: Some problems (quadratics, generalized linear models) have analytically simple Hessians6e:T934,def logistic_regression_newton( X: np.ndarray, y: np.ndarray, max_iter: int = 20, tol: float = 1e-8 ) -> np.ndarray: """ Logistic regression using Newton's method (IRLS). For binary classification: P(y=1|x) = σ(wᵀx) = 1/(1+exp(-wᵀx)) Parameters: ----------- X : np.ndarray Feature matrix (N, d) y : np.ndarray Binary labels (N,) with values in {0, 1} max_iter : int Maximum Newton iterations tol : float Convergence tolerance on weight change Returns: -------- w : np.ndarray Optimal weights (d,) """ N, d = X.shape w = np.zeros(d) # Initialize weights for iteration in range(max_iter): # Compute predictions z = X @ w p = 1 / (1 + np.exp(-z)) # Sigmoid # Gradient of negative log-likelihood # ∇L = Xᵀ(p - y) gradient = X.T @ (p - y) # Hessian of negative log-likelihood # H = Xᵀ W X, where W = diag(p * (1-p)) W = p * (1 - p) # Diagonal weights # Efficient Hessian-vector product # H = Xᵀ diag(W) X XtW = X.T * W # Broadcasting: (d, N) * (N,) = (d, N) H = XtW @ X # (d, N) @ (N, d) = (d, d) # Add small regularization for numerical stability H += 1e-8 * np.eye(d) # Newton direction: d = -H⁻¹g # Solve: H @ delta = -gradient delta = np.linalg.solve(H, -gradient) # Update weights w_new = w + delta # Check convergence if np.linalg.norm(delta) < tol: print(f"Converged in {iteration + 1} iterations") return w_new w = w_new print(f"Max iterations ({max_iter}) reached") return w # Example usage if __name__ == "__main__": from sklearn.datasets import make_classification # Generate synthetic binary classification data X, y = make_classification( n_samples=1000, n_features=20, n_informative=10, n_redundant=5, random_state=42 ) w_optimal = logistic_regression_newton(X, y) # Check accuracy predictions = (X @ w_optimal > 0).astype(int) accuracy = np.mean(predictions == y) print(f"Training accuracy: {accuracy:.2%}")6f:T40a,**Gaussian Process Posterior:** GP inference requires solving K⁻¹y where K is the kernel matrix. Newton's method efficiently optimizes hyperparameters by leveraging closed-form gradients and Hessians of the marginal likelihood. **Natural Gradient Methods:** In probabilistic models, the **natural gradient** preconditions the standard gradient with the Fisher information matrix (expected Hessian under the current distribution). This achieves approximate Newton-like adaptation without computing the true Hessian. **Second-order optimization for neural networks:** Recent advances approximate the Hessian efficiently: - **K-FAC (Kronecker-Factored Approximate Curvature)**: Approximates layer-wise Hessians using Kronecker products - **Shampoo**: Uses matrix-valued preconditioners in the spirit of Newton's method - **Hessian-free optimization**: Computes Hessian-vector products without forming the full Hessian These methods achieve better convergence than SGD for certain problems while maintaining tractable complexity.63:{"json":{"id":"seed_129b6de46a4416acb104338467c6e832","slug":"advanced-optimization","category":"Machine Learning","topic":"Advanced Optimization","title":"Advanced Optimization Methods","description":"Master sophisticated optimization algorithms that go beyond gradient descent, including Newton's method, quasi-Newton methods, coordinate descent, and proximal methods—essential tools for training complex machine learning models efficiently.","difficulty":"ADVANCED","estimatedMinutes":90,"pages":[{"id":"seed_589a6b98aff004e39ca82a7c1e99bdb6","title":"Newton's Method","order":0,"content":[{"id":"opening-hook","title":"Beyond First-Order Optimization","content":[{"data":{"content":"Gradient descent, while foundational, is a **first-order method**—it uses only gradient information (first derivatives) to navigate the loss landscape. This limitation manifests as slow convergence near optima, sensitivity to learning rates, and difficulty with ill-conditioned problems.\n\n**Newton's method** fundamentally changes the optimization paradigm by incorporating **second-order information**—the curvature of the objective function. Rather than taking small steps proportional to the gradient, Newton's method uses a local quadratic approximation to jump directly toward the minimum. This single insight transforms optimization from iterative trial to geometric precision."},"type":"text"},{"data":{"icon":"Target","title":"What You Will Learn","content":"By the end of this page, you will understand the mathematical derivation of Newton's method from Taylor expansion, analyze its quadratic convergence properties, implement the algorithm from first principles, and critically evaluate when Newton's method excels versus when simpler methods suffice.","variant":"info"},"type":"callout"}]},{"id":"historical-foundation","title":"Historical and Conceptual Foundation","content":[{"data":{"content":"Newton's method for optimization is an extension of **Newton-Raphson iteration** for root finding, first described by Isaac Newton in 1669 and later refined by Joseph Raphson in 1690. The connection is elegant: finding where a function achieves its minimum is equivalent to finding where its gradient equals zero.\n\n**The root-finding perspective:**\n\nIf we want to minimize a scalar function f(x), we seek x* where f'(x*) = 0. Newton-Raphson finds roots of g(x) = 0 by iterating:\n\nx_{n+1} = x_n - g(x_n) / g'(x_n)\n\nApplying this to g(x) = f'(x), we obtain:\n\nx_{n+1} = x_n - f'(x_n) / f''(x_n)\n\nThis is Newton's method for scalar optimization. The multivariate generalization replaces the scalar derivative with the gradient and the second derivative with the **Hessian matrix**."},"type":"text"},{"data":{"icon":"Lightbulb","title":"The Geometric Intuition","content":"At each iteration, Newton's method fits a parabola (or paraboloid in higher dimensions) to the local shape of the function and jumps directly to that parabola's minimum. If the function is actually quadratic, Newton's method finds the exact minimum in one step. For non-quadratic functions, repeated quadratic approximations rapidly hone in on the true minimum.","variant":"tip"},"type":"callout"},{"data":{"content":"**Why curvature matters:**\n\nConsider two scenarios where the gradient has the same magnitude but different curvatures:\n\n1. **High curvature**: The function is changing rapidly. We're close to a minimum and should take a smaller step.\n2. **Low curvature**: The function is relatively flat. We're far from the minimum and should take a larger step.\n\nGradient descent cannot distinguish these scenarios—it takes steps proportional only to gradient magnitude. Newton's method naturally incorporates this distinction through the Hessian. When curvature is high (large Hessian eigenvalues), the Hessian inverse shrinks the step. When curvature is low (small eigenvalues), the step is amplified.\n\nThis **curvature-adaptive step sizing** is the fundamental innovation of second-order methods."},"type":"text"}],"titleType":"number","titleNumber":1},{"id":"mathematical-derivation","title":"Mathematical Derivation from Taylor Expansion","content":[{"data":{"content":"The rigorous derivation of Newton's method proceeds from the second-order Taylor expansion of the objective function. Let f: ℝⁿ → ℝ be a twice continuously differentiable function. Around a point x_k, the Taylor expansion is:\n\nf(x_k + Δx) ≈ f(x_k) + ∇f(x_k)ᵀΔx + ½Δxᵀ H(x_k) Δx\n\nwhere:\n- ∇f(x_k) ∈ ℝⁿ is the gradient at x_k\n- H(x_k) ∈ ℝⁿˣⁿ is the Hessian matrix at x_k\n- Δx ∈ ℝⁿ is the step to take\n\nThis quadratic approximation defines a **local quadratic model** m_k(Δx) of the objective function."},"type":"text"},{"data":{"code":"# Local quadratic model at point x_k\ndef quadratic_model(delta_x, f_k, grad_k, hessian_k):\n \"\"\"\n Computes the quadratic approximation:\n m_k(Δx) = f(x_k) + ∇f(x_k)ᵀΔx + ½Δxᵀ H(x_k) Δx\n \n Parameters:\n -----------\n delta_x : np.ndarray\n Step direction (n,)\n f_k : float\n Function value at x_k\n grad_k : np.ndarray\n Gradient at x_k (n,)\n hessian_k : np.ndarray\n Hessian matrix at x_k (n, n)\n \n Returns:\n --------\n float : Quadratic model value\n \"\"\"\n linear_term = np.dot(grad_k, delta_x)\n quadratic_term = 0.5 * np.dot(delta_x, hessian_k @ delta_x)\n return f_k + linear_term + quadratic_term","fileName":"quadratic_model.py","showLineNumbers":true},"type":"code"},{"data":{"content":"**Minimizing the quadratic model:**\n\nTo find the optimal step Δx*, we set the gradient of the quadratic model to zero:\n\n∇_{Δx} m_k(Δx) = ∇f(x_k) + H(x_k)Δx = 0\n\nSolving for Δx:\n\n**Δx* = -H(x_k)⁻¹ ∇f(x_k)**\n\nThis is the **Newton direction**. The full Newton update is:\n\n**x_{k+1} = x_k - H(x_k)⁻¹ ∇f(x_k)**\n\nThis formula encapsulates the essence of Newton's method: the step is the gradient, preconditioned (rescaled and rotated) by the inverse Hessian."},"type":"text"},{"data":{"icon":"Info","title":"The Hessian's Role","content":"The inverse Hessian H⁻¹ acts as a **variable metric**—it stretches space along directions of low curvature and compresses it along directions of high curvature. This normalization makes Newton's method invariant to affine transformations of the coordinate system, a property first-order methods lack.","variant":"note"},"type":"callout"},{"data":{"content":"**Second-order optimality conditions:**\n\nFor the Newton direction to be a descent direction (moving toward a minimum rather than a maximum), the Hessian must be **positive definite**. Recall:\n\n- A matrix H is positive definite if vᵀHv > 0 for all nonzero v\n- This means all eigenvalues of H are positive\n- Geometrically, the function curves upward in all directions (local convexity)\n\nIf H is not positive definite at x_k, the Newton direction may point toward a maximum or saddle point. This is a fundamental challenge addressed by modifications we'll discuss later."},"type":"text"}],"titleType":"number","titleNumber":2},{"id":"convergence-analysis","title":"Convergence Analysis: Quadratic Convergence","content":[{"data":{"content":"The remarkable property of Newton's method is its **quadratic convergence rate** near a solution. This means the number of correct digits roughly doubles at each iteration—exponentially faster than the linear convergence of gradient descent.\n\n**Formal statement (Local Convergence Theorem):**\n\nLet f: ℝⁿ → ℝ be twice continuously differentiable with Lipschitz continuous Hessian near a local minimum x*. If:\n1. ∇f(x*) = 0 (first-order necessary condition)\n2. H(x*) is positive definite (second-order sufficient condition)\n3. x₀ is sufficiently close to x*\n\nThen Newton's method converges quadratically:\n\n‖x_{k+1} - x*‖ ≤ C · ‖x_k - x*‖²\n\nfor some constant C > 0."},"type":"text"},{"data":{"rows":[["Gradient Descent","Linear: O(ρᵏ)","Constant factor ρ < 1","~120 (for ρ = 0.9)"],["Accelerated GD (Nesterov)","Linear: O(ρ^√k)","Improved constant","~50-80"],["Newton's Method","Quadratic: O(ε²ᵏ)","Error squared each step","~4-6 (near convergence)"],["Conjugate Gradient","Super-linear","Depends on spectrum","~20-40"]],"caption":"Convergence Rate Comparison","headers":["Method","Convergence Rate","Error Reduction per Iteration","Iterations for 10⁻¹² Error"]},"type":"table"},{"data":{"content":"**Understanding quadratic convergence:**\n\nConsider a concrete example. Suppose at iteration k, the error ‖x_k - x*‖ = 10⁻². With quadratic convergence:\n\n- Iteration k+1: error ≈ (10⁻²)² = 10⁻⁴\n- Iteration k+2: error ≈ (10⁻⁴)² = 10⁻⁸\n- Iteration k+3: error ≈ (10⁻⁸)² = 10⁻¹⁶\n\nThree iterations take us from 1% error to machine precision! This dramatic acceleration explains why Newton's method, despite its per-iteration cost, often achieves convergence in far fewer total operations.\n\n**The catch—local convergence:**\n\nThe quadratic rate only holds when starting \"sufficiently close\" to x*. From far away, Newton's method may:\n- Diverge entirely\n- Oscillate without converging\n- Converge to a different stationary point (maximum or saddle)\n\nThis **local convergence** property motivates numerous modifications and globalizations."},"type":"text"},{"data":{"icon":"AlertCircle","title":"The Convergence Region","content":"The \"sufficiently close\" requirement is not just theoretical pedantry. In practice, pure Newton's method can fail spectacularly from poor starting points. The convergence basin can be surprisingly small for non-convex functions, and identifying whether you're inside it is generally impossible without already knowing x*.","variant":"warning"},"type":"callout"}],"titleType":"number","titleNumber":3},{"id":"algorithm-implementation","title":"Algorithm and Implementation","content":[{"data":{"content":"Implementing Newton's method correctly requires careful attention to numerical stability, especially in computing the Newton direction. We'll develop the algorithm from first principles.\n\n**The Pure Newton Algorithm:**\n\n```\nInput: Starting point x₀, tolerance ε, max iterations K\nOutput: Approximate minimizer x*\n\n1. k ← 0\n2. While ‖∇f(x_k)‖ > ε and k < K:\n a. Compute gradient g_k = ∇f(x_k)\n b. Compute Hessian H_k = ∇²f(x_k)\n c. Solve linear system: H_k d_k = -g_k for Newton direction d_k\n d. Update: x_{k+1} = x_k + d_k\n e. k ← k + 1\n3. Return x_k\n```\n\nThe critical computational step is solving the linear system in step 2c. **Never compute H⁻¹ explicitly**—this is numerically unstable and computationally wasteful."},"type":"text"},{"data":{"code":"$6a","fileName":"newtons_method.py","showLineNumbers":true},"type":"code"},{"data":{"icon":"Lightbulb","title":"Numerical Best Practices","content":"Always solve H·d = -g rather than computing d = H⁻¹·(-g). Use Cholesky factorization (cho_solve) when H is positive definite—it's twice as fast as LU and numerically more stable. For indefinite H, use LDL factorization or regularization techniques.","variant":"tip"},"type":"callout"}],"titleType":"number","titleNumber":4},{"id":"handling-indefinite-hessians","title":"Handling Non-Positive-Definite Hessians","content":[{"data":{"content":"When the Hessian is not positive definite, the Newton direction may not be a descent direction, leading to convergence toward saddle points or maxima. Several strategies address this challenge:\n\n**1. Hessian Modification (Regularization):**\n\nAdd a multiple of the identity matrix to ensure positive definiteness:\n\nH̃_k = H_k + τI\n\nwhere τ > 0 is chosen large enough that H̃_k is positive definite. This is equivalent to adding L2 regularization to the local quadratic model.\n\n**Selection strategies for τ:**\n- Fixed small τ (simple but may over-regularize)\n- Adaptive τ based on smallest eigenvalue: τ = max(0, -λ_min + δ) for small δ > 0\n- Multiple of gradient norm: τ = β‖∇f(x_k)‖\n\n**2. Modified Cholesky Factorization:**\n\nMore sophisticated algorithms (like the Gill-Murray modification) compute a modified Cholesky factorization directly, perturbing the factorization minimally to ensure positive definiteness while preserving the Cholesky structure."},"type":"text"},{"data":{"code":"$6b","fileName":"modified_newton.py","showLineNumbers":true},"type":"code"},{"data":{"content":"**3. Trust Region Methods:**\n\nRather than modifying the Hessian directly, trust region methods constrain the Newton step to a region where the quadratic model is trusted:\n\nminimize m_k(p) = f_k + ∇f_kᵀp + ½pᵀH_kp\nsubject to ‖p‖ ≤ Δ_k\n\nwhere Δ_k is the trust region radius. This subproblem always has a well-defined solution, even when H_k is indefinite. The trust region radius adapts based on how well the quadratic model predicts actual function decrease.\n\n**4. Line Search with Negative Curvature Direction:**\n\nWhen H_k has a negative eigenvalue, we can extract the corresponding eigenvector and use it as a descent direction. The function decreases along directions of negative curvature, potentially escaping saddle points."},"type":"text"},{"data":{"icon":"Info","title":"Practice vs Theory","content":"For strictly convex problems (like linear regression with regularization), the Hessian is always positive definite. For neural networks and other non-convex models, indefinite Hessians are the norm. Production optimizers like L-BFGS inherently maintain positive definite approximations, sidestepping this issue entirely.","variant":"note"},"type":"callout"}],"titleType":"number","titleNumber":5},{"id":"damped-newton","title":"Damped Newton and Line Search","content":[{"data":{"content":"Pure Newton's method takes full steps in the Newton direction, which can lead to divergence from poor starting points. **Damped Newton methods** introduce a step size α_k ∈ (0, 1]:\n\nx_{k+1} = x_k + α_k d_k\n\nwhere d_k = -H_k⁻¹∇f_k is the Newton direction.\n\n**Backtracking line search (Armijo rule):**\n\nThe most common approach finds the largest α ≤ 1 satisfying the **sufficient decrease condition**:\n\nf(x_k + αd_k) ≤ f(x_k) + c₁ α ∇f_kᵀd_k\n\nwhere c₁ ∈ (0, 1) is typically 10⁻⁴. This ensures the function actually decreases, providing global convergence guarantees while preserving fast local convergence."},"type":"text"},{"data":{"code":"$6c","fileName":"damped_newton.py","showLineNumbers":true},"type":"code"},{"data":{"content":"**Convergence properties of damped Newton:**\n\nWith appropriate line search:\n1. **Global convergence**: From any starting point, the sequence converges to a stationary point\n2. **Local quadratic convergence**: Near a strict local minimum with positive definite Hessian, the step size becomes α = 1, recovering pure Newton's quadratic rate\n3. **Transition**: The algorithm seamlessly transitions from \"careful\" globally convergent steps to \"aggressive\" locally optimal steps\n\nThis combination makes damped Newton practical for real-world optimization where initialization quality cannot be guaranteed."},"type":"text"}],"titleType":"number","titleNumber":6},{"id":"computational-complexity","title":"Computational Complexity and Cost Analysis","content":[{"data":{"content":"Newton's method's superior convergence rate comes at significant computational cost. Understanding when this tradeoff is favorable requires careful complexity analysis.\n\n**Per-iteration costs:**"},"type":"text"},{"data":{"rows":[["Gradient computation","O(n) to O(n·m)","Small n, large data (m samples)"],["Hessian computation","O(n²) to O(n²·m)","Moderate n"],["Hessian storage","O(n²) memory","Always for dense Hessian"],["Linear system solve","O(n³) time","Large n (thousands)"],["Total per iteration","O(n³) or O(n²·m)","Depends on problem structure"]],"caption":"Computational Cost Breakdown","headers":["Operation","Cost (n parameters)","Dominant When"]},"type":"table"},{"data":{"content":"$6d"},"type":"text"},{"data":{"icon":"Lightbulb","title":"The Deep Learning Exception","content":"Neural networks typically have n > 10⁶ parameters, making O(n²) Hessian storage alone prohibitive. This explains why deep learning universally uses first-order methods (SGD, Adam). Quasi-Newton methods like L-BFGS bridge the gap by approximating second-order information with O(n) storage, as we'll see in the next page.","variant":"tip"},"type":"callout"}],"titleType":"number","titleNumber":7},{"id":"newton-in-ml","title":"Newton's Method in Machine Learning","content":[{"data":{"content":"Despite its cubic complexity, Newton's method remains valuable in specific machine learning contexts where its advantages outweigh computational costs.\n\n**Logistic Regression (IRLS):**\n\nLogistic regression is typically trained using **Iteratively Reweighted Least Squares (IRLS)**, which is exactly Newton's method for maximum likelihood estimation. For N samples with d features:\n\n- Hessian: H = XᵀWX where W is diagonal (N×N)\n- Cost: O(Nd² + d³) per iteration\n- Convergence: 5-10 iterations typical\n\nFor d << N (typical), this is highly efficient."},"type":"text"},{"data":{"code":"$6e","fileName":"logistic_newton.py","showLineNumbers":true},"type":"code"},{"data":{"content":"$6f"},"type":"text"},{"data":{"icon":"Info","title":"The Gauss-Newton Connection","content":"For least-squares problems f(x) = ½‖r(x)‖², the Gauss-Newton method approximates the Hessian as H ≈ JᵀJ (ignoring second derivatives of residuals). This simplification is exact when residuals are small near the optimum and forms the basis for Levenberg-Marquardt and trust-region methods widely used in nonlinear regression.","variant":"note"},"type":"callout"}],"titleType":"number","titleNumber":8},{"id":"summary-key-takeaways","title":"Summary and Key Takeaways","content":[{"data":{"content":"Newton's method represents a fundamental paradigm in optimization—exploiting second-order information to achieve rapid convergence. While its computational cost limits direct application to large-scale ML, it provides the theoretical foundation for many practical algorithms."},"type":"text"},{"data":{"type":"ordered","items":["**Newton's method uses second-order information** — The Hessian matrix captures curvature, enabling curvature-adaptive step sizes that dramatically accelerate convergence.","**Quadratic convergence near optima** — Error squares at each iteration, reaching machine precision in just a few steps once in the convergence basin.","**Requires positive definite Hessian** — For non-convex problems, modifications like regularization, trust regions, or line search ensure descent.","**O(n³) per-iteration cost** — Computing and factorizing the Hessian is prohibitive for large n, motivating quasi-Newton approximations.","**Ideal for small-to-medium convex problems** — Logistic regression, GP hyperparameter optimization, and other structured problems benefit from Newton's method.","**Foundation for advanced methods** — Understanding Newton's method is essential for comprehending L-BFGS, natural gradients, and second-order deep learning optimizers."],"title":{"icon":"CheckCircle","text":"Key Takeaways","type":"icon","color":"primary","variant":"outline-circle"},"variant":"accent"},"type":"list"},{"data":{"content":"**What's next:**\n\nThe next page explores **Quasi-Newton methods**—especially BFGS and L-BFGS—which achieve Newton-like convergence without computing the Hessian explicitly. These methods approximate the inverse Hessian using only gradient information, making second-order optimization practical for problems with thousands of parameters."},"type":"text"},{"data":{"icon":"CheckCircle","title":"Page Complete","content":"You now understand Newton's method: its derivation from Taylor expansion, quadratic convergence properties, implementation with line search and regularization, and application to machine learning problems. This foundation prepares you for quasi-Newton methods that scale these ideas to larger problems.","variant":"success"},"type":"callout"}],"titleType":"number","titleNumber":9}],"isLocked":false},{"id":"seed_096687b06e5533b1372eeb8ae5616c71","title":"Quasi-Newton Methods (BFGS, L-BFGS)","order":1,"content":null,"isLocked":true},{"id":"seed_d959397675821e4e63a0cdd5679cef77","title":"Coordinate Descent","order":2,"content":null,"isLocked":true},{"id":"seed_c01022a29f9e92d7e5134127d547526b","title":"Proximal Methods","order":3,"content":null,"isLocked":true},{"id":"seed_cc611f28262e9e8cb78954194121d39b","title":"Optimization Landscapes","order":4,"content":null,"isLocked":true}]}}