# Automatic differentiation, adjoint methods, and linear systems of equations

**URL:** <https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979>\
**Category:** General Usage\
**Created:** [May 7, 2019, 9:28pm UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979 "2019-05-07T21:28:35Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![omalled](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omalled/32/12532_2.png) [@omalled](https://discourse.julialang.org/u/omalled)\
**Post date:** [May 7, 2019, 9:28pm UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979/1 "2019-05-07T21:28:35Z")

</div>

I am interested in differentiating things that are along the lines of what is described in section 2 of [these notes](https://math.mit.edu/~stevenj/18.336/adjoint.pdf) where the linear systems are sparse. For several years, I have occasionally been checking out Julia’s AD ecosystem to see if there is a package that can do this sort of differentiation automatically. I haven’t yet had success. The code below shows a couple of my feeble attempts to get [Zygote](https://github.com/FluxML/Zygote.jl) to do this for a very simple linear system. Am I doing something wrong, or is Zygote not ready for this yet? Is there another package that can automatically differentiate the code below? I am especially interested in the `g2` example which is more representative of real problems that I am working with.

```julia
import LinearAlgebra
import SparseArrays
import Zygote

function g1(p)
	n = 3
	A = LinearAlgebra.diagm(-1=>fill(-p, n - 1), 0=>fill(2p, n), 1=>fill(-p, n - 1))
	b = ones(n)
	return (A \ b)[1]
end

@show g1(1.0)
@show g1'(1.0)#Mutating arrays is not supported

function g2(p::T) where {T}
	n = 3
	I = Int[]
	J = Int[]
	V = T[]
	for i = 1:n
		if i > 1
			push!(I, i)
			push!(J, i - 1)
			push!(V, -p)
		end
		push!(I, i)
		push!(J, i)
		push!(V, 2p)
		if i < n
			push!(I, i)
			push!(J, i + 1)
			push!(V, -p)
		end
	end
	b = ones(n)
	return (SparseArrays.sparse(I, J, V) \ b)[1]
end

@show g2(1.0)
@show g2'(1.0)#segfault

```

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [May 7, 2019, 9:40pm UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979/2 "2019-05-07T21:40:31Z")

</div>

> [@omalled](#):
>
> along the lines of what is described in section 2 of [these notes](https://math.mit.edu/~stevenj/18.336/adjoint.pdf)

@stevengj , those are your notes 😄

---

<div class="post-metadata">

**Author:** ![jeff.bezanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff.bezanson/32/48_2.png) [@jeff.bezanson](https://discourse.julialang.org/u/jeff.bezanson)\
**Post date:** [May 7, 2019, 9:46pm UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979/3 "2019-05-07T21:46:57Z")

</div>

@MikeInnes any insights on if/when Zygote might be able to support this?

---

<div class="post-metadata">

**Author:** ![MikeInnes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikeinnes/32/3656_2.png) [@MikeInnes](https://discourse.julialang.org/u/MikeInnes)\
**Post date:** [May 8, 2019, 10:38am UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979/4 "2019-05-08T10:38:01Z")

</div>

Zygote should be able to handle this on the `mutate` branch (though it’ll work better if you can rewrite this in a more functional style).

If something crashes, it’s just a bug, and usually one that’s easy to fix; best to open an issue with this test case, simplified as much as you can, and I can take a look.

---

<div class="post-metadata">

**Author:** ![omalled](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omalled/32/12532_2.png) [@omalled](https://discourse.julialang.org/u/omalled)\
**Post date:** [May 8, 2019, 7:43pm UTC](https://discourse.julialang.org/t/automatic-differentiation-adjoint-methods-and-linear-systems-of-equations/23979/5 "2019-05-08T19:43:17Z")

</div>

Thanks for taking a look!

I tried these on the `mutate` branch. `g1` now gives a `MethodError: no method matching *(::Zygote.Stack{Any}, ::Int64)` and `g2` still segfaults. I also created a `g3` (see below) which reformulates `g2` in a more functional style (I think), but it still segfaults. I’ll open an issue on GitHub.

```julia
function g3(p)
	n = 3
	I = [2:n; 1:n; 1:n - 1]
	J = [1:n - 1; 1:n; 2:n]
	V = [fill(-p, n - 1); fill(2p, n); fill(-p, n - 1)]
	b = ones(n)
	return (SparseArrays.sparse(I, J, V) \ b)[1]
end

@show g3(1.0)
@show g3'(1.0)#segfault on mutate branch

```