# Are computations with Dual numbers twice as slow as non Duals?

**URL:** <https://discourse.julialang.org/t/are-computations-with-dual-numbers-twice-as-slow-as-non-duals/92502>\
**Category:** Performance\
**Tags:** forwarddiff, dual\
**Created:** [January 4, 2023, 12:26pm UTC](https://discourse.julialang.org/t/are-computations-with-dual-numbers-twice-as-slow-as-non-duals/92502 "2023-01-04T12:26:43Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![nrontsis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nrontsis/32/38535_2.png) [@nrontsis](https://discourse.julialang.org/u/nrontsis)\
**Post date:** [January 4, 2023, 12:26pm UTC](https://discourse.julialang.org/t/are-computations-with-dual-numbers-twice-as-slow-as-non-duals/92502/1 "2023-01-04T12:26:43Z")

</div>

After reading [this](https://book.sciml.ai/notes/08-Forward-Mode_Automatic_Differentiation_(AD)_via_High_Dimensional_Algebras/) I got the impression that computations with dual numbers are typically expected to be twice as slow as non-duals. The reasoning was that both the normal and dual computations have to calculated at each operation, much like in complex arithmetic.

However I fail to observe this in the following minimal example. The example is motivated by the computation of parts of the gradient of a function `g(x::AbstractVector)`. Assuming forward mode AD, I can either split the function’s input into two vectors (one operating with dual numbers and the other not) or use dual numbers for everything. In doing the latter I don’t observe the nearly as much as `(2*5 + 2*20)/(2* 5 + 20)` = 1.667x slowdown that I would expect.

```julia
using ForwardDiff
using BenchmarkTools
using Base: Iterators

function f(x::AbstractVector{T}, y::AbstractVector) where {T}
    result = zero(T)
    for j in 1:1000 # 1000 "computations" to simulate a situation where f or g dominates allocation costs in override
        for value in x
            result += (0.9 * value)^j
        end
        for value in y
            result += (0.9 * value)^j
        end
    end
    return result
end

const x_values = ones(5)
const y_values = ones(20)

@btime ForwardDiff.gradient(x -> f(x, y_values), x_values)

const combined_values = vcat(x_values, y_values)

function g(combined::AbstractVector{T}) where {T}
    result = zero(T)
    for j in 1:1000 # 1000 "computations" to simulate a situation where f or g dominates allocation costs in override
        for value in combined
            result += (0.9 * value)^j
        end
    end
    return result
end

function override(x, y::T) where {T}
    result = T(x)
    for i in 1:length(y)
        result[i] = y[i]
    end
    return result
end

override(combined_values, x_values)
@btime ForwardDiff.gradient(x -> g(override(combined_values, x)), x_values)

```

Which gives:

```julia
  587.081 μs (4 allocations: 848 bytes)
  657.692 μs (5 allocations: 2.16 KiB)

```

---

<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:** [January 4, 2023, 1:51pm UTC](https://discourse.julialang.org/t/are-computations-with-dual-numbers-twice-as-slow-as-non-duals/92502/2 "2023-01-04T13:51:04Z")

</div>

> [@nrontsis](#):
>
> After reading [this](https://book.sciml.ai/notes/08-Forward-Mode_Automatic_Differentiation_(AD)_via_High_Dimensional_Algebras/) I got the impression that computations with dual numbers are typically expected to be twice as slow as non-duals. The reasoning was that both the normal and dual computations have to calculated at each operation, much like in complex arithmetic.

Normally you can SIMD the primal and dual parts because there’s a lot of related computations.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [January 4, 2023, 2:39pm UTC](https://discourse.julialang.org/t/are-computations-with-dual-numbers-twice-as-slow-as-non-duals/92502/3 "2023-01-04T14:39:08Z")

</div>

`^` the slow part is calculating the power of the values. The fast part is updating partials via multiplication.

```julia
if y == zero(y) || iszero(partials(x))
    new_partials = zero(partials(x))
else
    new_partials = partials(x) * y * ($f)(v, y - 1)
end

```

If you have non-zero partials with respect to `x` in `x ^ y`, it would calculate `^` twice – I assume for accuracy reasons instead of `x ^ (y-1)` and multiplying for `x ^ y`.

In your chunk mode gradient, the `partials(x)` will sometimes be 0, letting it skip the second `^`.  
Take a linear combination before looping and `^` (e.g. multiply by a dense square matrix) to ensure none of the partials are zero, and you’ll probably get closer to the expected 2x.