# How to use ForwardDiff.jl with inv and matrix exponential?

**URL:** <https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880>\
**Category:** General Usage\
**Tags:** question, package, forwarddiff, autodiff\
**Created:** [June 30, 2026, 5:38pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880 "2026-06-30T17:38:01Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [June 30, 2026, 5:38pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/1 "2026-06-30T17:38:01Z")

</div>

I have a question listed here: [Why does ForwardDiff.jl give me all-zero Jacobian matrix?](https://discourse.julialang.org/t/why-does-forwarddiff-jl-give-me-all-zero-jacobian-matrix/137714), and I found it might be due to part of the function which contains `inv` and `exp` (matrix exponential…) from `ExponentialUtilities`. I am not sure whether ForwardDiff.jl could fail with `inv` of matrix or `exp`?

one example of my function could see below:

```julia-auto
# ------------------------------------------------------------------
# solve_odes_water: matrix-exponential solution of dW/dt = A × W + b.
#
# Analytical solution over one timestep ΔT:
# W(t+ΔT) = e^(A·ΔT) · (W₀ + A⁻¹·b) − A⁻¹·b
#
# Derivation: substituting u = W + A⁻¹b → du/dt = A·u → u(t) = e^(At)·u₀
# → W(t+ΔT) = e^(A·ΔT)·(W₀ + A⁻¹·b) − A⁻¹·b
#
# Falls back to W(t+ΔT) = W₀ when A is singular (isSingular_water = true),
# which avoids NaN propagation at very low carbon / dry spinup.
#
# Mathematically identical to solve_odes_3_water in original.
# ------------------------------------------------------------------
function solve_odes_water(new_vegW, deriv_mat, deriv_const, init_cond, deltaT, z_zero, helpers)
    if isSingular_water(deriv_mat, z_zero)
        for ix ∈ eachindex(new_vegW)
            @rep_elem init_cond[ix] ⇒ (new_vegW, ix, :vegW)
        end
    else
        inv_deriv = inv(deriv_mat)
        new_vegW = real(-inv_deriv * deriv_const +
                         exp(deltaT * deriv_mat) * (init_cond + inv_deriv * deriv_const))
    end
    return new_vegW
end

```

Btw, in my code/models, I am using static arrays and static matrix to reduce allocations, and at some places I forced to use Float32. In the final cost function, I set it as dual type. I guess this is not the problem? because I run same codes with different parameter values and input time series data (parameter values coming from optimization over different sites using different input data) the model over 40 independent sites (40 independent runs), only 5-8 sites return zero-valued Jacobians…others are fine…

---

<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 30, 2026, 9:19pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/2 "2026-06-30T21:19:07Z")

</div>

> [@Xu\_Shan](#):
>
> I am not sure whether [ForwardDiff.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ForwardDiff) could fail with `inv` of matrix or `exp`?

It’s easy to try — if you set `A = rand(3,3)`, you can easily see that `ForwardDiff.jacobian(inv, A)` succeeds, and indeed matches the [analytical result](https://github.com/mitmath/matrixcalc/) (-A^T \otimes A):

```julia-auto
julia> ForwardDiff.jacobian(inv, A) ≈ -kron(inv(A)', inv(A))
true

```

but with the matrix exponential it throws a `MethodError`:

```julia-auto
julia> ForwardDiff.jacobian(exp, A)
ERROR: MethodError: no method matching exp!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{typeof(exp), Float64}, Float64, 9}})

```

and indeed there is an open issue for this: [support for `Base.expm` (need advice implementing) · Issue #174 · JuliaDiff/ForwardDiff.jl · GitHub](https://github.com/JuliaDiff/ForwardDiff.jl/issues/174)

(ChainRules.jl has rules for the matrix exponential, however, so if you use something like Zygote.jl or Enzyme.jl that supports ChainRules it should be easier to get it working. Or use a different implementation of `exp` like `expv` in ExponentialUtilities.jl)

---

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [July 1, 2026, 7:34am UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/3 "2026-07-01T07:34:41Z")

</div>

> [@stevengj](#):
>
> Or use a different implementation of `exp` like in [ExponentialUtilities.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ExponentialUtilities)

Yes I am using `exp` from `ExponentialUtilities.jl`…but seems to be broken too…

```julia-auto
julia> ForwardDiff.jacobian(ExponentialUtilities.exp, Float32.(A))
ERROR: MethodError: no method matching exp!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{typeof(exp), Float32}, Float32, 9}})
The function `exp!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  exp!(::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ LinearAlgebra /viper/ptmp1/xshan/julia/depot/juliaup/julia-1.11.6+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:677

Stacktrace:
 [1] exp(A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{typeof(exp), Float32}, Float32, 9}})
   @ LinearAlgebra /viper/ptmp1/xshan/julia/depot/juliaup/julia-1.11.6+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:622
 [2] vector_mode_dual_eval!(f::typeof(exp), cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float32, 9, Matrix{…}}, x::Matrix{Float32})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/z3jRk/src/apiutils.jl:24
 [3] vector_mode_jacobian(f::typeof(exp), x::Matrix{Float32}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float32, 9, Matrix{…}})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/z3jRk/src/jacobian.jl:129
 [4] jacobian(f::typeof(exp), x::Matrix{Float32}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float32, 9, Matrix{…}}, ::Val{true})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/z3jRk/src/jacobian.jl:22
 [5] jacobian(f::typeof(exp), x::Matrix{Float32}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(exp), Float32}, Float32, 9, Matrix{ForwardDiff.Dual{…}}})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/z3jRk/src/jacobian.jl:19
 [6] jacobian(f::typeof(exp), x::Matrix{Float32})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/z3jRk/src/jacobian.jl:19
 [7] top-level scope
   @ REPL[13]:1
 [8] top-level scope
   @ none:1
Some type information was truncated. Use `show(err)` to see complete types.

```

---

<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:** [July 1, 2026, 11:21am UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/4 "2026-07-01T11:21:48Z")

</div>

> [@Xu\_Shan](#):
>
> Yes I am using `exp` from `ExponentialUtilities.jl`…but seems to be broken too…

For ExponentialUtilities you would use `expv` for the exponential–vector product. But that doesn’t seem to support ForwardDiff either.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [July 1, 2026, 3:19pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/5 "2026-07-01T15:19:07Z")

</div>

Also open in issue [#51008](https://github.com/JuliaLang/julia/issues/51008).

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [July 1, 2026, 4:48pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/6 "2026-07-01T16:48:34Z")

</div>

I’ve needed this before too so I slammed together my own implementation of \phi\_1(x) = \frac{\exp(x)-1}{x}. It should work with `StaticArrays` and `ForwardDiff`. No warranty but here it is:

```julia
using LinearAlgebra

"""
    expm1o(A)

Compute `(exp(A) - I) / A` via a division-free calculation that is suitable for singular `A`.
"""
function expm1o(A)
    # note: this functionality is available as `last(ExponentialUtilities.phi(A, 1))`
    # see https://www.maths.uq.edu.au/expokit/paper.pdf or https://www.sciencedirect.com/science/article/pii/S0377042720300492 for allegedly better ways to compute it

    function maybematrixpoly(x::Union{AbstractMatrix,UniformScaling},p)
        if eltype(p) <: Number
            return evalpoly(x,UniformScaling.(p))
        else
            return evalpoly(x,p)
        end
    end
    maybematrixpoly(x,p) = evalpoly(x,p)

	nA = norm(A,2) # frobenius norm
    T = typeof(nA)
	nA = nA > floatmin(T) ? nA : floatmin(T) # avoid issue with degenerate nA
	ex = ceil(Int, log2(nA))
	logscale = max(4 + ex, 0)
    p = inv.(cumprod(ntuple(T, Val(9)))) # a Pade approximation would probably be faster/more accurate
	scale = exp2(T(-logscale))
	sA = A * scale
	expmsA = maybematrixpoly(sA, p)
	expmA = expmsA
	for _ in 1:logscale
		expmA = (expmA * expmA) * A * (scale / 2) + expmA
		scale *= 2
	end
	return expmA
end

```

You should be able to write your implementation in terms of this (and probably don’t need to branch on singular matrices any more).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 2, 2026, 10:20am UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/7 "2026-07-02T10:20:06Z")

</div>

> [@Xu\_Shan](#):
>
> Yes I am using `exp` from `ExponentialUtilities.jl`…but seems to be broken too…
> 
> ```julia-auto
> 
> ```

Use `exp_generic`. That is tested to work for this.

> <https://github.com/SciML/ExponentialUtilities.jl/blob/v1.30.1/test/basictests.jl#L157-L227>

---

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [July 2, 2026, 11:24am UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/8 "2026-07-02T11:24:53Z")

</div>

> [@ChrisRackauckas](#):
>
> `Matrix(reshape((1.0:(23.0^2)) ./ 700, (23, 23)))`

Not working…I am using `StaticMatrix` and `StaticArray`…this is required by the model to create zero allocations…

```julia-auto
julia> B = dt * A
4×4 SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16} with indices SOneTo(4)×SOneTo(4):
 Dual{Nothing}(-55.1042,2.67628e20) Dual{Nothing}(28.9682,4.77492e22) … Dual{Nothing}(0.0,0.0)
  Dual{Nothing}(55.1042,-2.67628e20) Dual{Nothing}(-57.9383,-9.55121e22) Dual{Nothing}(0.0,0.0)
   Dual{Nothing}(0.0,0.0) Dual{Nothing}(28.9701,4.77628e22) Dual{Nothing}(33.8189,-6.4245e19)
   Dual{Nothing}(0.0,0.0) Dual{Nothing}(0.0,0.0) Dual{Nothing}(-67.6425,5.67798e19)

julia> E_base = exp(B)
4×4 SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16} with indices SOneTo(4)×SOneTo(4):
 Dual{Nothing}(0.00445987,NaN) Dual{Nothing}(0.00414118,NaN) Dual{Nothing}(0.00325967,NaN) Dual{Nothing}(0.00173045,NaN)
 Dual{Nothing}(0.00787748,NaN) Dual{Nothing}(0.00731459,NaN) Dual{Nothing}(0.00575757,NaN) Dual{Nothing}(0.0030565,NaN)
 Dual{Nothing}(0.00906899,NaN) Dual{Nothing}(0.00842096,NaN) Dual{Nothing}(0.00662844,NaN) Dual{Nothing}(0.00351882,NaN)
 Dual{Nothing}(0.00282028,NaN) Dual{Nothing}(0.00261875,NaN) Dual{Nothing}(0.00206131,NaN) Dual{Nothing}(0.00109428,NaN)

julia> E_generic = exp_generic(B)
ERROR: setindex!(::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] setindex!(a::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}, value::ForwardDiff.Dual{Nothing, Float64, 1}, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/0cEwi/src/indexing.jl:3
  [3] macro expansion
    @ ~/.julia/packages/StaticArrays/0cEwi/src/broadcast.jl:159 [inlined]
  [4] _broadcast!
    @ ~/.julia/packages/StaticArrays/0cEwi/src/broadcast.jl:143 [inlined]
  [5] _copyto!
    @ ~/.julia/packages/StaticArrays/0cEwi/src/broadcast.jl:70 [inlined]
  [6] copyto!
    @ ~/.julia/packages/StaticArrays/0cEwi/src/broadcast.jl:63 [inlined]
  [7] materialize!
    @ ./broadcast.jl:883 [inlined]
  [8] materialize!
    @ ./broadcast.jl:880 [inlined]
  [9] exp_gen!(cache::Vector{MMatrix{4, 4, ForwardDiff.Dual{…}, 16}}, A::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}, ::Val{10})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/mb4U7/src/exp_generated/exp_10.jl:22
 [10] exponential!(A::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}, method::ExpMethodHigham2005, _cache::T)
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/mb4U7/src/exp_noalloc.jl:106
 [11] exponential!(A::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}, method::ExpMethodHigham2005)
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/mb4U7/src/exp_noalloc.jl:88
 [12] exponential!(A::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/mb4U7/src/exp.jl:15
 [13] exp_generic(args::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16}; kwargs::@Kwargs{})
    @ ExponentialUtilities ./deprecated.jl:116
 [14] exp_generic(args::SMatrix{4, 4, ForwardDiff.Dual{Nothing, Float32, 1}, 16})
    @ ExponentialUtilities ./deprecated.jl:113
 [15] top-level scope
    @ REPL[343]:1
 [16] top-level scope
    @ none:1
Some type information was truncated. Use `show(err)` to see complete types.

```

does `exp_generic` work on static array/matrix? @ChrisRackauckas

---

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [July 2, 2026, 11:26am UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/9 "2026-07-02T11:26:31Z")

</div>

Thank you for your suggestions and link to another issues…but how to solve it for `cos`? our model also have `cos` and `sin` in other places…

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 2, 2026, 12:25pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/10 "2026-07-02T12:25:47Z")

</div>

Okay, so a few things here, @Xu_Shan you say

> [@Xu\_Shan](#):
>
> Yes I am using `exp` from `ExponentialUtilities.jl`…but seems to be broken too…
> 
> ```julia-auto
> julia> ForwardDiff.jacobian(ExponentialUtilities.exp, Float32.(A))
> ERROR: MethodError: no method matching exp!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{typeof(exp), Float32}, Float32, 9}})
> The function `exp!` exists, but no method is defined for this combination of argument types.
> 
> ```

but please note that `ExponentialUtilities.exp` is just `Base.exp`:

```julia-auto
julia> ExponentialUtilities.exp === Base.exp
true

```

Next, you mention that you are using StaticArrays, but on my machine at least, `ForwardDiff.exp` works totally fine on `SMatrix`:

```julia-auto
julia> using ForwardDiff, StaticArrays

julia> let A = rand(SMatrix{3, 3})
           ForwardDiff.jacobian(exp, A)
       end
9×9 SMatrix{9, 9, Float64, 81} with indices SOneTo(9)×SOneTo(9):
 1.76523 0.0467088 0.327616 0.806547 0.013858 0.100769 0.310702 0.00533908 0.0388228
 0.806547 1.59729 0.682086 0.245686 0.755658 0.207098 0.0946549 0.291097 0.0797884
 0.310702 0.0723638 1.67116 0.0946549 0.0224092 0.776938 0.0364675 0.0086335 0.299295
 0.0467088 0.0007806 0.00568669 1.59729 0.0438392 0.306742 0.0723638 0.00126505 0.00918839
 0.013858 0.0438392 0.0116795 0.755658 1.43989 0.639192 0.0224092 0.0677215 0.0188915
 0.00533908 0.00126505 0.0450359 0.291097 0.0677215 1.50936 0.0086335 0.00204291 0.0696659
 0.327616 0.00568669 0.0413231 0.682086 0.0116795 0.084947 1.67116 0.0450359 0.315479
 0.100769 0.306742 0.084947 0.207098 0.639192 0.174568 0.776938 1.50936 0.657123
 0.0388228 0.00918839 0.315479 0.0797884 0.0188915 0.657123 0.299295 0.0696659 1.58065

```

You also mention that you’re “forced to use `Float32`” at various points, and I suppose this is the actual problem. It seems the jacobian of `exp(::SMatrix{N, N, Float32})` isn’t very well behaved and can sometimes be `NaN`:

```julia-auto
julia> let A = rand(SMatrix{3, 3, Float32})
           ForwardDiff.jacobian(exp, A)
       end
9×9 SMatrix{9, 9, Float32, 81} with indices SOneTo(9)×SOneTo(9):
  0.610376 0.656253 0.547119 0.377831 -0.708216 0.067911 0.293016 0.0872581 -0.724079
 -0.097091 1.93593 0.744129 0.0534829 -0.0192576 0.107209 0.0449941 0.375067 -0.380432
 -0.0226828 0.280314 1.40075 0.0476514 -0.289407 0.378518 0.0397389 0.035158 -0.0227381
 -0.114842 0.178618 0.147041 1.93864 0.0291169 0.672426 0.277805 0.0731433 -0.725562
 -1.30146 0.801856 0.1931 0.458767 1.17613 0.894499 0.0250352 0.341649 -1.3225
 -0.274758 0.0704489 0.668607 0.374683 -0.0232118 1.93445 0.0328772 0.0287411 -0.0800439
 -0.11172 0.14564 0.119864 0.745587 -0.468221 0.161845 1.39804 0.670251 -0.111745
 -0.800691 0.671748 0.162845 0.109852 0.0235237 0.219857 0.375415 1.93717 -0.135987
 -0.709717 0.0463567 0.547772 0.0949469 -0.714231 0.73353 0.293344 0.282309 0.616263

julia> let A = rand(SMatrix{3, 3, Float32})
           ForwardDiff.jacobian(exp, A)
       end
9×9 SMatrix{9, 9, Float32, 81} with indices SOneTo(9)×SOneTo(9):
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN
 NaN NaN NaN NaN NaN NaN NaN NaN NaN

```

In this case, it seems that the easiest path forward would be to simply not use `Float32`, and eat the modest performance loss, but not suffer a massive loss of accuracy and reliability.

---

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [July 2, 2026, 12:43pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/11 "2026-07-02T12:43:59Z")

</div>

> [@Mason](#):
>
> It seems the jacobian of `exp(::SMatrix{N, N, Float32})` isn’t very well behaved and can sometimes be `NaN`

Yes…we also found it…we run the model optimization and ForwardDiff.jl (on optimized parameter) by 40 independent runs, and 8 or 5 of them broke and in each time it broke at different sites (also different numbers, i.e. “8 or 5” varies..)…

For the choice of `Float32`, if I remember it correctly, it is required by ForwardDiff.jl, right? for `StaticArrays` or static matrix, they are forced to use in the model to reduce the allocations and improve the model running time…otherwise each optimization would run at least 2 times slower…we did extensive efforts to make operations on static arrays to avoid excessive allocations…

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 2, 2026, 12:47pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/12 "2026-07-02T12:47:23Z")

</div>

> [@Xu\_Shan](#):
>
> For the choice of `Float32`, if I remember it correctly, it is required by [ForwardDiff.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ForwardDiff), right?

No? I literally show an example of it working on Float64. I suspect it was _only_ tested at 64 bit precision, and that Float32 was an afterthought if it was even thought about at all (hence the `NaN` issues).

> [@Xu\_Shan](#):
>
> for `StaticArrays` or static matrix, they are forced to use in the model to reduce the allocations and improve the model running time…otherwise each optimization would run at least 2 times slower…we did extensive efforts to make operations on static arrays to avoid excessive allocations…

Again, I literally show an example in my comment of it working fine on StaticArrays of Float64 numbers. It’s regular arrays that don’t have a `ForwardDiff` `exp` method.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [July 2, 2026, 2:19pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/13 "2026-07-02T14:19:02Z")

</div>

> [@Xu\_Shan](#):
>
> but how to solve it for `cos`

If you can get a working `exp` then one option is to use the Euler identity `exp(x*im) == cos(x) + sin(x)*im` (i.e., `cos` and `sin` are the real and imag parts of a complex exponential). If you don’t have a working `exp`:

- you can modify `expm1o` to compute `exp` directly by changing the lines `p = inv.(cumprod(ntuple(T, Val(9))))` to `p = inv.( (one(T), cumprod(ntuple(T, Val(9)))... ) )` and `expmA = (expmA * expmA) * A * (scale / 2) + expmA` to `expmA = expmA * expmA` (test this as I haven’t verified it).
- use `exp(x) == expm1o(x) * x + I`

If your `cos` and `sin` terms can actually be combined with adjacent terms into `(1 - cos(x)) / x` and `sin(x) / x` then you should use `expm1o(x*im)` to compute these expressions directly to avoid issues with singular `x`.

* * *

> [@Mason](#):
>
> It seems the jacobian of `exp(::SMatrix{N, N, Float32})` isn’t very well behaved

The frequent overflows in `exp(::SMatrix)` (especially the 3x3 case) are tracked in issue [#785](https://github.com/JuliaArrays/StaticArrays.jl/issues/785) which links to some other issues that are relevant as well.

The reason `exp(::SMatrix{3,3,Float32})` overflows so frequently is because the Pade approximant coefficients used are large (like `1f16`, if I recall correctly) and the 3x3 solve/inverse uses the `det` of the matrix which leads to overflow. If one downscaled the coefficients uniformly AND/OR one used `lu` to compute the Pade approximant the issue would largely disappear. I swear I had a comment on there talking about this but I can’t find it. The issue is present for 2x2 also but less of an issue because the `det` has fewer products. For 4x4+ it’s mostly fine because `lu` is used, so the large terms aren’t multiplied to overflow.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 3, 2026, 1:07pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/14 "2026-07-03T13:07:33Z")

</div>

A few updates after digging into this:

1. exp\_generic no longer exists as a name — the current API in ExponentialUtilities.jl is exponential!(A, method). The setindex!(::SMatrix, …) error comes from the default method selection: exponential!(A) with no method argument picks the in-place ExpMethodHigham2005, which mutates and therefore can’t work on an SMatrix.

2. The good news: exponential!(A, ExpMethodGeneric()) already works on static arrays today, including through ForwardDiff. ExpMethodGeneric has a non-mutating code path for immutable array types (Padé + scaling-and-squaring built only out of _, /, and + c_I), so:

```julia-auto
using ExponentialUtilities, StaticArrays, ForwardDiff
A = @SMatrix rand(3, 3)
exponential!(A, ExpMethodGeneric()) 
ForwardDiff.jacobian(x -> exponential!(x, ExpMethodGeneric()), A) 

```

Works

1. On the Float32 NaN issue @Mason found and @mikmoore diagnosed: confirmed, and ExpMethodGeneric avoids it. StaticArrays’ own exp (which is what exp(::SMatrix) calls, and which does work with dual numbers since its generic path is non-mutating) uses the Base-style Padé formulation with large integer coefficients (up to ~3.2e16), and those overflow Float32 on the way through the 3×3 cofactor-based solve. ExpMethodGeneric fixes both the mutation error and the Float32 instability at once.

2. I’ve opened a PR to ExponentialUtilities.jl making the default exponential!(A) automatically fall back to ExpMethodGeneric for immutable matrices, so plain exponential!(A) on an SMatrix will just work: [Fall back to ExpMethodGeneric for immutable matrices in exponential! - Pull Request #237 - SciML/ExponentialUtilities.jl - GitHub](https://github.com/SciML/ExponentialUtilities.jl/pull/237)

So with that, “just” using ExponentialUtilities.jl will do what you expected now.

---

<div class="post-metadata">

**Author:** ![Xu\_Shan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xu_shan/32/214900_2.png) [@Xu\_Shan](https://discourse.julialang.org/u/Xu_Shan)\
**Post date:** [July 3, 2026, 6:25pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/15 "2026-07-03T18:25:33Z")

</div>

> [@ChrisRackauckas](#):
>
> 1. [Fall back to ExpMethodGeneric for immutable matrices in exponential! - Pull Request #237 - SciML/ExponentialUtilities.jl - GitHub](https://github.com/SciML/ExponentialUtilities.jl/pull/237)

Thanks…but it seems that this one only sticks to Float64, right? it might use polynomial of Float64 internal?

```julia-auto
julia> A = SMatrix{3, 3, Float32}(rand(Float32, 3, 3))
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
 0.276386 0.954482 0.714192
 0.574739 0.3517 0.418869
 0.555151 0.3695 0.359785

julia> exponential!(A, ExpMethodGeneric()) 
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 2.11407 1.78511 1.49504
 1.13127 2.05394 1.03571
 1.08678 1.04209 1.95143

julia> ForwardDiff.jacobian(x -> exponential!(x, ExpMethodGeneric()), A) 
9×9 SMatrix{9, 9, Float64, 81} with indices SOneTo(9)×SOneTo(9):
 1.8236 0.838051 0.680035 0.523368 … 0.503557 0.158246 0.126138
 0.523368 1.82291 0.451652 0.101791 0.0980227 0.506292 0.0816883
 0.503557 0.441118 1.79091 0.0980227 0.0943939 0.0782926 0.501925
 0.838051 0.265264 0.211498 1.82291 0.441118 0.131362 0.104497
 0.164332 0.842565 0.137019 0.526213 0.0812944 0.443656 0.0674744
 0.158246 0.131362 0.835204 0.506292 … 0.0782926 0.0645242 0.440187
 0.680035 0.211498 0.16852 0.451652 1.79091 0.835204 0.677983
 0.130985 0.683771 0.109075 0.0848232 0.521681 1.7899 0.450526
 0.126138 0.104497 0.677983 0.0816883 0.501925 0.440187 1.75782

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 4, 2026, 8:31pm UTC](https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880/16 "2026-07-04T20:31:45Z")

</div>

Handled.