We continue our study of methods to solve systems of equations, now looking at systems that contain nonlinear equations
Unlike linear systems, with nonlinear systems we don’t have guarantees of closed-form solutions, so solving them is harder
We will study two types of methods: derivative-free and derivative-based
0.3 Main references for today
Miranda & Fackler (2002), Ch. 3
Judd (1998), Ch. 5
Nocedal & Wright (2006), Ch. 11
Lecture notes for Ivan Rudik’s Dynamic Optimization (Cornell)
1 Types of nonlinear equation problems
1.1 Types of nonlinear equation problems
Problems involving nonlinear equations are very common in economics
We can classify them into two types
Systems of nonlinear equations
Complementarity problems
1.2 System of nonlinear equations
These come in two forms
Rootfinding problems
For a function f: \mathbb{R}^n \rightarrow \mathbb{R}^n, we want to find a n-vector x that satisfies f(x) = 0
\rightarrow Any x that satisfies this condition is called a root of f
Examples?
. . .
Market equilibrium in general: market clearing conditions
No-arbitrage conditions (pricing models)
Solving first-order optimality conditions
1.3 System of nonlinear equations
Fixed-point problems
For a function g: \mathbb{R}^n \rightarrow \mathbb{R}^n, we want to find a n-vector x that satisfies x = g(x)\rightarrow Any x that satisfies this condition is called a fixed point of g
Examples?
. . .
Best-response functions
Many equilibrium concepts in game theory
A Nash equilibrium is a fixed point in the strategy space
1.4 System of nonlinear equations
Rootfinding and fixed-point problems are equivalent
We can easily convert one into another
Rootfinding\rightarrowfixed-point
. . .
Define a new g(x) = x - f(x)
. . .
Fixed-point\rightarrowrootfinding
Define a new f(x) = x - g(x)
1.5 Complementarity problems
In these problems, we have
A function f: \mathbb{R}^n \rightarrow \mathbb{R}^n
n-vectors a and b, with a<b
And we are looking for an n-vector x \in [a,b] such that for all i = 1,\dots,n
If the constraints on x_i do not bind ( a_i < x_i < b_i) , then the first-order condition is precisely zero
. . .
But suppose the upper bound binds ( x_i = b_i ). Then f_i(x) \geq 0 since x_i > a_i
But we can’t guarantee that f_i(x) = 0 because f_i(x) might still be increasing at that point
1.7 Complementarity problems
Rootfinding is a special case of complementarity problems: a = -\infty and b = \infty
But complementarity problems are not just about finding a root within[a, b]
Remember: if some x_i is at the boundary ( x_i = a_i or x_i = b_i ), some element of f(x) can be non-zero!
. . .
NOTE: We will complementary problems next lecture (separate slide deck)
2 Derivative-free methods
2.1 Rootfinding problems
Let’s start simple: we have a continuous function f:[a,b]\in\mathbb{R} \rightarrow\mathbb{R} and we know that f(a)<0 and f(b)>0
What does the Intermediate Value Theorem says here?
. . .
If f is continuous and f(a) \neq f(b), then f must assume all values in between f(a) and f(b)
So if f(a)<0 and f(b)>0 \Rightarrow there must be at least one root x\in[a,b] such that f(x) = 0
. . .
How would you go about finding a root?
2.2 Bisection method
Basic idea: split the search interval in two parts and check whether there’s a root in each part
How do we check that?
. . .
By looking at the signs of f(x) at the boundaries of each interval
If they are different, there’s a root \Rightarrow we keep looking there
Let’s see an illustration
2.3 Bisection method
We start with l = a, u = b and find m = (u+l)/2
Let’s say f(l)<0, f(u)>0, and f(m)>0. What do we do next?
2.4 Bisection method
Since f(l)<0 and f(m)>0 have different signs, we continue our search in [l,m]
We set u \leftarrow m
. . .
Then we calculate the new midpoint m
. . .
Now say f(l)<0, f(u)>0, and f(m)<0. What do we do next?
2.5 Bisection method
Since f(m)<0 and f(u)>0 have different signs, we continue our search in [m,u]
We set l \leftarrow m
And the search continues until we are satisfied with the precision
2.6 Try it! Bisection method
Your turn!
Write a function that takes (f, a, b) and returns a root of f using the bisection method.
Then, use it to find the root of f(x) = -x^{-2} + x - 1 between [0.2, 4]
Here are the basic steps:
Start with a lower ( l = a ) and an upper ( u = b ) bounds
Get the midpoint m = (u+l)/2
Check the sign of f(m)
If sign(f(m)) = sign(f(l)), move lower bound up: l \leftarrow m
If sign(f(m)) = sign(f(u)), move upper bound down: u \leftarrow m
Repeat 2 and 3 until our interval is short enough ( (u-l)/2 < tol ) and return x = m
2.7 Bisection method
Show the code
functionbisection(f, lo, up) tolerance =1e-3# tolerance for solution mid = (lo + up)/2# initial guess, bisect the interval difference = (up - lo)/2# initialize bound differencewhile difference > tolerance # loop until convergenceprintln("Intermediate guess: $mid")ifsign(f(lo)) ==sign(f(mid)) # if the guess has the same sign as the lower bound lo = mid # a solution is in the upper half of the interval mid = (lo + up)/2else# else the solution is in the lower half of the interval up = mid mid = (lo + up)/2end difference = (up - lo)/2# update the difference endprintln("The root of f(x) is $mid")end;
So it’s a good idea to check if you have the right boundaries
Code
f(2.0)
0.75
Code
f(4.0)
2.9375
These are both positive. So by the IVT, we can’t know for sure if there’s a root here
2.12 Bisection method
The bisection method is incredibly robust: if a function f satisfies the IVT, it is guaranteed to converge in a finite number of iterations
. . .
A root can be calculated to arbitrary precision \epsilon in a maximum of log_2\frac{b-a}{\epsilon} iterations
But robustness comes with drawbacks:
It only works in one dimension
It is slow because it only uses information about the function’s level but not its variation
2.13 Function iteration
For the next method, we recast rootfinding as a fixed-point problem
. . .
f(x) = 0 \Rightarrow g(x) = x - f(x)
Then, we start with an initial guess x^{(0)} and iterate x^{(k+1)}\leftarrow g(x^{k}) until convergence: |x^{(k+1)}-x^{(k)}| \approx 0
. . .
Let’s try it!
Write a function that takes (f, initial_guess) and returns a root of f using function iteration
Then, use it to find the root of f(x) = -x^{-2} + x - 1 with initial_guess = 1
2.14 Function iteration
Show the code
functionfunction_iteration(f, initial_guess) tolerance =1e-3# tolerance for solution difference =Inf# initialize difference x = initial_guess # initialize current valuewhileabs(difference) > tolerance # loop until convergenceprintln("Intermediate guess: $x") x_prev = x # store previous value x = x_prev -f(x_prev) # calculate next guess difference = x - x_prev # update differenceendprintln("The root of f(x) is $x")end;
A fixed point x = g(x) is at the intersection between g(x) and the 45o line
Starting from x^{(0)}, we calculate g(x^{(0)}) and find the corresponding point on the 45o line for x^{(1)}
We keep iterating until we (approximately) find the fixed point
2.17 Function iteration
Function iteration is guaranteed to converge to a fixed point x^* if
g is differentiable, and
the initial guess is “sufficiently close” to an x^* at which \|g^\prime (x^{*}) \| < 1
It may also converge when these conditions are not met
Since this is an easy method to implement, it’s worth trying it before switching to more complex methods
2.18 Function iteration
But wait: What is “sufficiently close”?
. . .
Good question! There is no practical formula. As Miranda and Fackler (2002) put it
Typically, an analyst makes a reasonable guess for the root f and counts his/her blessings if the iterates converge. If the iterates do not converge, then the analyst must look more closely at the properties of f to find a better starting value, or change to another rootfinding method.
This is where science also becomes a bit of an art
3 Derivative-based methods
3.1 Newton’s method
Newton’s method and variants are the workhorses of solving n-dimensional non-linear problems
. . .
Key idea: take a hard non-linear problem and replace it with a sequence of linear problems \rightarrowsuccessive linearization
3.2 Newton’s method
Start with an initial guess of the root at x^{(0)}
Approximate the non-linear function with its first-order Taylor expansion around x^{(0)}
This is just the tangent line at x^0
Solve for the root of this linear approximation, call it x^{(1)}
3.3 Newton’s method
Repeat starting at x^{(1)} until we converge to x^*
This can be applied to a function with an arbitrary number of dimensions
3.4 Newton’s method
Formally: begin with some initial guess of the root vector \mathbf{x^{(0)}}
. . .
Given \mathbf{x^{(k)}}, our new guess \mathbf{x^{(k+1)}} is obtained by approximating f(\mathbf{x}) using a first-order Taylor expansion around \mathbf{x^{(k)}}
functionnewtons_method(f, f_prime, initial_guess) tolerance =1e-3# tolerance for solution difference =Inf# initialize difference x = initial_guess # initialize current valuewhileabs(difference) > tolerance # loop until convergenceprintln("Intermediate guess: $x") x_prev = x # store previous value x = x_prev -f(x_prev)/f_prime(x_prev) # calculate next guess# ^ this is the only line that changes from function iteration difference = x - x_prev # update differenceendprintln("The root of f(x) is $x")end;
3.7 Newton’s method
Code
f(x) =-x^(-2) + x -1;f_prime(x) =2x^(-3) +1;newtons_method(f, f_prime, 1.0)
Intermediate guess: 1.0
Intermediate guess: 1.3333333333333333
Intermediate guess: 1.4576271186440677
Intermediate guess: 1.4655459327062879
The root of f(x) is 1.465571231622256
3.8 Newton’s method
Newton’s method has nice properties regarding convergence and speed
It converges if
If f(x) is continuously differentiable,
The initial guess is “sufficiently close” to the root, and
f(x) is invertible near the root
. . .
We need f(x) to be invertible so the algorithm above is well defined
. . .
If f'(x) is ill-conditioned we can run into problems with rounding error
3.9 Newton’s method: a duopoly example
Inverse demand: P(q) = q^{-1/\epsilon}
Two firms with costs: C_i(q_i) = \frac{1}{2}c_i q_i^2
Visually, this is the path Newton’s method followed
3.21 Newton’s method: analytic vs numerical derivatives
It was tedious but in the previous example we could calculate the Jacobian analytically. Sometimes it’s much harder and we can make mistakes.
Here is a challenge for you: redo the previous example but, instead of defining the Jacobian analytically, use the ForwardDiff package to do the derivatives for you. Then, compare your solution with mine.
. . .
An alternative is to use the Symbolics.jl package. For low-dimensional derivative and Jacobians, it does a good job!
3.22 Quasi-Newton methods
We usually don’t want to deal with analytic derivatives unless we have access to autodifferentiation
Why?
. . .
It can be difficult to do the analytic derivation
Coding a complicate Jacobian is prone to errors and takes time
Can actually be slower to evaluate than finite differences for a nonlinear problem
. . .
Alternative \rightarrowfinite differences instead of analytic derivatives
3.23 Quasi-Newton: Secant method
Using our current root guess x^{(k)} and our previous root guess x^{(k-1)}:
We must initially provide a guess of the root, x^{(0)}, but also a guess of the Jacobian, A_{(0)}
A good guess for A_{(0)} is to calculate it numerically at our chosen x^{(0)}
3.26 Quasi-Newton: Broyden’s method
The iteration rule is the same as before but with our guess of the Jacobian substituted in for the actual Jacobian (or the finite difference approximation)
Any reasonable guess for the Jacobian should satisfy this condition
. . .
But this gives n conditions with n^2 elements to solve for in A
3.27 Quasi-Newton: Broyden’s method
Broyden’s methods solves this under-determined problem with an assumption that focuses on the direction we are most interested in: d^{(k)} = \mathbf{x^{(k+1)}} - \mathbf{x^{(k)}}
For any direction q orthogonal to d^{(k)}, it assumes that A^{(k+1)}q = A^{(k)}q
In other words, our next guess is as good as the current ones for any changes in x that are orthogonal to the one we are interested right now
. . .
Jointly, the secant condition and the orthogonality assumption give the iteration rule for the Jacobian:
functionf(q) Q =sum(q) F = Q^(-1/epsilon) .- (1/epsilon)Q^(-1/epsilon-1) .*q .- c .*qend;
Do you see any potential inefficiency?
. . .
We allocate a new F every time we call this function!
Instead, we can be more efficient by writing a function that modifies a pre-allocated vector
3.38 Quick detour: functions that modify arguments
By convention, in Julia we name functions that modify arguments with a ! at the end. For our f, we can define
Code
functionf!(F, q) F[1] =sum(q)^(-1/epsilon) - (1/epsilon)sum(q)^(-1/epsilon-1)*q[1] - c[1]*q[1] F[2] =sum(q)^(-1/epsilon) - (1/epsilon)sum(q)^(-1/epsilon-1)*q[2] - c[2]*q[2]end;F =zeros(2) # This allocates a 2-vector with elements equal to zerof!(F, [0.2; 0.2]); # Note the first argument is a pre-allocated vectorF
Rootfinding algorithms will converge at different speeds in terms of the number of operations
. . .
A sequence of iterates x^{(k)} is said to converge to x^* at a rate of order p if there is a constant C such that
\|x^{(k+1)}-x^{*}\|\leq C \|x^{(k)} - x^{*}\|^p
for sufficiently large k
4.2 Convergence speed
\|x^{(k+1)}-x^{*}\|\leq C \|x^{(k)} - x^{*}\|^p
If C < 1 and p = 1: linear convergence
If 1 < p < 2: superlinear convergence
If p = 2: quadratic convergence
The higher order the convergence rate, the faster it converges
4.3 Convergence speed
How fast do the methods we’ve seen converge?
. . .
Bisection: linear rate with C = 0.5
. . .
Function iteration: linear rate with C = ||f'(x^*)||
. . .
Secant and Broyden: superlinear rate with p \approx 1.62
. . .
Newton: p = 2
4.4 Convergence speed
Consider an example where f(x) = x - \sqrt(x) = 0
This is how the 3 main approaches converge in terms of the L^1-norm for an initial guess x^{(0)} = 0.5
4.5 Choosing a solution method
Convergence rates only account for the number of iterations of the method
The steps taken in a given iteration of each solution method may vary in computational cost because of differences in the number of arithmetic operations
Although an algorithm may take more iterations to solve, each iteration may be solved faster and the overall algorithm takes less time
4.6 Choosing a solution method
Bisection method only requires a single function evaluation during each iteration
Function iteration only requires a single function evaluation during each iteration
Broyden’s method requires both a function evaluation and matrix multiplication
Newton’s method requires a function evaluation, a derivative evaluation, and solving a linear system
. . .
\rightarrow Bisection and function iteration are usually slow
\rightarrow Broyden’s method can be faster than Newton’s method if derivatives are costly to compute
4.7 Choosing a solution method
Besides convergence rates and algorithm speed, you should also factor development time
Newton’s method is fastest to converge
If deriving and programming the Jacobian is relatively simple and not too costly to compute, this method is a good choice
If derivatives are complex, quasi-Newton methods are good candidates
Bisection and function iteration are generally dominated options but are easy to program/debug, so they have value as a quick proof-of-concept
Bisection is often used in hybrid methods, such as Dekker’s and Brent’s. Hybrid methods select between bisection, secant, or other basic solution methods every iteration depending on a set of criteria
5 Common problems and solutions
5.1 When solvers fail: domain errors
A very common problem with Newton and Broyden: the solver tries values outside the function’s domain
. . .
Consider a market equilibrium with isoelastic supply S_i = c_i p_i^{\sigma_i} and demand D_i = y\, p_i^{\epsilon_{ii}} p_j^{\epsilon_{ij}}
# From an asymmetric guess, it crashes:try NLsolve.nlsolve(f!, [0.05; 2.0], method=:broyden)catcheprintln("Error: $e")end
Error: DomainError(-3.57491732087526, "Exponentiation yielding a complex result requires a complex argument.\nReplace x^y with (x+0im)^y, Complex(x)^y, or similar.")
. . .
What happened? Broyden’s approximate Jacobian led to a step that pushed p_1 below zero. Julia tried to compute a negative number raised to a fractional power and threw a DomainError
5.4 Why this happens
Newton and Broyden are unconstrained methods: they search over all of \mathbb{R}^n
. . .
This is problematic when:
Variables must be positive (prices, quantities, income)
You use isoelastic functions: x^{\alpha} with non-integer \alpha requires x > 0
Your initial guess is far from the solution \rightarrow large steps \rightarrow overshooting
. . .
Broyden is more prone to this than Newton because its Jacobian approximation can degrade, leading to poorly directed steps
. . .
Three common strategies to deal with this:
Start with a better initial guess
Try a different method
Transform variables so the solver’s domain is unrestricted
5.5 Strategy 1: Better initial guesses
Idea: Solve a simpler version of the problem first, then use that solution as the starting point
. . .
For a system with cross-price effects, we can ignore them and solve each market in isolation
This has a closed-form solution and gives a reasonable starting point for the coupled system
5.6 Strategy 1: Better initial guesses
Code
# Solve each market in isolation (set cross-price elasticities to zero)p0 = [(y₀/c[i])^(1/(σ[i] - ϵ_own[i])) for i in1:2];println("Simplified guess: $p0");
This is the same system, parameterized differently
Same equilibrium: if z^* solves g = 0, then p^* = e^{z^*} solves f = 0
But now the solver can explore any z \in \mathbb{R}^2 without ever triggering a domain error
5.10 Strategy 3: Variable transformation
Code
functiong!(G, z) G[1] = c[1]*exp(σ[1]*z[1]) -y₀*exp(ϵ_own[1]*z[1] + ϵ_cross[1]*z[2]) G[2] = c[2]*exp(σ[2]*z[2]) -y₀*exp(ϵ_cross[2]*z[1] + ϵ_own[2]*z[2])end# Same starting point that failed before: p = (0.05, 2.0) → z = log.((0.05, 2.0))result = NLsolve.nlsolve(g!, log.([0.05; 2.0]), method=:broyden)p_solution =exp.(result.zero)println("Solution: p = $p_solution")
Solution: p = [1.877722541843727, 1.5818822321046384]
No DomainError. The solver explored negative z values, and e^z kept prices positive throughout
5.11 When to use each strategy
Better initial guess: simple, requires no code changes to f, but needs problem-specific insight
Different method: Newton with autodiff is more robust than Broyden with finite differences, at the cost of more computation per step
Log transformation: most robust, works from distant guesses, but adds a layer of indirection and can lead to scaling issues if the solution is very large/small
. . .
Note: A related transformation for variables in (0,1) (e.g., probabilities or shares):
z = \log\!\left(\frac{x}{1-x}\right), \quad x = \frac{1}{1+e^{-z}}
This is the logit transformation: z \in \mathbb{R} maps to x \in (0,1)
5.12 When solvers fail: convergence problems
We’ve seen that Newton’s method converges if the initial guess is sufficiently close to a root where f' is invertible
But what does failure look like in practice?
Three classic failure modes:
Cycling: iterates bounce back and forth without approaching the root
Divergence: iterates move farther and farther from the root
Slow convergence: iterates approach the root, but VERY slowly
The method enters a perfect 2-cycle: the tangent at x = 0 hits zero at x = 1, and the tangent at x = 1 hits zero at x = 0
The actual root is at x^* \approx -1.77, which Newton never reaches from this starting point
5.14 Failure mode 1: Cycling
Let’s verify with a modified version of our newtons_method function
Code
f(x) = x^3-2*x +2;f_prime(x) =3*x^2-2;
We need to add a maximum iteration count so the loop terminates
Code
functionnewtons_method_maxiter(f, f_prime, initial_guess; maxiter=20) tolerance =1e-3# tolerance for solution difference =Inf# initialize difference x = initial_guess # initialize current valuefor k in1:maxiterprintln("Iteration $k: x = $(round(x, digits=6))") x_prev = x x = x_prev -f(x_prev)/f_prime(x_prev) difference = x - x_prevabs(difference) < tolerance &&return xendprintln("Did not converge after $maxiter iterations")return xend;
5.15 Failure mode 1: Cycling
The iterates alternate between 0.0 and 1.0 indefinitely
Iteration 1: x = 0.0
Iteration 2: x = 1.0
Iteration 3: x = 0.0
Iteration 4: x = 1.0
Iteration 5: x = 0.0
Iteration 6: x = 1.0
Iteration 7: x = 0.0
Iteration 8: x = 1.0
Iteration 9: x = 0.0
Iteration 10: x = 1.0
Did not converge after 10 iterations
5.16 Failure mode 1: Cycling
Show the code
usingPlots; gr()f(x) = x^3-2*x +2;f_prime(x) =3*x^2-2;xgrid =range(-2.5, 2.0, length=300)p =plot(xgrid, f.(xgrid), label="f(x) = x³ - 2x + 2", color=:black, linewidth=2)hline!([0], color=:gray, linestyle=:dash, label=nothing)# Newton steps: tangent from (0, f(0)) to (1, 0), then (1, f(1)) to (0, 0)xs_cycle = [0.0, 1.0, 0.0, 1.0, 0.0]colors_cycle = [:royalblue, :crimson]for k in1:4 xk = xs_cycle[k]; xk1 = xs_cycle[k+1]plot!([xk, xk1], [f(xk), 0.0], color=colors_cycle[mod1(k,2)], linestyle=:dot, linewidth=1.5, label=nothing)plot!([xk1, xk1], [0.0, f(xk1)], color=colors_cycle[mod1(k,2)], linestyle=:dot, linewidth=1.5, label=nothing)endscatter!([0.0, 1.0], [f(0.0), f(1.0)], color=:red, markersize=6, label="Iterates")title!("Cycling: Newton bounces between x = 0 and x = 1")xlabel!("x"); ylabel!("f(x)")p
5.17 Failure mode 2: Divergence
Consider f(x) = \text{sign}(x)|x|^{1/3} (the cube root)
The root is at x^* = 0, with f'(x) = \frac{1}{3}|x|^{-2/3}. The Newton update simplifies to:
Iteration 1: x = 0.5
Iteration 2: x = -1.0
Iteration 3: x = 2.0
Iteration 4: x = -4.0
Iteration 5: x = 8.0
Iteration 6: x = -16.0
Iteration 7: x = 32.0
Iteration 8: x = -64.0
Did not converge after 8 iterations
Iteration 1: x = 1.0
Iteration 2: x = 0.666667
Iteration 3: x = 0.444444
Iteration 4: x = 0.296296
Iteration 5: x = 0.197531
Iteration 6: x = 0.131687
Iteration 7: x = 0.087791
Iteration 8: x = 0.058528
Iteration 9: x = 0.039018
Iteration 10: x = 0.026012
Iteration 11: x = 0.017342
Iteration 12: x = 0.011561
Iteration 13: x = 0.007707
Iteration 14: x = 0.005138
Iteration 15: x = 0.003425
Iteration 16: x = 0.002284
5.22 Failure mode 3: Slow convergence
Compare with a well-behaved case: f(x) = x^3 - 1 has a simple root at x^* = 1
Show the code
# Collect errors for x^3 (repeated root)f_slow(x) = x^3; fp_slow(x) =3*x^2;x =1.0; errors_slow = [abs(x)]for k in1:50 x = x -f_slow(x)/fp_slow(x)push!(errors_slow, abs(x))end# Collect errors for x^3 - 1 (simple root)f_fast(x) = x^3-1; fp_fast(x) =3*x^2;x =2.0; errors_fast = [abs(x -1.0)]for k in1:15 x = x -f_fast(x)/fp_fast(x)push!(errors_fast, abs(x -1.0))endp =plot(0:length(errors_slow)-1, errors_slow, marker=:circle, markersize=3, color=:crimson, label="x³ (repeated root, linear rate)", yscale=:log10, xlabel="Iteration", ylabel="Error (log scale)", ylims=(1e-16, 10), linewidth=1.5)plot!(0:length(errors_fast)-1, errors_fast, marker=:diamond, markersize=3, color=:royalblue, label="x³ - 1 (simple root, quadratic rate)", linewidth=1.5)title!("Convergence comparison: repeated vs simple root")p
5.23 Failure mode 3: Slow convergence
The simple root reaches machine precision in about 7 iterations
The repeated root (level and 1st derivative) is still crawling after 50 iterations. Why?
