# ForwardDiff with Matrix of Vector as input argument
**URL:** https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049
**Category:** Specific Domains
**Tags:** question
**Created:** [October 8, 2024, 6:46am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049 "2024-10-08T06:46:19Z")
**Posts on this page:** 18
**Page:** 2
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 9, 2024, 9:47pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/21 "2024-10-09T21:47:11Z")
Also note that my package DifferentiationInterface.jl is meant to make your life easier in such situations, by hiding away the details of preallocation and performance optimization. Take a look at the [tutorial](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorials/basic/) if you’re interested.
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 11, 2024, 2:21pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/22 "2024-10-11T14:21:38Z")
First tests with foo!() instead of foo() fail. Resulting gradient is a zero matrix.
Should ForwardDiff.jacobian() be replaced by ForwardDiff.jacobian!() as well?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 11, 2024, 3:21pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/23 "2024-10-11T15:21:05Z")
Can you provide an MWE?
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 11, 2024, 5:32pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/24 "2024-10-11T17:32:29Z")
Here is an attempt
```julia
function foo!(y,x)
y = sum(x)*ones(length(x))
return y
end
x1 = ones(10); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
ForwardDiff.jacobian(foo!, y1, x1)
```
yields as output
```julia
10×10 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 11, 2024, 5:43pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/25 "2024-10-11T17:43:59Z")
The issue here is that you’re not modifying the argument `y` passed to the function in-place. You’re creating a new object (which happens to have the same name) which you fill with the desired output. What you need to do is mutate the very `y` that the function receives.
A simple fix is the following:
```julia
function foo!(y,x)
# notice the dot to denote elementwise assignment
y .= sum(x)*ones(length(x))
return y
end
```
but a more clever one would be
```julia
function foo!(y, x)
y .= sum(x)
end
```
You may want to take a look at the docs page on performance tips to improve such in-place functions and maybe avoid allocations altogether
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 11, 2024, 7:21pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/26 "2024-10-11T19:21:21Z")
Thx again.
Now obtaining correct Jacobian with in-place functions.
Improvements in terms of memory allocations remain only marginal.
**Version 1 - Using Jacobian (no pound)**
```julia
function foo!(y,x)
# need to mutate the very y that the function receives
# notice the dot to denote elementwise assignment
y .= sum(x)*ones(length(x))
return y
end
```
```julia
x1 = ones(10); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
cfg = ForwardDiff.JacobianConfig(foo!,y1,x1)
@btime ForwardDiff.jacobian(foo!, y1, x1, cfg)
```
resulting in
```julia
528.361 ns (4 allocations: 2.88 KiB)
```
**Version 2 - Using Jacobian! (with pound)**
```julia
x1 = ones(10); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
cfg = ForwardDiff.JacobianConfig(foo!,y1,x1)
@btime ForwardDiff.jacobian!(DiffResults.JacobianResult(y1,x1),foo!,y1,x1,cfg)
```
resulting in
```julia
600.000 ns (6 allocations: 3.05 KiB)
```
More soon.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [October 11, 2024, 7:35pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/27 "2024-10-11T19:35:17Z")
> [@Domenico\_Lahaye](#):
>
> Improvements in terms of memory allocations remain only marginal.
>
> ```julia
> function foo!(y,x)
> # need to mutate the very y that the function receives
> # notice the dot to denote elementwise assignment
> y .= sum(x)*ones(length(x))
> return y
> end
>
> ```
The call to `ones` allocates a new array. It is unnecessary here, as broadcasting the assignment takes care of writing into all of `y` anyways:
```julia
function baz!(y, x)
y .= sum(x)
end
```
```julia-repl
# Better to interpolate global variables when benchmarking to avoid spurious allocations
julia> @btime foo!($y1,$x1);
143.241 ns (2 allocations: 288 bytes)
julia> @btime baz!($y1,$x1);
16.323 ns (0 allocations: 0 bytes)
```
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 11, 2024, 7:51pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/28 "2024-10-11T19:51:18Z")
Thx!
That drastically reduces the number of allocations (hurray!).
The CPU time, however, remains roughly equal to earlier more naive implementations (see Post # 18 above).
```julia
function foo!(y,x)
# need to mutate the very y that the function receives
# notice the dot to denote elementwise assignment
y .= sum(x)
return y
end
```
```julia
x1 = ones(10); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
cfg = ForwardDiff.JacobianConfig(foo!,y1,x1)
@btime ForwardDiff.jacobian(foo!, $y1, $x1, cfg)
```
```julia
432.161 ns (2 allocations: 1.75 KiB)
```
```julia
x1 = ones(100); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
cfg = ForwardDiff.JacobianConfig(foo!,y1,x1)
@btime ForwardDiff.jacobian(foo!, $y1, $x1, cfg)
```
```julia
24.792 μs (5 allocations: 80.30 KiB)
```
```julia
x1 = ones(1000); y1 = zeros(length(x1))
y1 = foo!(y1,x1)
cfg = ForwardDiff.JacobianConfig(foo!,y1,x1)
@btime ForwardDiff.jacobian(foo!, $y1, $x1, cfg)
```
```julia
2.425 ms (3 allocations: 7.63 MiB)
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 11, 2024, 9:52pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/29 "2024-10-11T21:52:00Z")
Try it again while interpolating (putting a `$` on) the `cfg` object too, since it is also a global variable
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 12, 2024, 5:32am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/30 "2024-10-12T05:32:05Z")
Hardly changes
```julia
287.826 ns (1 allocation: 896 bytes) # problem size 10
24.666 μs (4 allocations: 79.08 KiB) # problem size 100
2.399 ms (2 allocations: 7.63 MiB) # problem size 1000
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 12, 2024, 10:58am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/31 "2024-10-12T10:58:18Z")
For optimal performance, you also want to pre-allocate the output matrix:
```julia
using ForwardDiff, BenchmarkTools
foo!(y, x) = y .= sum(x)
function benchmark_trivial_jacobian(n)
x = ones(n)
y = similar(x)
J = similar(y, length(y), length(x))
cfg = ForwardDiff.JacobianConfig(foo!, y, x)
@btime ForwardDiff.jacobian(foo!, $y, $x, $cfg)
@btime ForwardDiff.jacobian!($J, foo!, $y, $x, $cfg)
return nothing
end
```
But it doesn’t change much for larger problems, except for memory use:
```julia
julia> benchmark_trivial_jacobian(10)
370.061 ns (1 allocation: 896 bytes)
304.813 ns (0 allocations: 0 bytes)
julia> benchmark_trivial_jacobian(100)
27.706 μs (19 allocations: 89.31 KiB)
26.636 μs (17 allocations: 11.14 KiB)
julia> benchmark_trivial_jacobian(1000)
3.515 ms (169 allocations: 7.73 MiB)
3.423 ms (167 allocations: 104.89 KiB)
```
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 13, 2024, 6:43am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/32 "2024-10-13T06:43:46Z")
Thanks. I much appreciate all your input.
I would like to better understand why the version with all the bells and whistles (buffer on input, preallocation of output) is not faster than the more naive implementation.
Stated differently, I wonder why the reduction of number of allocations does not results in a reduction of CPU time.
Is this due to the Jacobian being dense? Is looking into a MWE with a sparse Jacobian worthwhile here?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 13, 2024, 6:55am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/33 "2024-10-13T06:55:45Z")
At some point you’re limited by
- the time it takes your function to run
- the overhead caused by autodiff
If you know for a fact that your Jacobian is sparse, then yes you can accelerate it significantly by making use of that knowledge. See the [DifferentiationInterface sparse tutorial](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorials/advanced/#Sparsity) for an example.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 13, 2024, 7:04am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/34 "2024-10-13T07:04:12Z")
It may also be worth pondering whether you need a full Jacobian at all. For instance, if you compute `J` only to take its products with some vector, there are faster ways to achieve this
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 13, 2024, 7:00pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/35 "2024-10-13T19:00:37Z")
Thanks again.
The Jacobian in the application I target is dense. The application requires to solve a linear system with the Jacobian.
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 16, 2024, 8:35pm UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/36 "2024-10-16T20:35:11Z")
Does using ForwardDiff on a function that calls hcubature with buffer as argument cause particular challenges?
Without buffer as argument things works fine.
With buffer I get a type conversion error.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [October 17, 2024, 6:24am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/37 "2024-10-17T06:24:39Z")
You’re gonna need PreallocationTools.jl to automatically generate buffers that have the right type
---
### Author:  [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)
#### Post date: [October 17, 2024, 6:52am UTC](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049/38 "2024-10-17T06:52:23Z")
Aha! Will check it out. Thx again.
[Previous page](https://discourse.julialang.org/t/forwarddiff-with-matrix-of-vector-as-input-argument/121049.md?page=1)