# Help with derivatives of matrix functions

**URL:** <https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437>\
**Category:** Machine Learning\
**Tags:** question\
**Created:** [February 14, 2022, 2:39pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437 "2022-02-14T14:39:55Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![raktim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raktim/32/11188_2.png) [@raktim](https://discourse.julialang.org/u/raktim)\
**Post date:** [February 14, 2022, 2:39pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/1 "2022-02-14T14:39:55Z")

</div>

Hi I have a matrix function G(x): \mathbb{R}^n \mapsto \mathbb{R}^{n\times n}, where x\in\mathbb{R}^n is a vector, and would like to use AD to compute derivatives \frac{\partial^2 G\_{ij}(x)}{\partial x\_i \partial x\_j}, where G\_{ij}(x) is the element of G(x).

This is how I am computing it:

```julia
∂(f, x, i) = ForwardDiff.partials(f([ForwardDiff.Dual{}(x[j], float(i==j)) for j in eachindex(x)]))[]; 

```

The j loop seems wasteful and I get slightly slower results with the following terrible approach:

```julia
X = [ForwardDiff.hessian(z->G(z)[i,j],x) for i = 1:2, j=1:2]

```

and picking the derivatives I need.

Is there a better way to do this? Any help here would be greatly appreciated.

Thanks!

---

<div class="post-metadata">

**Author:** ![StevenWhitaker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevenwhitaker/32/9749_2.png) [@StevenWhitaker](https://discourse.julialang.org/u/StevenWhitaker)\
**Post date:** [February 14, 2022, 4:10pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/2 "2022-02-14T16:10:30Z")

</div>

~~Sounds like you want [`ForwardDiff.jacobian`](https://juliadiff.org/ForwardDiff.jl/stable/user/api/#Jacobians-of-f(x::AbstractArray)::AbstractArray).~~

EDIT: I misread the question; `jacobian` is for first derivatives.

---

<div class="post-metadata">

**Author:** ![raktim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raktim/32/11188_2.png) [@raktim](https://discourse.julialang.org/u/raktim)\
**Post date:** [February 14, 2022, 4:20pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/3 "2022-02-14T16:20:44Z")

</div>

I want second derivative, which is different for each element. I fixed my question to indicate \frac{\partial^2 G\_{ij}(x)}{\partial x\_i \partial x\_j}.

---

<div class="post-metadata">

**Author:** ![CarloLucibello](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carlolucibello/32/3278_2.png) [@CarloLucibello](https://discourse.julialang.org/u/CarloLucibello)\
**Post date:** [February 14, 2022, 6:26pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/4 "2022-02-14T18:26:36Z")

</div>

Ideally you would like to construct a scalar-valued function out of G such that the hessian gives you the desired derivatives. Not sure how to do that though or if it’s even feasible

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 14, 2022, 6:35pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/5 "2022-02-14T18:35:33Z")

</div>

Here is how I’d do it:

```julia
using ForwardDiff
using BenchmarkTools

f(x) = x.^2 * transpose(x.^2) 

@show [ForwardDiff.hessian(x->f(x)[i, j], [3.0, 2.0, 1.0]) for i in 1:3, j in 1:3]
@show ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), [3.0, 2.0, 1.0])

@btime [ForwardDiff.hessian(x->f(x)[i, j], [3.0, 2.0, 1.0]) for i in 1:3, j in 1:3]
@btime ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), [3.0, 2.0, 1.0])

```

yielding

```julia
  8.800 μs (130 allocations: 39.62 KiB)
  2.256 μs (16 allocations: 5.41 KiB)

```

Interestingly enough the same implementation for `ReverseDiff`

```julia
ReverseDiff.jacobian(x -> ReverseDiff.jacobian(f, x), [3.0, 2.0, 1.0])

```

fails with

```julia
ERROR: LoadError: DimensionMismatch("matrix A has dimensions (3,3), matrix B has dimensions (1,3)")

```

---

<div class="post-metadata">

**Author:** ![raktim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raktim/32/11188_2.png) [@raktim](https://discourse.julialang.org/u/raktim)\
**Post date:** [February 15, 2022, 3:26am UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/6 "2022-02-15T03:26:44Z")

</div>

Thanks. That could work … but it still computes derivatives I do not need. For example, for a 2\times 2 problem, I am looking to efficiently compute \frac{\partial^2 G\_{11}(x)}{\partial x\_1^2}, \frac{\partial^2 G\_{12}(x)}{\partial x\_1\partial x\_2}, \frac{\partial^2 G\_{21}(x)}{\partial x\_2\partial x\_1}, and \frac{\partial^2 G\_{22}(x)}{\partial x\_2^2} only, and not the other derivatives. If the derivatives are computed over large number of grid points, the difference in computational time could be significant.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 15, 2022, 7:06am UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/7 "2022-02-15T07:06:40Z")

</div>

> [@raktim](#):
>
> That could work … but it still computes derivatives I do not need.

This computes fewer derivatives

```julia
g(i, j) = x -> (x.^2 * transpose(x.^2))[i, j]
dg(i, j, k) = x -> ForwardDiff.gradient(g(i, j), x)[k]
d2g(i, j, k, l) = x -> ForwardDiff.gradient(dg(i, j, k), x)[l]

@btime d2g(1, 1, 1, 1)([3.0, 2.0, 1.0])

h(i, j) = x -> x[i]^2 * x[j]^2
dh(i, j, k) = x -> ForwardDiff.gradient(h(i, j), x)[k]
d2h(i, j, k, l) = x -> ForwardDiff.gradient(dh(i, j, k), x)[l]

@btime d2h(1, 1, 1, 1)([3.0, 2.0, 1.0]) 

```

but yields a small improvement only

```julia
  1.940 μs (18 allocations: 3.97 KiB)
  1.890 μs (15 allocations: 1.88 KiB)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 25, 2022, 7:52am UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/8 "2022-02-25T07:52:58Z")

</div>

Hm, there might be a way do it somewhat faster in your special case rolling your own AD and using @elrod`s [lazy\_e](https://discourse.julialang.org/t/linear-algebra-unit-vectors-e-i/10409/15):

```julia
using StaticArrays, Random, Test, BenchmarkTools

struct Dual{T<:Number, G<:Number} <: Number
    val::T
    grad::G
    Dual{T, G}(val) where {T<:Number, G<:Number} = begin
        new{T, G}(val, zero(T))
    end
    Dual{T, G}(val, grad) where {T<:Number, G<:Number} = begin
        new{T, G}(val, grad)
    end
end

Dual(val::T) where {T<:Number} = Dual{T, T}(val)
Dual(val::T, grad::G) where {T<:Number, G<:Number} = Dual{T, G}(val, grad)

Base.convert(::Type{T}, d::Dual{T, G}) where {T<:Number, G<:Number} = d.val

import Base.+
+(d1::Dual{T, G1}, d2::Dual{T, G2}) where {T<:Number, G1<:Number, G2<:Number} = Dual(d1.val + d2.val, d1.grad + d2.grad)
+(d::Dual{T, G}, t::T) where {T<:Number, G<:Number} = Dual(d.val + t, d.grad)
+(t::T, d::Dual{T, G}) where {T<:Number, G<:Number} = d + t

import Base.-
-(d::Dual{T, G}) where {T<:Number, G<:Number} = Dual(-d.val, -d.grad)
-(d1::Dual{T, G1}, d2::Dual{T, G2}) where {T<:Number, G1<:Number, G2<:Number} = d1 + (-d2)
-(d::Dual{T, G}, t::T) where {T<:Number, G<:Number} = d + (-t)
-(t::T, d::Dual{T, G}) where {T<:Number, G<:Number} = t+ (-d)

import Base.*
*(d1::Dual{T, G1}, d2::Dual{T, G2}) where {T<:Number, G1<:Number, G2<:Number} = Dual(d1.val * d2.val, d1.val * d2.grad + d2.val * d1.grad)
*(d::Dual{T, G}, t::T) where {T<:Number, G<:Number} = Dual(d.val * t, t * d.grad)
*(t::T, d::Dual{T, G}) where {T<:Number, G<:Number} = d * t

import Base.inv
inv(d::Dual{T, G}) where {T<:Number, G<:Number} = Dual(one(T) / d.val, - d.grad / (d.val * d.val))

import Base./
/(d1::Dual{T, G1}, d2::Dual{T, G2}) where {T<:Number, G1<:Number, G2<:Number} = d1 * inv(d2)
/(d::Dual{T, G}, t::T) where {T<:Number, G<:Number} = d * inv(t)
/(t::T, d::Dual{T, G}) where {T<:Number, G<:Number} = t * inv(d)

struct lazy_e{T} <: AbstractVector{T}
    i::Int
    n::Int
end
lazy_e(i::Int, n::Int, ::Type{T} = Int) where T = lazy_e{T}(i, n)
function Base.getindex(e::lazy_e{T}, i) where T
    @boundscheck @assert i <= e.n
    ifelse(i == e.i, one(T), zero(T))
end
Base.size(e::lazy_e) = (e.n,)

derivative(f, x::T) where T = f(Dual(x, one(T))).grad
partial(f, i, x::SVector{N, T}) where {N, T} = f(map((x, d) -> Dual(x, d), x, lazy_e{T}(i, N))).grad

p(i, j) = x -> x[i]^2 * x[j]^2
dp(i, j, k) = x -> partial(p(i, j), k, x)
d2p(i, j, k, l) = x -> partial(dp(i, j, k), l, x)

@btime dp(1, 1, 1)(SVector(3.0, 2.0, 1.0))
@btime d2p(1, 1, 1, 1)(SVector(3.0, 2.0, 1.0))

```

yielding

```julia
  1.200 ns (0 allocations: 0 bytes)
  1.300 ns (0 allocations: 0 bytes)

```

which looks a bit like compile time AD to me 😉

---

<div class="post-metadata">

**Author:** ![raktim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raktim/32/11188_2.png) [@raktim](https://discourse.julialang.org/u/raktim)\
**Post date:** [February 25, 2022, 2:02pm UTC](https://discourse.julialang.org/t/help-with-derivatives-of-matrix-functions/76437/9 "2022-02-25T14:02:39Z")

</div>

This is fantastic! Perhaps should be a Package? Thank you @goerch!