# Non-linear root-finding: Works in Matlab, but not in Julia

**URL:** <https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858>\
**Category:** General Usage\
**Tags:** roots, nlsolve\
**Created:** [September 12, 2018, 1:37pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858 "2018-09-12T13:37:58Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![IljaK91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iljak91/32/44301_2.png) [@IljaK91](https://discourse.julialang.org/u/IljaK91)\
**Post date:** [September 12, 2018, 1:37pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/1 "2018-09-12T13:37:58Z")

</div>

Hey everyone,

I have a small numerical problem here. The main idea is that I try to approximate a univariate function, that needs to satisfy some constraint, with a polynomial. I had this already working with an older version of [NLsolve](https://github.com/JuliaNLSolvers/NLsolve.jl) and I run this code under Julia 0.6.4 using the package [CompEcon](https://github.com/JuliaPackageMirrors/CompEcon.jl). However, after letting the code rest for a while and updating packages, I stopped finding a solution to my problem. Here is the code:

```julia
using CompEcon
# Parameters
global γ = 2; # Risk aversion
global α = 0.4; # Capital share
global β = 0.9; # Discount factor

# First, create a grid using the QuantEcon toolbox
# Set the endpoints of approximation interval:
a = 2.0 # left endpoint
b = 3.0 # right endpoint
# Note: They have to be both integer or real. Important in Julia.

# Choose an approximation scheme. In this case, let us use an order 10
n = 3 # order of approximation
fspace = fundefn(:cheb, n, a, b) # define fspace
nodes = funnode(fspace)[1] # safe grid/nodes
Phi = funbase(fspace,nodes) # derive basis matrix Φ

using NLsolve

function f!(res, c_guess, fspace, nodes) # we want to find roots for f
  global γ,α,β

  K = nodes
  con = funeval(c_guess,fspace,K)[1]
  K_next = K.^α - con
  con_next = funeval(c_guess,fspace,K_next)[1]
  res = con.^(-γ) - β*con_next.^(-γ).*α.*K_next.^(α-1)
end

g!(res, c_guess) = f!(res, c_guess, fspace, nodes)
c_guess = Phi\nodes.^0.3 # some initial guess for parameters c
result = nlsolve(g!,c_guess)

```

I have put equivalent code in MATLAB, where it works and delivers a sensible solution. What may be the reason for this not to work in Julia and what can I do about it?

The Matlab code is the following (before running this code, you need to install the [CompEcon Toolbox](http://www4.ncsu.edu/~pfackler/compecon/toolbox.html)).

```julia

function [res] = ExampleProblem(c_guess,fspace,nodes)

gamma = 2;
alpha = 0.4;
beta = 0.9;

% Unscreened Investment

K = nodes;
con = funeval(c_guess,fspace,K);
K_next = K.^alpha - con;
con_next = funeval(c_guess,fspace,K_next);

% Save the residual
res = con.^(-gamma) - beta*con_next.^(-gamma).*alpha.*K_next.^(alpha-1);
end

```

```julia

n = 3;
a = 2;
b = 3;
fspace = fundefn('cheb', n, a, b);   
nodes = funnode(fspace);       
Phi = funbas(fspace,nodes);  

c_guess = Phi\nodes.^0.3;

options = optimset('display','iter');

y = fsolve(@(x)ExampleProblem(x,fspace,nodes),c_guess,options);

```

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [September 12, 2018, 2:00pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/2 "2018-09-12T14:00:05Z")

</div>

In your `f!`, the last line isn’t actually updating `res`. Try this:

```julia
res[:] = con.^(-γ) - β*con_next.^(-γ).*α.*K_next.^(α-1)

```

---

<div class="post-metadata">

**Author:** ![IljaK91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iljak91/32/44301_2.png) [@IljaK91](https://discourse.julialang.org/u/IljaK91)\
**Post date:** [September 12, 2018, 2:06pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/3 "2018-09-12T14:06:42Z")

</div>

Thanks! That helped. Can you give me an explanation why using `res[:]` is different? That would be very helpful 🙂

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [September 12, 2018, 2:08pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/4 "2018-09-12T14:08:18Z")

</div>

`res = ` creates a new variable. It has the same name as the function argument, so it looks confusing.

`res[:]` is a slicing operation that puts the result from the right-hand side into `res`.

---

<div class="post-metadata">

**Author:** ![IljaK91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iljak91/32/44301_2.png) [@IljaK91](https://discourse.julialang.org/u/IljaK91)\
**Post date:** [September 12, 2018, 2:10pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/5 "2018-09-12T14:10:52Z")

</div>

Makes sense, thanks a lot!

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [September 12, 2018, 2:49pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/6 "2018-09-12T14:49:15Z")

</div>

If the shapes match, in-place assignment with `res .= ` might be better than `res[:] =`

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [September 12, 2018, 2:57pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/7 "2018-09-12T14:57:05Z")

</div>

I think I tried that, and I don’t think the shapes matched.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 12, 2018, 3:19pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/8 "2018-09-12T15:19:35Z")

</div>

By the way, there are a lot of pretty simple things you can do to make this code much more performant (see [the performance tips](https://docs.julialang.org/en/stable/manual/performance-tips/#Avoid-global-variables-1) for more).

For reference the original version of your code (with the `res[:]` fix), runs in:

```julia
julia> @btime nlsolve($g!, $c_guess)
  3.965 ms (14859 allocations: 751.09 KiB)

```

1. Don’t use global variables inside your inner loop. Instead, you can create a function which takes in all of your parameters (gamma, alpha, nodes, etc.) and returns a new `g!` that will use those values internally (or, in other words, a closure):

```julia
function create_cost_function(fspace, nodes, γ, α, β)
   function f!(res, c_guess, fspace, nodes) # we want to find roots for f
     K = nodes
     con = funeval(c_guess,fspace,K)[1]
     K_next = K.^α - con
     con_next = funeval(c_guess,fspace,K_next)[1]
     res[:] = con.^(-γ) - β*con_next.^(-γ).*α.*K_next.^(α-1)
   end

   g!(res, c_guess) = f!(res, c_guess, fspace, nodes)
   return g!
end

g! = create_cost_function(fspace, nodes, γ, α, β)

```

which gives:

```julia
julia> @btime nlsolve($g!, $c_guess)
  2.836 ms (12824 allocations: 680.64 KiB)

```

or about a 30% improvement. Not much, but we’re just getting started.

1. Take advantage of broadcast fusion to update `res` completely in-place. For an explanation of what’s going on, check out [More Dots: Syntactic Loop Fusion in Julia](https://julialang.org/blog/2017/01/moredots) . In this case, I’m changing `res[:] = ` into `res .=` and adding dots everywhere on that line so that the entire assignment fuses into a single efficient loop with no extra allocated memory. To do this, I also have to make `con` and `con_next` into vectors, since `CompEcon` makes the very Matlab-like choice to return them as Nx1 matrices instead.

```julia
function create_cost_function(fspace, nodes, γ, α, β)
    function f!(res, c_guess, fspace, nodes) # we want to find roots for f
      K = nodes
      con = vec(funeval(c_guess,fspace,K)[1])
      K_next = K.^α - con
      con_next = vec(funeval(c_guess,fspace,K_next)[1])
      res .= con.^(-γ) .- β .* con_next.^(-γ) .* α .* K_next.^(α-1)
    end

    g!(res, c_guess) = f!(res, c_guess, fspace, nodes)
    return g!
end

g! = create_cost_function(fspace, nodes, γ, α, β)

```

Performance is now up to:

```julia
julia> @btime nlsolve($g!, $c_guess)
  415.725 μs (6004 allocations: 412.52 KiB)

```

or a solid 10X faster than the original code.

I suspect there’s a lot more that could be done, but it might require updating CompEcon’s internals.

---

<div class="post-metadata">

**Author:** ![IljaK91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iljak91/32/44301_2.png) [@IljaK91](https://discourse.julialang.org/u/IljaK91)\
**Post date:** [September 12, 2018, 3:35pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/9 "2018-09-12T15:35:30Z")

</div>

Wow! That is super useful. I will have a good look and update my code accordingly. I think the time has come when I switch to Julia for solving my root-finding problems 🙂

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 12, 2018, 3:38pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/10 "2018-09-12T15:38:27Z")

</div>

Great! If you want further performance improvements, it would be useful to try to pre-allocate your `con` and `con_next` variables, as in:

```julia
function create_cost(....)
  con = zeros(...)
  function f!(...)
    # Update `con` in-place here instead of creating a new vector
  end
end

```

Unfortunately, I don’t know how to tell `CompEcon.funeval` to update `con` in-place rather than creating a brand-new vector (that’s what I was referring to about possibly requiring some internal changes to CompEcon).

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [September 12, 2018, 3:44pm UTC](https://discourse.julialang.org/t/non-linear-root-finding-works-in-matlab-but-not-in-julia/14858/11 "2018-09-12T15:44:59Z")

</div>

This was a really helpful illustration. Thanks!