Multivariate Optimization: Gradient Descent to Newton-Raphson

Part 26 of 12-Week BusinessMath Series

What You’ll Learn

• Understanding numerical differentiation for gradients and Hessians
• Using gradient descent with momentum and Nesterov acceleration
• Applying Newton-Raphson and BFGS quasi-Newton methods
• Choosing the right optimization algorithm for your problem
• Using AdaptiveOptimizer for automatic algorithm selection
• Understanding convergence rates and algorithm trade-offs

The Problem

Real-world optimization problems have multiple variables:

• Parameter fitting: Minimize error across 10+ parameters
• Cost minimization: Optimize production mix across multiple products/facilities
• Portfolio construction: Find optimal weights for N assets
• Machine learning: Fit models with thousands of parameters

Single-variable methods (like goal-seeking) don’t extend to multivariate problems—you need algorithms designed for N dimensions.

The Solution

BusinessMath provides a progression of multivariate optimizers, from simple gradient descent to sophisticated second-order methods. All work generically with any VectorSpace type.

Numerical Differentiation

When you can’t compute derivatives analytically, BusinessMath computes them numerically:

import BusinessMath

// Define f(x,y) = x² + 2y²
let function: (VectorN) -> Double = { v in
    let x = v[0]
    let y = v[1]
    return x*x + 2*y*y
}

// Compute gradient at (1, 2)
let point = VectorN([1.0, 2.0])
let gradient = try numericalGradient(function, at: point)
print("Gradient: \(gradient.toArray())")  // ≈ [2.0, 8.0]
// Analytical: ∂f/∂x = 2x = 2, ∂f/∂y = 4y = 8 ✓

// Compute Hessian (curvature matrix)
let hessian = try numericalHessian(function, at: point)
print("Hessian:")
for row in hessian {
    print(row.map { $0.number(1) })
}
// [[2.0, 0.0], [0.0, 4.0]]

Output:

Gradient: [2.0, 8.0]
Hessian:
[2.0, 0.0]
[0.0, 4.0]

How it works:

• Gradient: Central finite differences (f(x+h) - f(x-h)) / 2h
• Hessian: Second-order finite differences (N² function evaluations)

Gradient Descent: The Workhorse

Gradient descent iteratively moves in the direction of steepest descent:

import BusinessMath

// Minimize f(x,y) = x² + 4y²
let function: (VectorN) -> Double = { v in
    v[0]*v[0] + 4*v[1]*v[1]
}

// Basic gradient descent
let optimizer = MultivariateGradientDescent>(
    learningRate: 0.01,
    maxIterations: 1000,
    tolerance: 1e-6
)

let result = try optimizer.minimize(
    function: function,
    gradient: { x in try numericalGradient(function, at: x) },
    initialGuess: VectorN([5.0, 5.0])
)

print("Minimum at: \(result.solution.toArray().map({ $0.number(3) }))")
print("Value: \(result.objectiveValue.number(6))")
print("Iterations: \(result.iterations)")
print("Converged: \(result.converged)")

Output:

Minimum at: [0.0, 0.0]
Value: 0.000000
Iterations: 247
Converged: true

Gradient Descent with Momentum

Momentum accelerates convergence and reduces oscillations:

import BusinessMath

// Rosenbrock function: classic test problem
let rosenbrock: (VectorN) -> Double = { v in
    let x = v[0], y = v[1]
    let a = 1 - x
    let b = y - x*x
    return a*a + 100*b*b  // Minimum at (1, 1)
}

// Gradient descent with momentum (default 0.9)
let optimizerWithMomentum = GradientDescentOptimizer(
	learningRate: 0.01,
	maxIterations: 5000,
	momentum: 0.9,
	useNesterov: false
)

// Note: Using scalar optimizer for demonstration
// For VectorN, use MultivariateGradientDescent

let result = optimizerWithMomentum.optimize(
    objective: { x in (x - 5) * (x - 5) },
    constraints: [],
    initialGuess: 0.0,
    bounds: nil
)

print("Converged to: \(result.optimalValue.number(4))")
print("Iterations: \(result.iterations)")

Nesterov Acceleration (look-ahead gradient) often converges even faster:

let nesterovOptimizer = GradientDescentOptimizer(
    learningRate: 0.01,
    momentum: 0.9,
    useNesterov: true  // Nesterov acceleration
)

Newton-Raphson: Quadratic Convergence

Newton-Raphson uses second-order information (Hessian) for much faster convergence:

import BusinessMath

// Quadratic function: f(x,y) = x² + 4y² + 2xy
let quadratic: (VectorN) -> Double = { v in
    let x = v[0], y = v[1]
    return x*x + 4*y*y + 2*x*y
}

// Full Newton-Raphson (uses exact Hessian)
let newtonOptimizer = MultivariateNewtonRaphson>(
    maxIterations: 100,
    tolerance: 1e-8,
    useLineSearch: true
)

let result = try newtonOptimizer.minimize(
    function: quadratic,
    gradient: { try numericalGradient(quadratic, at: $0) },
    hessian: { try numericalHessian(quadratic, at: $0) },
    initialGuess: VectorN([10.0, 10.0])
)

print("Solution: \(result.solution.toArray().map({ $0.number(6) }))")
print("Converged in: \(result.iterations) iterations")

Output:

Solution: [0.000000, 0.000000]
Converged in: 3 iterations

The power: Newton-Raphson found the minimum in 3 iterations vs. 247 for gradient descent!

BFGS: Quasi-Newton Sweet Spot

BFGS approximates the Hessian, giving Newton-like speed without expensive Hessian computation:

import BusinessMath

let rosenbrock: (VectorN) -> Double = { v in
    let x = v[0], y = v[1]
    let a = 1 - x
    let b = y - x*x
    return a*a + 100*b*b
}

// BFGS quasi-Newton
let bfgsOptimizer = MultivariateNewtonRaphson>()

let result = try bfgsOptimizer.minimizeBFGS(
    function: rosenbrock,
    gradient: { try numericalGradient(rosenbrock, at: $0) },
    initialGuess: VectorN([0.0, 0.0])
)

print("Solution: \(result.solution.toArray().map({ $0.rounded(toPlaces: 4) }))")
print("Iterations: \(result.iterations)")
print("Final value: \(result.objectiveValue.rounded(toPlaces: 8))")

Output:

Solution: [1.0000, 1.0000]
Iterations: 24
Final value: 0.00000001

Comparison:

The trade-off: BFGS balances speed and computational cost—best for most practical problems.

AdaptiveOptimizer: Automatic Algorithm Selection

Don’t know which algorithm to use? Let AdaptiveOptimizer decide:

import BusinessMath

// AdaptiveOptimizer chooses the best algorithm automatically
let optimizer = AdaptiveOptimizer>()

let rosenbrock: (VectorN) -> Double = { v in
    let x = v[0], y = v[1]
    return (1-x)*(1-x) + 100*(y-x*x)*(y-x*x)
}

let result = try optimizer.optimize(
    objective: rosenbrock,
    initialGuess: VectorN([0.0, 0.0]),
    constraints: []
)

print("Solution: \(result.solution.toArray().map({ $0.rounded(toPlaces: 4) }))")
print("Algorithm used: \(result.algorithmUsed ?? "N/A")")
print("Reason: \(result.selectionReason ?? "N/A")")

Output:

Solution: [1.0000, 1.0000]
Algorithm used: Newton-Raphson
Reason: Small problem (2 variables) - using Newton-Raphson for fast convergence

How it works:

• Analyzes problem size, constraints, gradient availability
• Selects: Gradient Descent, Newton-Raphson, BFGS, or Constrained optimizer
• Reports which algorithm was chosen and why

Choosing the Right Algorithm

Rule of thumb:

• < 100 variables + smooth: Use BFGS
• 100-10,000 variables: Use Momentum/Nesterov
• > 10,000 variables: Use basic Gradient Descent
• Constraints: Use ConstrainedOptimizer or InequalityOptimizer (Week 8 Wednesday)

Real-World Example: Parameter Fitting

Fit a curve to noisy data:

import BusinessMath

// Data: y = a*x² + b*x + c + noise
let xData = VectorN.linearSpace(from: 0.0, to: 10.0, count: 50)
let yData = xData.map { x in
    2.0 * x * x + 3.0 * x + 5.0 + Double.random(in: -5...5)
}

// Objective: Minimize sum of squared errors
let objective: (VectorN) -> Double = { params in
    let a = params[0], b = params[1], c = params[2]
    var sse = 0.0
    for i in 0..>()
let result = try optimizer.minimizeBFGS(
    function: objective,
    gradient: { try numericalGradient(objective, at: $0) },
	initialGuess: VectorN([1.0, 2.0, 3.0])
)

print("Fitted parameters:")
print("  a = \(result_params.solution[0].number(2)) (true: 2.0)")
print("  b = \(result_params.solution[1].number(2)) (true: 3.0)")
print("  c = \(result_params.solution[2].number(2)) (true: 5.0)")
print("SSE: \(result_params.objectiveValue.number(1))")

Output:

Fitted parameters:
  a = 1.98 (true: 2.0)
  b = 3.17 (true: 3.0)
  c = 4.82 (true: 5.0)
SSE: 311.7

Try It Yourself

import BusinessMath

// MARK: - Numerical Differentiation

// Define f(x,y) = x² + 2y²
let function_nd: (VectorN) -> Double = { v in
	let x = v[0]
	let y = v[1]
	return x*x + 2*y*y
}

// Compute gradient at (1, 2)
let point_nd = VectorN([1.0, 2.0])
let gradient_nd = try numericalGradient(function_nd, at: point_nd)
print("Gradient: \(gradient_nd.toArray())")  // ≈ [2.0, 8.0]
// Analytical: ∂f/∂x = 2x = 2, ∂f/∂y = 4y = 8 ✓

// Compute Hessian (curvature matrix)
let hessian = try numericalHessian(function_nd, at: point_nd)
print("Hessian:")
for row in hessian {
	print(row.map { $0.number(1) })
}
// [[2.0, 0.0], [0.0, 4.0]]

// MARK: - Gradient Descent

// Minimize f(x,y) = x² + 4y²
let function_gd: (VectorN) -> Double = { v in
	v[0]*v[0] + 4*v[1]*v[1]
}

// Basic gradient descent
let optimizer_gd = MultivariateGradientDescent>(
	learningRate: 0.01,
	maxIterations: 1000,
	tolerance: 1e-6
)

let result_gd = try optimizer_gd.minimize(
	function: function_gd,
	gradient: { x in try numericalGradient(function_gd, at: x) },
	initialGuess: VectorN([5.0, 5.0])
)

print("Minimum at: \(result_gd.solution.toArray().map({ $0.number(3)  }))")
print("Value: \(result_gd.objectiveValue.number(6))")
print("Iterations: \(result_gd.iterations)")
print("Converged: \(result_gd.converged)")

// MARK: Gradient Descent with Momentum

// Rosenbrock function: classic test problem
let rosenbrock: (VectorN) -> Double = { v in
	let x = v[0], y = v[1]
	let a = 1 - x
	let b = y - x*x
	return a*a + 100*b*b  // Minimum at (1, 1)
}

// Gradient descent with momentum (default 0.9)
let optimizerWithMomentum = GradientDescentOptimizer(
	learningRate: 0.01,
	maxIterations: 5000,
	momentum: 0.9,
	useNesterov: false
)

// Note: Using scalar optimizer for demonstration
// For VectorN, use MultivariateGradientDescent

let result_gdm = optimizerWithMomentum.optimize(
	objective: { x in (x - 5) * (x - 5) },
	constraints: [],
	initialGuess: 0.0,
	bounds: nil
)

print("Converged to: \(result_gdm.optimalValue.number(1))")
print("Iterations: \(result_gdm.iterations)")

// MARK: - Newton-Raphson: Quadratic Convergence

	// Quadratic function: f(x,y) = x² + 4y² + 2xy
	let quadratic: (VectorN) -> Double = { v in
		let x = v[0], y = v[1]
		return x*x + 4*y*y + 2*x*y
	}

	// Full Newton-Raphson (uses exact Hessian)
	let newtonOptimizer = MultivariateNewtonRaphson>(
		maxIterations: 100,
		tolerance: 1e-8,
		useLineSearch: true
	)

	let result_newton = try newtonOptimizer.minimize(
		function: quadratic,
		gradient: { try numericalGradient(quadratic, at: $0) },
		hessian: { try numericalHessian(quadratic, at: $0) },
		initialGuess: VectorN([10.0, 10.0])
	)

	print("Solution: \(result_newton.solution.toArray().map({ $0.number(6) }))")
	print("Converged in: \(result_newton.iterations) iterations")

// MARK: - BFGS: Quasi-Newton Sweet Spot

// BFGS quasi-Newton
let bfgsOptimizer = MultivariateNewtonRaphson>()

let result_bfgs = try bfgsOptimizer.minimizeBFGS(
	function: rosenbrock,
	gradient: { try numericalGradient(rosenbrock, at: $0) },
	initialGuess: VectorN([0.0, 0.0])
)

print("Solution: \(result_bfgs.solution.toArray().map({ $0.number(4) }))")
print("Iterations: \(result_bfgs.iterations)")
print("Final value: \(result_bfgs.objectiveValue.number(8))")

// MARK: - Adaptive Optimizer

// AdaptiveOptimizer chooses the best algorithm automatically
let optimizer_adaptive = AdaptiveOptimizer>()

let result_adaptive = try optimizer_adaptive.optimize(
	objective: rosenbrock,
	initialGuess: VectorN([0.0, 0.0]),
	constraints: []
)

print("Solution: \(result_adaptive.solution.toArray().map({ $0.number(4) }))")
print("Algorithm used: \(result_adaptive.algorithmUsed)")
print("Reason: \(result_adaptive.selectionReason)")


// MARK: - Parameter Fitting Example

// Data: y = a*x² + b*x + c + noise
let xData = VectorN.linearSpace(from: 0.0, to: 10.0, count: 50)
let yData = xData.map { x in
	2.0 * x * x + 3.0 * x + 5.0 + Double.random(in: -5...5)
}

// Objective: Minimize sum of squared errors
let objective_params: (VectorN) -> Double = { params in
	let a = params[0], b = params[1], c = params[2]
	var sse = 0.0
	for i in 0..>()
let result_params = try optimizer_params.minimizeBFGS(
	function: objective_params,
	gradient: { try numericalGradient(objective_params, at: $0) },
	initialGuess: VectorN([1.0, 2.0, 3.0])
)

print("Fitted parameters:")
print("  a = \(result_params.solution[0].number(2)) (true: 2.0)")
print("  b = \(result_params.solution[1].number(2)) (true: 3.0)")
print("  c = \(result_params.solution[2].number(2)) (true: 5.0)")
print("SSE: \(result_params.objectiveValue.number(1))")

→ Full API Reference: BusinessMath Docs – 5.5 Multivariate Optimization

Modifications to try:

1. Compare convergence rates: plot iteration vs. objective value for each algorithm
2. Fit a 10-parameter model (polynomial, exponential, custom function)
3. Optimize a high-dimensional problem (100+ variables) with Momentum
4. Test robustness: start from different initial guesses and compare results

Real-World Application

• Machine learning: Train models by minimizing loss functions
• Engineering: Optimize design parameters (aerodynamics, materials, structures)
• Finance: Calibrate option pricing models to market data
• Operations: Optimize production schedules, routing, inventory

Data scientist use case: “I need to fit a pricing model with 15 parameters to historical transaction data. Manual tuning is infeasible—I need automated optimization.”

BFGS converges in seconds, not hours of manual tweaking.

★ Insight ─────────────────────────────────────

Why BFGS Beats Full Newton for Most Problems

Full Newton requires computing the Hessian (N² second derivatives). For a 100-variable problem:

• Full Newton: 10,000 Hessian elements per iteration
• BFGS: 0 Hessian evaluations (approximated from gradients)

BFGS maintains a Hessian approximation updated using gradient changes:

H_{k+1} = H_k + correction terms based on (∇f_{k+1} - ∇f_k)

Result: BFGS gets ~90% of Newton’s convergence speed with ~10% of the cost.

When Full Newton wins: Very small problems (< 10 variables) where Hessian computation is cheap and you need extreme precision.

─────────────────────────────────────────────────

📝 Development Note

The hardest challenge was making numerical differentiation robust across all Real types (Float, Double, Decimal).

Problem: The step size h for (f(x+h) - f(x-h)) / 2h must be:

• Small enough to approximate the true derivative
• Large enough to avoid catastrophic cancellation (subtracting nearly equal floating-point numbers)

Solution: Adaptive step size h = √ε × max(|x|, 1) where ε is machine epsilon:

• Float (ε ≈ 10⁻⁷): h ≈ 10⁻³
• Double (ε ≈ 10⁻¹⁶): h ≈ 10⁻⁸
• Decimal (ε ≈ 10⁻²⁸): h ≈ 10⁻¹⁴

This automatically adjusts to the numeric type’s precision.

Related Methodology: Numerical Stability (Week 2) - Covered machine epsilon and catastrophic cancellation.

Next Steps

Coming up Wednesday: Constrained Optimization - Lagrange multipliers, augmented Lagrangian methods, and handling equality/inequality constraints.

Series Progress:

• Week: 8/12
• Posts Published: 26/~48
• Playgrounds: 22 available