Constrained Optimization: Lagrange Multipliers and Real-World Constraints

BusinessMath Quarterly Series

13 min read

Part 27 of 12-Week BusinessMath Series


What You’ll Learn


The Problem

Real-world optimization has constraints:

Unconstrained optimization finds solutions that violate real-world constraints. Post-hoc normalization (e.g., dividing by sum) doesn’t minimize the original objective.


The Solution

BusinessMath provides constrained optimization via augmented Lagrangian methods. Constraints are first-class citizens, satisfied throughout optimization—not normalized after the fact.

Type-Safe Constraint Infrastructure

The MultivariateConstraint enum provides type-safe constraint specification:

import BusinessMath

// Equality constraint: x + y = 1
let equality: MultivariateConstraint
            
              > = .equality { v in v[0] + v[1] - 1.0 } // Inequality constraint: x ≥ 0 → -x ≤ 0 let inequality: MultivariateConstraint
              
                > = .inequality { v in -v[0] } // Check if satisfied let point = VectorN([0.5, 0.5]) print("Equality satisfied: \(equality.isSatisfied(at: point))") // true print("Inequality satisfied: \(inequality.isSatisfied(at: point))") // true 
              
            

Pre-Built Constraint Helpers

BusinessMath provides common constraint patterns:

import BusinessMath

// Budget constraint: weights sum to 1
let budget = MultivariateConstraint
            
              >.budgetConstraint // Non-negativity: all components ≥ 0 (long-only) let longOnly = MultivariateConstraint
              
                >.nonNegativity(dimension: 5) // Position limits: each weight ≤ 30% let positionLimits = MultivariateConstraint
                
                  >.positionLimit(0.30, dimension: 5) // Box constraints: 5% ≤ wᵢ ≤ 40% let box = MultivariateConstraint
                  
                    >.boxConstraints( min: 0.05, max: 0.40, dimension: 5 ) // Combine multiple constraints let allConstraints = [budget] + longOnly + positionLimits 
                  
                
              
            

Equality-Constrained Optimization

Problem: Minimize f(x) subject to h(x) = 0

Example: Minimize portfolio risk subject to weights summing to 100%

import BusinessMath

// Minimize x² + y² subject to x + y = 1
let objective: (VectorN
            
              ) -> Double = { v in v[0]*v[0] + v[1]*v[1] } let constraints = [ MultivariateConstraint
              
                >.equality { v in v[0] + v[1] - 1.0 // x + y = 1 } ] let optimizer = ConstrainedOptimizer
                
                  >() let result = try optimizer.minimize( objective, from: VectorN([0.0, 1.0]), subjectTo: constraints ) print("Solution: \(result.solution.toArray().map({ $0.number(4) }))") print("Objective: \(result.objectiveValue.number(6))") print("Constraint satisfied: \(constraints[0].isSatisfied(at: result.solution))") 
                
              
            

Output:

Solution: [0.5000, 0.5000]
Objective: 0.500000
Constraint satisfied: true

The insight: The optimal solution is where both variables equal 0.5, balancing the objective (minimize sum of squares) with the constraint (sum to 1).


Shadow Prices (Lagrange Multipliers)

The Lagrange multiplier λ tells you how much the objective improves if you relax the constraint by one unit.

let result = try optimizer.minimize(objective, from: initial, subjectTo: constraints)

if let multipliers = result.lagrangeMultipliers {
    for (i, λ) in multipliers.enumerated() {
        print("Constraint \(i): λ = \(λ.number(3))")
        print("  Marginal value of relaxing: \(λ.number(3)) per unit")
    }
}

Output:

Constraint 0: λ = -0.999
  Marginal value of relaxing: -0.999 per unit

Interpretation: If we relax “x + y = 1” to “x + y = 1.01”, the objective improves by ~0.005 (λ × 0.01).

Applications:


Inequality-Constrained Optimization

Problem: Minimize f(x) subject to g(x) ≤ 0

Example: Portfolio optimization with no short-selling and position limits

import BusinessMath
import Foundation

// Portfolio variance
let covariance = [
    [0.04, 0.01, 0.02],
    [0.01, 0.09, 0.03],
    [0.02, 0.03, 0.16]
]

let portfolioVariance: (VectorN
            
              ) -> Double = { w in var variance = 0.0 for i in 0..<3 { for j in 0..<3 { variance += w[i] * w[j] * covariance[i][j] } } return variance } // Constraints let constraints: [MultivariateConstraint
              
                >] = [ // Budget: weights sum to 1 .equality { w in w.reduce(0, +) - 1.0 }, // Long-only: wᵢ ≥ 0 → -wᵢ ≤ 0 .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, // Position limits: wᵢ ≤ 0.5 → wᵢ - 0.5 ≤ 0 .inequality { w in w[0] - 0.5 }, .inequality { w in w[1] - 0.5 }, .inequality { w in w[2] - 0.5 } ] let optimizer = InequalityOptimizer
                
                  >() let result = try optimizer.minimize( portfolioVariance, from: VectorN([1.0/3, 1.0/3, 1.0/3]), subjectTo: constraints ) print("Optimal weights: \(result.solution.toArray().map({ $0.percent(1) }))") print("Portfolio risk: \(sqrt(result.objectiveValue).percent(2))") print("All constraints satisfied: \(constraints.allSatisfy { $0.isSatisfied(at: result.solution) })") 
                
              
            

Output:

Optimal weights: ["50.0%", "36.8%", "13.2%"]
Portfolio risk: 18.50%
All constraints satisfied: true

The result: Asset 1 (lowest variance) gets the highest allocation, but capped at position limit. Constraint-aware optimization finds the true optimum.


Real-World Example: Target Return Portfolio

Minimize risk subject to achieving a target return:

import BusinessMath

let expectedReturns = VectorN([0.08, 0.10, 0.12, 0.15])
let covarianceMatrix = [
    [0.0400, 0.0100, 0.0080, 0.0050],
    [0.0100, 0.0625, 0.0150, 0.0100],
    [0.0080, 0.0150, 0.0900, 0.0200],
    [0.0050, 0.0100, 0.0200, 0.1600]
]

// Objective: Minimize variance
func portfolioVariance(_ weights: VectorN
            
              ) -> Double { var variance = 0.0 for i in 0..
              
                >() let result = try optimizer.minimize( portfolioVariance, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ // Fully invested .equality { w in w.reduce(0, +) - 1.0 }, // Target return ≥ 12% .inequality { w in let ret = w.dot(expectedReturns) return 0.12 - ret // ≤ 0 means ret ≥ 12% }, // Long-only .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] } ] ) print("Optimal weights: \(result.solution.toArray().map({ $0.percent(1) }))") let optimalReturn = result.solution.dot(expectedReturns) let optimalRisk = sqrt(portfolioVariance(result.solution)) print("Expected return: \(optimalReturn.percent(2))") print("Volatility: \(optimalRisk.percent(2))") print("Sharpe ratio (rf=3%): \((optimalReturn - 0.03) / optimalRisk)") 
              
            

Output:

Optimal weights: ["11.0%", "25.9%", "31.2%", "31.9%"]
Expected return: 12.00%
Volatility: 19.81%
Sharpe ratio (rf=3%): 0.4542157498481902

The solution: The optimizer found the minimum-risk portfolio that achieves exactly 12% return. Asset 4 (highest return but highest risk) gets only 31.9% because we’re minimizing risk, not maximizing return.


Comparing Constrained vs. Unconstrained

// Unconstrained: Minimize variance (allows short-selling, arbitrary weights)
let unconstrainedOptimizer = MultivariateNewtonRaphson
            
              >() let unconstrained = try unconstrainedOptimizer.minimizeBFGS( function: portfolioVariance, gradient: { try numericalGradient(portfolioVariance, at: $0) }, initialGuess: VectorN([0.25, 0.25, 0.25, 0.25]) ) print("Unconstrained solution: \(unconstrained.solution.toArray().map({ $0.percent(1) }))") print("Sum of weights: \((unconstrained.solution.reduce(0, +)).percent(1))") print("\n=== Impact of Constraints ===\n") let constrainedOptimizer = InequalityOptimizer
              
                >() // Budget-only: Minimum variance with just the budget constraint (allows shorting) let budgetOnly = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 } // Only budget constraint ] ) print("Budget-only (allows shorting):") print(" Weights: \(budgetOnly.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(budgetOnly.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(budgetOnly.solution)).percent(2))") // Long-only: Add non-negativity constraints let longOnly_option = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 }, .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] } ] ) print("\nLong-only (no short positions):") print(" Weights: \(longOnly_option.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(longOnly_option.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(longOnly_option.solution)).percent(2))") // Position limits: Add 40% maximum per position let positionLimited = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 }, .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] }, .inequality { w in w[0] - 0.40 }, .inequality { w in w[1] - 0.40 }, .inequality { w in w[2] - 0.40 }, .inequality { w in w[3] - 0.40 } ] ) print("\nPosition-limited (max 40% per asset):") print(" Weights: \(positionLimited.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(positionLimited.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(positionLimited.solution)).percent(2))") print("\n💡 Note: More constraints → higher variance (constraints limit optimization)") print(" But constraints reflect real-world limitations (no shorting, diversification rules, etc.)") 
              
            

Output:

Unconstrained solution: [150.2%, -25.3%, -18.7%, -6.2%]
Sum of weights: 100.0%

Constrained solution: [62.5%, 25.3%, 12.2%, 0.0%]
Sum of weights: 100.0%

The difference: Unconstrained allows short-selling (negative weights), which may be unrealistic. Constrained enforces real-world requirements.


Try It Yourself

Click to expand full playground code
import BusinessMath
import Foundation

// MARK: - Basic Constraint Infrastructure

// Equality constraint: x + y = 1
let equality: MultivariateConstraint
              
                > = .equality { v in let x = v[0], y = v[1] return x + y - 1.0 } // Inequality constraint: x ≥ 0 → -x ≤ 0 let inequality: MultivariateConstraint
                
                  > = .inequality { v in -v[0] } // Check if satisfied let point = VectorN([0.5, 0.5]) print("Equality satisfied: \(equality.isSatisfied(at: point))") // true print("Inequality satisfied: \(inequality.isSatisfied(at: point))") // true // MARK: - Pre-Built Helpers // Budget constraint: weights sum to 1 let budget = MultivariateConstraint
                  
                    >.budgetConstraint // Non-negativity: all components ≥ 0 (long-only) let longOnly = MultivariateConstraint
                    
                      >.nonNegativity(dimension: 5) // Position limits: each weight ≤ 30% let positionLimits = MultivariateConstraint
                      
                        >.positionLimit(0.30, dimension: 5) // Box constraints: 5% ≤ wᵢ ≤ 40% let box = MultivariateConstraint
                        
                          >.boxConstraints( min: 0.05, max: 0.40, dimension: 5 ) // Combine multiple constraints let allConstraints = [budget] + longOnly + positionLimits // MARK: - Equality-Constrained Optimization // Minimize x² + y² subject to x + y = 1 let objective_eqConst: (VectorN
                          
                            ) -> Double = { v in let x = v[0], y = v[1] return x*x + y*y } let constraints_eqConst = [ MultivariateConstraint
                            
                              >.equality { v in v[0] + v[1] - 1.0 // x + y = 1 } ] let optimizer_eqConst = ConstrainedOptimizer
                              
                                >() let result_eqConst = try optimizer_eqConst.minimize( objective_eqConst, from: VectorN([0.0, 1.0]), subjectTo: constraints_eqConst ) print("Solution: \(result_eqConst.solution.toArray().map({ $0.number(4) }))") print("Objective: \(result_eqConst.objectiveValue.number(6))") print("Constraint satisfied: \(constraints_eqConst[0].isSatisfied(at: result_eqConst.solution))") for (i, λ) in result_eqConst.lagrangeMultipliers.enumerated() { print("Constraint \(i): λ = \(λ.number(3))") print(" Marginal value of relaxing: \(λ.number(3)) per unit") } // MARK: Inequality-Constrained Example // Portfolio variance let covariance_portfolio = [ [0.04, 0.01, 0.02], [0.01, 0.09, 0.03], [0.02, 0.03, 0.16] ] let portfolioVariance_portfolio: (VectorN
                                
                                  ) -> Double = { w in var variance = 0.0 for i in 0..<3 { for j in 0..<3 { variance += w[i] * w[j] * covariance_portfolio[i][j] } } return variance } // Constraints let constraints_portfolio: [MultivariateConstraint
                                  
                                    >] = [ // Budget: weights sum to 1 .equality { w in w.reduce(0, +) - 1.0 }, // Long-only: wᵢ ≥ 0 → -wᵢ ≤ 0 .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, // Position limits: wᵢ ≤ 0.5 → wᵢ - 0.5 ≤ 0 .inequality { w in w[0] - 0.5 }, .inequality { w in w[1] - 0.5 }, .inequality { w in w[2] - 0.5 } ] let optimizer_portfolio = InequalityOptimizer
                                    
                                      >() let result_portfolio = try optimizer_portfolio.minimize( portfolioVariance_portfolio, from: VectorN([1.0/3, 1.0/3, 1.0/3]), subjectTo: constraints_portfolio ) print("Optimal weights: \(result_portfolio.solution.toArray().map({ $0.percent(1) }))") print("Portfolio risk: \(sqrt(result_portfolio.objectiveValue).percent(2))") print("All constraints satisfied: \(constraints_portfolio.allSatisfy { $0.isSatisfied(at: result_portfolio.solution) })") // MARK: - Target Return Portfolio let expectedReturns_targetP = VectorN([0.08, 0.10, 0.12, 0.15]) let covarianceMatrix_targetP = [ [0.0400, 0.0100, 0.0080, 0.0050], [0.0100, 0.0625, 0.0150, 0.0100], [0.0080, 0.0150, 0.0900, 0.0200], [0.0050, 0.0100, 0.0200, 0.1600] ] // Objective: Minimize variance func portfolioVariance_targetP(_ weights: VectorN
                                      
                                        ) -> Double { var variance = 0.0 for i in 0..
                                        
                                          >() let result_targetP = try optimizer_targetP.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ // Fully invested .equality { w in w.reduce(0, +) - 1.0 }, // Target return ≥ 12% .inequality { w in let ret = w.dot(expectedReturns_targetP) return 0.12 - ret // ≤ 0 means ret ≥ 12% }, // Long-only .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] } ] ) print("Optimal weights: \(result_targetP.solution.toArray().map({ $0.percent(1) }))") let optimalReturn_targetP = result_targetP.solution.dot(expectedReturns_targetP) let optimalRisk_targetP = sqrt(portfolioVariance_targetP(result_targetP.solution)) print("Expected return: \(optimalReturn_targetP.percent(2))") print("Volatility: \(optimalRisk_targetP.percent(2))") print("Sharpe ratio (rf=3%): \((optimalReturn_targetP - 0.03) / optimalRisk_targetP)") // MARK: - Comparing Constrained vs Fewer Constraints print("\n=== Impact of Constraints ===\n") let constrainedOptimizer = InequalityOptimizer
                                          
                                            >() // Budget-only: Minimum variance with just the budget constraint (allows shorting) let budgetOnly = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 } // Only budget constraint ] ) print("Budget-only (allows shorting):") print(" Weights: \(budgetOnly.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(budgetOnly.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(budgetOnly.solution)).percent(2))") // Long-only: Add non-negativity constraints let longOnly_option = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 }, .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] } ] ) print("\nLong-only (no short positions):") print(" Weights: \(longOnly_option.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(longOnly_option.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(longOnly_option.solution)).percent(2))") // Position limits: Add 40% maximum per position let positionLimited = try constrainedOptimizer.minimize( portfolioVariance_targetP, from: VectorN([0.25, 0.25, 0.25, 0.25]), subjectTo: [ .equality { w in w.reduce(0, +) - 1.0 }, .inequality { w in -w[0] }, .inequality { w in -w[1] }, .inequality { w in -w[2] }, .inequality { w in -w[3] }, .inequality { w in w[0] - 0.40 }, .inequality { w in w[1] - 0.40 }, .inequality { w in w[2] - 0.40 }, .inequality { w in w[3] - 0.40 } ] ) print("\nPosition-limited (max 40% per asset):") print(" Weights: \(positionLimited.solution.toArray().map({ $0.percent(1) }))") print(" Variance: \(portfolioVariance_targetP(positionLimited.solution).number(6))") print(" Volatility: \(sqrt(portfolioVariance_targetP(positionLimited.solution)).percent(2))") print("\n💡 Note: More constraints → higher variance (constraints limit optimization)") print(" But constraints reflect real-world limitations (no shorting, diversification rules, etc.)") 
                                          
                                        
                                      
                                    
                                  
                                
                              
                            
                          
                        
                      
                    
                  
                
              

→ Full API Reference: BusinessMath Docs – 5.6 Constrained Optimization

Modifications to try:

  1. Add sector constraints (e.g., max 40% in any sector across multiple assets)
  2. Optimize a production mix with material, labor, and capacity constraints
  3. Build a portfolio with leverage constraints (130/30 strategy)
  4. Compare shadow prices: which constraints are binding?

Real-World Application

Portfolio manager use case: “I need to build a portfolio with:

Constrained optimization solves this exactly—no manual tweaking required.


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

Why Post-Hoc Normalization Doesn’t Work

Common mistake:

// ❌ Wrong: Normalize after unconstrained optimization
let weights = unconstrainedOptimizer.minimize(variance)
let normalized = weights / weights.sum()  // Not optimal!

Why it’s wrong:

  1. The unconstrained optimum is at a different point in parameter space
  2. Normalizing changes the objective value (variance ≠ variance after scaling)
  3. Violates constraint throughout optimization (no feedback to guide search)

Correct approach:

// ✅ Right: Constraints during optimization
let result = optimizer.minimize(
    variance,
    subjectTo: [.budgetConstraint, .longOnly]
)
// Constraint satisfied at every iteration

Rule: Constraints must be part of the optimization, not post-processing.

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


📝 Development Note

The hardest challenge was choosing the right constrained optimization algorithm. We evaluated:

  1. Penalty methods: Add constraint violations to objective
    • Simple but requires tuning penalty weights
    • Can be numerically unstable
  2. Augmented Lagrangian: Penalty + Lagrange multipliers
    • More robust than pure penalty
    • Self-adjusting penalties
    • What we chose
  3. Sequential Quadratic Programming (SQP): Second-order method
    • Fastest convergence
    • Complex implementation, requires Hessian

We chose Augmented Lagrangian because:

Related Methodology: Algorithm Selection (Week 1) - Covered how we evaluate trade-offs between implementation complexity and user benefit.


Next Steps

Coming up Friday: Case Study #4 - Real-world portfolio optimization combining everything from Weeks 7-8 (goal-seeking, multivariate optimization, constraints, and risk models).


Series Progress:


Tagged with: optimization