# Efficient Vector^T \* Jacobian matrix?

**URL:** <https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329>\
**Category:** Numerics\
**Created:** [March 1, 2019, 5:42am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329 "2019-03-01T05:42:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Christian\_Dengler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christian_dengler/32/20969_2.png) [@Christian\_Dengler](https://discourse.julialang.org/u/Christian_Dengler)\
**Post date:** [March 1, 2019, 5:42am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/1 "2019-03-01T05:42:23Z")

</div>

I’m trying to find a way to efficiently compute the product of a transposed vector with a Jacobian matrix. So given a vector function f and two large vectors v and x of same length, I want to compute v^T \* jacobian(f)(x) without explicitly computing the jacobian matrix. Does anyone know how this can be done? Apparently it can be done in some way similar to the jacobian-vector product jacobian(f)(x)\*v, but I’m unable to make the connection.

PS.: for the jacobian-vector product, the way to do it is shown in a paper called “Fast Exact Multiplication by the Hessian”, showing that jacobian(f)(x)_v = d/dr f(x+r_v)|r=0. Below my simple implementation of it, but as mentionned I need the other way around v^T \* jacobian(f)(x).

```julia
using ForwardDiff

# Compute jacobian(f)(x)*v without building a jacobian
function jvp(f::Function, x::AbstractVector, v::AbstractVector)
    g = r->f(x+r*v)
    return ForwardDiff.derivative(g, 0.0)
end

# I actually want v'*jacobian(f)(x)
vjp(f::Function, x::AbstractVector, v::AbstractVector) = ??

```

---

<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:** [March 1, 2019, 6:49am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/2 "2019-03-01T06:49:03Z")

</div>

This is exactly what reverse-mode automatic differentiation is (e.g. if f is R^n to R you’re asking for a O(1) gradient). There are a few packages doing it, but it’s more tricky than forward diff.

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [March 1, 2019, 6:50am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/3 "2019-03-01T06:50:05Z")

</div>

I think you can use something like

```julia
function vjp(f::Function, x::AbstractVector, v::AbstractVector)
    g(y) = dot(v, f(y))
    ForwardDiff.gradient(g, x)
end

```

Though reverse mode autodiff would be much more efficient unless your vector `x` is short.

[Edit: Although it looks nice, I think this formulation amounts to computing the full Jacobian internally (possibly in batches when `x` is long… but nevertheless).]

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [March 1, 2019, 7:25am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/4 "2019-03-01T07:25:09Z")

</div>

Yes it’s such a tricky problem that we have, in very rough order of age, and listed with comments / companion libraries:

- [ReverseDiffSource](https://github.com/JuliaDiff/ReverseDiffSource.jl) - via source transformation
- [ReverseDiff](https://github.com/JuliaDiff/ReverseDiff.jl) - via instruction tapes
- [ReverseDiffSparse / JuMP.Derivatives](https://github.com/mlubin/ReverseDiffSparse.jl) - JuMP
- [ReverseDiffTape](https://github.com/fqiang/ReverseDiffTape.jl) - JuMP?
- [AutoDiffSource](https://github.com/gaika/AutoDiffSource.jl) (?)
- [AutoGrad](https://github.com/denizyuret/AutoGrad.jl) - Knet
- [Flux.Tracker](https://github.com/FluxML/Flux.jl) - Flux
- [Nabla](https://github.com/invenia/Nabla.jl) - ML on models containing linear algebra
- [Zygote](https://github.com/FluxML/Zygote.jl) - Flux, next gen
- [Yötä](https://github.com/dfdx/Yota.jl) - ML model focused, backed by Cassette  
…
- [Capstan?](https://github.com/JuliaDiff/Capstan.jl) - Long awaited but stalled on compiler tech. Motivated awesome stuff like Cassette in the meantime.

Having listed all those, it’s hard to believe I’ve got them all…

---

<div class="post-metadata">

**Author:** ![Christian\_Dengler](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christian_dengler/32/20969_2.png) [@Christian\_Dengler](https://discourse.julialang.org/u/Christian_Dengler)\
**Post date:** [March 4, 2019, 5:41am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/5 "2019-03-04T05:41:56Z")

</div>

Hello, thanks for the tips! So as you posted the trick to use is that jacobian(f, x)\*v = d/dx (f(x)'\*v) .The ForwardDiff approach already gives a speedup for large vectors for some reason, and using AutoGrad for backward diff, the real speedup comes into play. So if someone else needs this:

```julia
using AutoGrad

# Compute jacobian(f, x)*v
function vjp(f::Function, x::AbstractVector, v::AbstractVector)
    x = Param(x)
    temp = @diff dot(v, f(x))
    return grad(temp, x)
end

```

I also benchmarked the functions agains the “create full jacobian and multiply” version for vectors of length 10000. and the vjp and jvp functions above.

```julia
using BenchmarkTools, ForwardDiff

f(x) = sum(x) .* x.^2
truejac(x) = ForwardDiff.jacobian(f, x)

x_test = randn(10000)
v_test = randn(10000)

# Hopefully slower functions, building the full jacobian
slow_jvp(x, v) = truejac(x)*v
slow_vjp(x, v) = vec(v'*truejac(x))

@benchmark slow_jvp($x_test, $v_test)
@benchmark jvp($f, $x_test, $v_test)

@benchmark slow_vjp($x_test, $v_test)
@benchmark vjp($f, $x_test, $v_test)

```

The results on my laptop are a 20000 times faster execution of jvp and of 7000 for vjp compared to the naive approach 🙂

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [March 4, 2019, 9:51am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/6 "2019-03-04T09:51:38Z")

</div>

> [@Christian\_Dengler](#):
>
> The ForwardDiff approach already gives a speedup for large vectors for some reason

I think it might just be that the jacobian doesn’t need to be stored all at once because the columns are computed in batches with `ForwardDiff` and each column can be reduced with `v` as it’s generated. So improved cache efficiency, etc.

---

<div class="post-metadata">

**Author:** ![richardr2926](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richardr2926/32/207243_2.png) [@richardr2926](https://discourse.julialang.org/u/richardr2926)\
**Post date:** [September 11, 2024, 4:01am UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/7 "2024-09-11T04:01:34Z")

</div>

Hello.

This is what pullback in Zygote is meant to do. So:

```julia
function f(x)
    output = [
        x[1]^6 * x[2]^4 * x[3]^9 * x[4]^2;
        x[1]^2 * x[2]^3 * x[3]^5 * x[4]^3;
        x[1]^5 * x[2]^7 * x[3]^7 * x[4]^6;
    ]
    return output
end

left_multiplication_point = [0.5; 0.8; 1.0]
evaluation_point = [0.2; 0.1; 1.0; 4.0]

function clever_vjp(func, primal, cotangent)
    _, func_pullback = pullback(func, primal)
    vjp_result, = func_pullback(cotangent)
    return vjp_result
end

clever_vjp(f, evaluation_point, left_multiplication_point)

```

---

<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:** [September 11, 2024, 5:41pm UTC](https://discourse.julialang.org/t/efficient-vector-t-jacobian-matrix/21329/8 "2024-09-11T17:41:37Z")

</div>

Thanks for your answer, but in general it’s best not to ressuscitate threads from 6 years ago, because everyone gets pinged in the process. If you have specific questions you can create a new topic for that 😉  
Also do check out DifferentiationInterface.jl with its `pullback` function. Nowadays that’s probably the most generic solution.