# Wrong gradient from Zygote?

**URL:** https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204
**Category:** Numerics
**Tags:** question, zygote
**Created:** [November 14, 2023, 12:43pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204 "2023-11-14T12:43:21Z")
**Posts on this page:** 14
**Page:** 1

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 12:43pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/1 "2023-11-14T12:43:21Z")

</div>

I have the following example where I get different gradient values for forward and backward:

```julia
using Random, ChainRulesCore, Zygote, Statistics
function simulate(N, drift, sigma, rng,i)
    ChainRulesCore.@ignore_derivatives randn!(rng,N)
    return @. sigma * N + drift
end

function simulate_bs(S0, r, T, sigma, d, nsim::Integer, nsteps::Integer, rng)
    mu = r - d
    zero_dual =ChainRulesCore.@ignore_derivatives 0 * mu * sigma * T
    dt = T / nsteps
    mu_adj = (mu - sigma^2 / 2) * dt
    sigma_adj = sigma * sqrt(dt)
    X =zeros(typeof(zero_dual), nsim)
    N =Array{Float64}(undef,nsim)
    ChainRulesCore.@ignore_derivatives Random.seed!(rng, 1)
    for i = 1:nsteps
        X += simulate(N, mu_adj, sigma_adj, rng,i)
    end
    @. S0 * exp(X)
end

function pricer_bs(S0, r, T, sigma, d, nsim::Integer, nsteps::Integer, rng)
    X_ = simulate_bs(S0, r, T, sigma, d, nsim, nsteps, rng)
    return mean(X_) * exp(-r * T)
end
function pricer_zygote_f(S0, r, T, sigma, d)
	Nsim1 = 100_000;
	Nstep1 = 3;
	rng1 = MersenneTwister()
	return pricer_bs(S0, r, T, sigma, d, Nsim1, Nstep1, rng1)
end

const S0 = 100.0;
const r = 0.01;
const T = 1.0;
const d = 0.01;
const sigma = 0.2;

fwd_grad=collect(Zygote.forward_jacobian(x->pricer_zygote_f(x...), [S0,r, T, sigma, d])[2])
rev_grad=collect(Zygote.gradient(pricer_zygote_f, S0,r, T, sigma, d))
@show fwd_grad.-rev_grad

```

and I get some non negligible mismatching on the derivative on T and sigma.  
Do you have any idea of what is going on?

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 12:44pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/2 "2023-11-14T12:44:00Z")

</div>

This is the difference I get between forward and backward:

```julia
-1.1102230246251565e-16
  2.842170943040401e-14
  0.014689082933301956
  0.1468908293330241
 -2.842170943040401e-14

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 14, 2023, 3:49pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/3 "2023-11-14T15:49:24Z")

</div>

Can you try with ForwardDiff? My guess is that `Zygote.forward_jacobian` is the issue. That function isn’t in the documentation for a reason.

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 3:53pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/4 "2023-11-14T15:53:30Z")

</div>

I tried with DualNumbers and I get the same result as forward\_jacobian

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 3:55pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/5 "2023-11-14T15:55:17Z")

</div>

Even with ForwardDiff (the difference is against fwd\_grad):

```julia
ForwardDiff.gradient(x->pricer_zygote_f(x...), [S0,r, T, sigma, d]).-fwd_grad

```

I get:

```julia
 0.0
 0.0
 0.0
 0.0
 0.0

```

---

<div class="post-metadata">

### Author: ![ToucheSir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/touchesir/32/14411_2.png) [@ToucheSir](https://discourse.julialang.org/u/ToucheSir)
#### Post date: [November 14, 2023, 4:38pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/6 "2023-11-14T16:38:36Z")

</div>

I suspect differences between how both libraries operate are leading to the RNG being called and possibly seeded a different number of times. e.g. the number of pushforward calls far exceeds the number of pullback ones. If you can try (just for testing) generating `N` once and fixing it for the duration of the simulation, you should be able to isolate if this number of RNG calls discrepancy is the culprit.

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 14, 2023, 4:48pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/7 "2023-11-14T16:48:24Z")

</div>

> [@ToucheSir](#):
>
> If you can try (just for testing) generating `N` once and fixing it for the duration of the simulation, you should be able to isolate if this number of RNG calls discrepancy is the culprit.

Yep, that appears to be the culprit:

```julia-repl
julia> const N = rand(100_000);

julia> function simulate_bs(S0, r, T, sigma, d, nsim::Integer, nsteps::Integer, rng)
           mu = r - d
           zero_dual =ChainRulesCore.@ignore_derivatives 0 * mu * sigma * T
           dt = T / nsteps
           mu_adj = (mu - sigma^2 / 2) * dt
           sigma_adj = sigma * sqrt(dt)
           X =zeros(typeof(zero_dual), nsim)
           #N =Array{Float64}(undef,nsim)
           ChainRulesCore.@ignore_derivatives Random.seed!(rng, 1)
           for i = 1:nsteps
               X += simulate(N, mu_adj, sigma_adj, rng,i)
           end
           @. S0 * exp(X)
       end;

julia> function simulate(N, drift, sigma, rng,i)
           #ChainRulesCore.@ignore_derivatives randn!(rng,N)
           return @. sigma * N + drift
       end;

julia> fwd_grad=collect(Zygote.forward_jacobian(x->pricer_zygote_f(x...), [S0,r, T, sigma, d])[2]);

julia> rev_grad=collect(Zygote.gradient(pricer_zygote_f, S0,r, T, sigma, d));

```

```julia-repl
julia> @show fwd_grad.-rev_grad
fwd_grad .- rev_grad = [2.220446049250313e-16; 1.4210854715202004e-14; -3.552713678800501e-15; -2.842170943040401e-14; -2.842170943040401e-14;;]
5×1 Matrix{Float64}:
  2.220446049250313e-16
  1.4210854715202004e-14
 -3.552713678800501e-15
 -2.842170943040401e-14
 -2.842170943040401e-14

```

---

<div class="post-metadata">

### Author: ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [November 14, 2023, 6:02pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/8 "2023-11-14T18:02:42Z")

</div>

What this means is that taking gradients of a function that includs randomness is tricky business. This article helped me understand things a bit better:

> **[Monte Carlo Gradient Estimation in Machine Learning](https://arxiv.org/abs/1906.10652)**
>
> This paper is a broad and accessible survey of the methods we have at our disposal for Monte Carlo gradient estimation in machine learning and across the statistical sciences: the problem of computing the gradient of an expectation of a function with...

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 14, 2023, 8:46pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/9 "2023-11-14T20:46:11Z")

</div>

> [@gdalle](#):
>
> What this means is that taking gradients of a function that includs randomness is tricky business.

See also [this lecture](https://ocw.mit.edu/courses/18-s096-matrix-calculus-for-machine-learning-and-beyond-january-iap-2023/resources/ocw_18s096_lecture07-part1_2023feb01_mp4/) and the [notes](https://ocw.mit.edu/courses/18-s096-matrix-calculus-for-machine-learning-and-beyond-january-iap-2023/pages/lecture-notes/#lecture-7) from our matrix-calculus class. For stochastic gradient descent (SGD) and similar algorithms (e.g. Adam), you only need a function whose _expected value_ is the gradient of the expected value of your function.

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 10:44pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/10 "2023-11-14T22:44:43Z")

</div>

it seems to me a more Zygote related issue. The number returned is simply wrong.

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 14, 2023, 10:50pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/11 "2023-11-14T22:50:23Z")

</div>

No, this is just a difference between forwards mode autodiff and reverse mode autodiff. The code path that calls `rand` is hit a different number of times when you run in reverse mode versus when you run in forwards mode.

This would happen with any autodiff system with both modes.

This is kinda the whole reason reverse mode exists in the first place. For a function with many inputs and few outputs, reverse mode can calculate the gradient with fewer function evaluations than forwards mode.

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 11:00pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/12 "2023-11-14T23:00:28Z")

</div>

The randn function should be ignored in reverse mode no?

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 14, 2023, 11:08pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/13 "2023-11-14T23:08:41Z")

</div>

In the derivative part, yes, but the primal (the non-gradient version) still gets called during the construction of the gradient, and the number of times it gets called is different depending on forwards or backwards.

Here’s a demonstration with a simple counter in a `@ignore_derivatives` block:

```julia-repl
julia> using Zygote, ChainRulesCore

julia> const counter = Ref(0);

julia> function foo(x::Vector)
           ChainRulesCore.@ignore_derivatives counter[] += 1
           sum(x)
       end;

julia> Zygote.gradient(foo, rand(100));

julia> counter[]
1

```

So when we calculate the gradient of `foo` here, we increment the counter only once.

Now lets reset the counter and look at forwards mode:

```julia-repl
julia> counter[] = 0;

julia> Zygote.forward_jacobian(foo, rand(100));

julia> counter[]
9

```

So you see `foo` has been called `9` times!

---

<div class="post-metadata">

### Author: ![rcalxrc08](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rcalxrc08/32/14973_2.png) [@rcalxrc08](https://discourse.julialang.org/u/rcalxrc08)
#### Post date: [November 14, 2023, 11:39pm UTC](https://discourse.julialang.org/t/wrong-gradient-from-zygote/106204/14 "2023-11-14T23:39:58Z")

</div>

Thank you!  
Shouldn’t my example be more similar to the following?

```julia
using Zygote, ChainRulesCore

const counter = Ref(0);

function incr_counter()
	ChainRulesCore.@ignore_derivatives counter[] += 1
end

function foo(x::Vector)
		ChainRulesCore.@ignore_derivatives counter[]=10#as a seed
		for i in x
			incr_counter()
		end
       sum(x)
       end

```
