# ReverseDiff.jl with iterative solvers

**URL:** <https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299>\
**Category:** General Usage\
**Created:** [June 16, 2017, 8:39am UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299 "2017-06-16T08:39:28Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 16, 2017, 8:39am UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/1 "2017-06-16T08:39:28Z")

</div>

I am trying to find the derivative of a function (that has input as the converged values of Jacobi iterator), with respect to the elements of the matrices used in Jacobi iterations. For example, if _ **f** (x)_ is a function of _x_, which in turn comes from **A** _x_ = **b** , I am trying to find the derivative of **f** with respect to element of **A** and **b**.

I am using **reverseDiff.jl** and **forwardDiff.jl** for this purpose, but all I can get is derivative of **f** with respect to the converged _x_, even when I pass **b** as parameter in the **GradientConfig** and **GradientTape** functions.

I have done this in C++ using CoDiPack library, but no success in Julia. Is this even possible in Julia? Any help would be appreciated. Thanks.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 16, 2017, 9:41am UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/2 "2017-06-16T09:41:16Z")

</div>

I’m not sure what you mean, but if x is the solution to Ax=b then you can compute analytically the derivative of x wrt A and b, which would probably be more efficient than reverse differentiation on the iteration. You can possibly even use more clever tricks like adjoint methods if you do this several times.

---

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 16, 2017, 12:16pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/3 "2017-06-16T12:16:09Z")

</div>

I can’t find it analytically, reason being, it is an aerodynamic application. There is fine distribution of points that define the shape of aerofoil. An iterative solver solves some differential equations iteratively with these points as inputs (which are floating point numbers) to get coefficients of drag and lift, and the sensitivities of these coefficients are computed with respect to each point on the surface of aerofoil.

I was using Jacobi as a test case, in which case, A and b would comprise of floating point numbers. Finite differences does not give resuts with desired precision.

It works with the C++ library I am using, just curious if it can be done in Julia.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 16, 2017, 12:33pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/4 "2017-06-16T12:33:43Z")

</div>

Your question is not very clear. Providing a MWE may help.

---

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 16, 2017, 12:53pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/5 "2017-06-16T12:53:31Z")

</div>

I am sorry for not being clear.  
If someone can shed some light on how to find derivative of f(x(b)) wrt b, with doing anything analytically.  
All I get is derivative of f(b) wrt b when using :  
b = 2  
x = b^2 + b;  
f(x) = sin(x) + x^2;  
result = ForwardDiff.derivative(f,b)

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [June 16, 2017, 1:36pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/6 "2017-06-16T13:36:20Z")

</div>

Chain rule?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 16, 2017, 1:39pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/7 "2017-06-16T13:39:05Z")

</div>

Possibly

```julia
using ForwardDiff
f(x) = sin(x) + x^2;
g(b) = b^2 + b;
b = 2
result = ForwardDiff.derivative(f ∘ g, b)

```

?

---

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 16, 2017, 2:51pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/8 "2017-06-16T14:51:29Z")

</div>

> [@Tamas\_Papp](#):
>
> f ∘ g, b

The symbol is throwing an error : ERROR: LoadError: UndefVarError: ∘ not defined.  
Am I missing something?

Anyways, here the code that I am actually trying to run, may give some insight :

```julia
using IterativeSolvers
using ReverseDiff: GradientConfig, GradientTape

A = [10.0 -1.0 2.0 0.0;
         -1.0 11.0 -1.0 3.0;
          2.0 -1.0 10.0 -1.0; 
          0.0 3.0 -1.0 8.0]
b = [6.0 ,25.0 ,-11.0 ,15.0]
x = similar(b)
f(x) = (x[1]^3) + log10((x[2]^2)) + exp(x[3]) + 4*x[4]
cfg = GradientConfig(b)
tape = GradientTape(f, b, cfg)

#x = jacobi(A, b)

jacobi!(x, A, b)

print(x)
println()
print(f(x))
println()

result = similar(x)
# ReverseDiff.gradient!(result, f, x, cfg)
# print(result)
# println()
ReverseDiff.gradient!(result, tape, b)
print(result)
println()

```

It gives derivative of f(b), not f(x(b)) with respect to b.  
Thanks.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 16, 2017, 2:54pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/9 "2017-06-16T14:54:23Z")

</div>

The composition symbol is defined only on Julia 0.6. You can use x → f(g(x)) instead.

---

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 16, 2017, 3:17pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/10 "2017-06-16T15:17:52Z")

</div>

Here is the same thing in forward mode :

```julia
using IterativeSolvers
using ForwardDiff: GradientConfig, gradient

A = [ 10.0 -1.0 2.0 0.0;
         -1.0 11.0 -1.0 3.0;
          2.0 -1.0 10.0 -1.0; 
          0.0 3.0 -1.0 8.0 ]
b = [6.0, 25.0, -11.0, 15.0]
x = [0,0,0,0]
f(x) = (x[1]^3) + log10((x[2]^2)) + exp(x[3]) + 4*x[4]
#cfg = GradientConfig(b)

x = jacobi(A,b)

print(x)
println()
print(f(x))
println()

result = ForwardDiff.gradient(f,b)
#result = ForwardDiff.derivative(f,b[1])

print(result)
println()

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 16, 2017, 3:29pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/11 "2017-06-16T15:29:19Z")

</div>

Independently from the technical aspects of getting this to work, my point was that if you want to compute the derivative of f with respect to b in this way, you have to compute the derivative of x with respect to b as an intermediate step, which is expensive. But because x solves a linear system, the derivative simplifies. Ie let F(b) = f(x(b)), you want to compute the gradient of F, which is J^T grad f, where J = A^-1 is the jacobian of x wrt b. Therefore, the correct way to compute this is to compute grad f, and solve the linear system A^T grad F = grad f by an iterative method.

---

<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:** [June 16, 2017, 3:59pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/12 "2017-06-16T15:59:53Z")

</div>

> [@antoine-levitt](#):
>
> Therefore, the correct way to compute this is to compute grad f, and solve the linear system A^T grad F = grad f by an iterative method.

Note that this is called the [the adjoint system of equations. The adjoint method](http://math.mit.edu/~stevenj/18.336/adjoint.pdf) is the same basic idea as reverse-mode auto-differentiation, but when you are using an iterative solver often auto-differentiation tools don’t work well (unless they know specifically about the solver being used), and then you have to apply the adjoint method “manually”.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 16, 2017, 4:32pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/13 "2017-06-16T16:32:30Z")

</div>

Thanks for the very nice write-up! It’s a useful trick that many people and communities use in some form or other and everybody should know about but is not often covered in a general setting in textbooks (Schur complements seem to be another instance of that)

---

<div class="post-metadata">

**Author:** ![jrevels](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jrevels/32/10393_2.png) [@jrevels](https://discourse.julialang.org/u/jrevels)\
**Post date:** [June 16, 2017, 6:37pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/14 "2017-06-16T18:37:35Z")

</div>

ForwardDiff/ReverseDiff differentiates native Julia functions at the input(s) you provide. When you call `ReverseDiff.gradient(g, v)` for some function `g` and some array `v`, you’re saying “return the gradient of `g` evaluated at `v`”.

> [@shubhamRanjan](#):
>
> It gives derivative of f(b), not f(x(b)) with respect to b.

Correct! It’s evaluating the gradient of `f` at the value `b`.

As @Tamas_Papp mentioned, it seems like you want to differentiate through some intermediary function that you’re not actually calling. Note that `x` (as you’ve defined it) isn’t a Julia function, it’s a `Vector{Float64}`.

Is the below closer to what you’d like to do?

```julia
using IterativeSolvers, ReverseDiff 

A = [10.0 -1.0 2.0 0.0;
     -1.0 11.0 -1.0 3.0;
      2.0 -1.0 10.0 -1.0; 
      0.0 3.0 -1.0 8.0]
b = [6.0 ,25.0 ,-11.0 ,15.0]

f(x) = (x[1]^3) + log10((x[2]^2)) + exp(x[3]) + 4*x[4]

# this is the function you seem to actually want to differentiate
g(A, b) = f(jacobi(A, b))

# get gradients w.r.t. `A` and `b`
dA, db = ReverseDiff.gradient(g, (A, b))

```

You can, of course, preallocate memory/prerecord the tape etc. to make the above more efficient, but before we go that far, I’m just making sure that this is closer to your intention.

---

<div class="post-metadata">

**Author:** ![shubhamRanjan](https://avatars.discourse-cdn.com/v4/letter/s/278dde/32.png) [@shubhamRanjan](https://discourse.julialang.org/u/shubhamRanjan)\
**Post date:** [June 23, 2017, 12:42pm UTC](https://discourse.julialang.org/t/reversediff-jl-with-iterative-solvers/4299/15 "2017-06-23T12:42:10Z")

</div>

Sorry for the late reply.  
This is actually what I wanted to do. Thanks for clearing it out. I was also able to do it slightly differently:

```julia
using IterativeSolvers
using ReverseDiff: GradientConfig, GradientTape

A = [10.0 -1.0 2.0 0.0;
         -1.0 11.0 -1.0 3.0;
          2.0 -1.0 10.0 -1.0; 
          0.0 3.0 -1.0 8.0]
b = [6.0 ,25.0 ,-11.0 ,15.0]
x = similar(b)
f(x) = (x[1]^3) + log10((x[2]^2)) + exp(x[3]) + 4*x[4]
cfg = GradientConfig(b)
tape = GradientTape(b -> f(jacobi(A,b)), b, cfg)

#x = jacobi(A, b)

jacobi!(x, A, b)

result = similar(x)

ReverseDiff.gradient!(result, tape, b)

```

Thanks everyone for helping out.