# GTPSA.jl: A Julia interface to the new GTPSA library for fast, high-order automatic differentiation

**URL:** <https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370>\
**Category:** Package Announcements\
**Tags:** package, announcement, differentiation, optimization, physics\
**Created:** [June 8, 2024, 8:45pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370 "2024-06-08T20:45:50Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![mattsignorelli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattsignorelli/32/221502_2.png) [@mattsignorelli](https://discourse.julialang.org/u/mattsignorelli)\
**Post date:** [June 8, 2024, 8:45pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/1 "2024-06-08T20:45:50Z")

</div>

[GTPSA.jl](https://github.com/bmad-sim/GTPSA.jl) is a new Julia package that provides a full-featured interface to the [Generalised Truncated Power Series Algebra library](https://mad.web.cern.ch/mad/releases/madng/html/mad_mod_diffalg.html#gtpsa) for computing Taylor expansions, or Truncated Power Series (TPSs) of real and complex multivariable functions to arbitrary orders.

Truncated Power Series Algebra (TPSA) performs forward-mode automatic differentation (AD) similar to the dual-number implementation as in ForwardDiff.jl. However, instead of nesting derivatives for higher orders, TPSA naturally extends to arbitrary orders by directly using the power series expansions. This, paired with a highly optimized monomial indexing function and storage for propagating the partial derivatives, makes GTPSA.jl significantly faster for 2nd-order calculations and above, and have similar performance at 1st-order. **In our [example](https://github.com/bmad-sim/GTPSA.jl/blob/main/benchmark/track.jl), GTPSA was x3.3 faster than ForwardDiff to 2nd order, and x18.5 faster to 3rd order.**

GTPSA provides several advantages over current Julia AD packages:

1. **Speed** : GTPSA.jl is significantly faster than ForwardDiff.jl for 2nd-order calculations and above, and has very similar performance at 1st-order
2. **Custom Orders in Individual Variables** : Other packages use a single maximum order for all variables. With GTPSA, the maximum order can be set differently for individual variables, as well as for a separate part of the monomial. For example, computing the Taylor expansion of f(x\_1,x\_2) to 1st order in x\_1 and 3rd order in x\_2 is possible
3. **Complex Numbers** : GTPSA.jl natively supports complex numbers and allows for mixing of complex and real truncated power series
4. **Distinction Between State Variables and Parameters** : Distinguishing between dependent variables and parameters in the solution of a differential equation expressed as a power series in the dependent variables/parameters can be advantageous in analysis

To read more about the inner-workings of the GTPSA C library, advanced users are referred to [this paper](https://inspirehep.net/files/286f2ab60e1e7c372cec485337ab5eb6), written by the creators of GTPSA.

As we’d like to prepare for a v1.0 release, we welcome any comments or criticisms of the package interface and documentation!

Because of the power of the GTPSA library, our long-term goal is to rewrite the library fully in Julia. If you are interested in this effort, feel free to reach out or add to the current [issue on Github](https://github.com/bmad-sim/GTPSA.jl/issues/96).

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 8, 2024, 9:20pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/2 "2024-06-08T21:20:01Z")

</div>

How does it compare to TaylorDiff.jl ? It seems like the same kind of ideas (leverage Faa di Bruno formula), but TaylorDiff does that at Julia’s compile time.

---

<div class="post-metadata">

**Author:** ![mattsignorelli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattsignorelli/32/221502_2.png) [@mattsignorelli](https://discourse.julialang.org/u/mattsignorelli)\
**Post date:** [June 8, 2024, 9:46pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/3 "2024-06-08T21:46:56Z")

</div>

For multivariable functions, TaylorDiff.jl only seems to give you the directional derivatives, whereas GTPSA.jl will give you all of the higher-order partial derivatives. Because of this, I haven’t been able to come up with an apples-to-apples benchmark to compare the two.

---

<div class="post-metadata">

**Author:** ![mattsignorelli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattsignorelli/32/221502_2.png) [@mattsignorelli](https://discourse.julialang.org/u/mattsignorelli)\
**Post date:** [June 8, 2024, 11:24pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/4 "2024-06-08T23:24:13Z")

</div>

Also, GTPSA does much more than just use the power series directly for calculating higher orders; the library plays a lot of tricks in both storing and propagating efficiently the many monomial coefficients, making it especially powerful for calculating derivatives w.r.t. to many variables to high orders.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 9, 2024, 1:11am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/5 "2024-06-09T01:11:59Z")

</div>

I opened two issues on your Github repo page. Some questions:

How does GTPSA.jl compare to Diffractor.jl?

How does GTPSA compare to the [ADOL-C](https://github.com/coin-or/ADOL-C) library? It also seems to have an in-development Julia wrapper package: [GitHub - TimSiebert1/ADOLC.jl](https://github.com/TimSiebert1/ADOLC.jl)

---

<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:** [June 9, 2024, 5:09am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/6 "2024-06-09T05:09:18Z")

</div>

Congrats on the release! Do you think it might make sense to add GPTSA.jl bindings to DifferentiationInterface.jl?  
I didn’t add any for TaylorDiff.jl because it does not support things like gradient or Hessian.

---

<div class="post-metadata">

**Author:** ![mattsignorelli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattsignorelli/32/221502_2.png) [@mattsignorelli](https://discourse.julialang.org/u/mattsignorelli)\
**Post date:** [June 9, 2024, 1:36pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/7 "2024-06-09T13:36:41Z")

</div>

> [@nsajko](#):
>
> I opened two issues on your Github repo page

Thank you! These are very helpful

> [@nsajko](#):
>
> How does [GTPSA.jl](https://juliahub.com/ui/Packages/General/GTPSA) compare to [Diffractor.jl](https://juliahub.com/ui/Packages/General/Diffractor)?

I just tried to add this as a benchmark in our [example](https://github.com/bmad-sim/GTPSA.jl/blob/main/benchmark/track.jl), however Diffractor.jl is returning an incorrect Jacobian. The benchmark uses a function m: \mathbb{R} ^{58} \rightarrow \mathbb{R}^6:

```julia-auto
julia> j = AD.jacobian(AD.ForwardDiffBackend(), m, zeros(58)) |> only
6×58 Matrix{Float64}:
 -0.699483 -3.08238 … 7.18128 4.7683 2.17832
  0.0406188 -1.25064 1.49797 1.36552 1.18237
  0.0 0.0 0.0 0.0 0.0
  0.0 0.0 0.0 0.0 0.0
  1.33161 -14.3074 -0.789128 -0.266338 0.0
  0.000133161 -0.00143074 … -7.89128e-5 -2.66338e-5 0.0

julia> j = AD.jacobian(DiffractorForwardBackend(), m, zeros(58)) |> only
6×58 Matrix{Float64}:
 -1.13036e-19 -2.52907e-21 … 0.170851 1.36948 2.17832
  3.21068e-22 -1.14018e-19 -0.133861 0.180479 1.18237
  0.0 0.0 0.0 0.0 0.0
  0.0 0.0 0.0 0.0 0.0
 -0.0115275 -0.113072 -0.17197 -0.12261 0.0
 -1.15275e-6 -1.13072e-5 … -1.7197e-5 -1.2261e-5 0.0

```

The Jacobian calculated with GTPSA.jl agrees exactly with ForwardDiff.jl. Therefore since Diffractor.jl is not calculating the correct Jacobian, it is probably not appropriate to add it as a benchmark yet. I will submit this as an issue to Diffractor.jl, though.

> [@nsajko](#):
>
> How does GTPSA compare to the [ADOL-C](https://github.com/coin-or/ADOL-C) library? It also seems to have an in-development Julia wrapper package: [GitHub - TimSiebert1/ADOLC.jl](https://github.com/TimSiebert1/ADOLC.jl)

I haven’t heard of this library, but I would be really interested to see how GTPSA compares. I’ll keep an eye on [GitHub - TimSiebert1/ADOLC.jl](https://github.com/TimSiebert1/ADOLC.jl), and once it is partially usable I can look into creating a benchmark.

---

<div class="post-metadata">

**Author:** ![mattsignorelli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattsignorelli/32/221502_2.png) [@mattsignorelli](https://discourse.julialang.org/u/mattsignorelli)\
**Post date:** [June 9, 2024, 1:44pm UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/8 "2024-06-09T13:44:55Z")

</div>

> [@gdalle](#):
>
> Congrats on the release! Do you think it might make sense to add [GPTSA.jl](https://juliahub.com/ui/Packages/General/GPTSA) bindings to [DifferentiationInterface.jl](https://juliahub.com/ui/Packages/General/DifferentiationInterface)?

Absolutely! We definitely want to use GTPSA.jl for AD with Julia’s many already-available optimizers, and this seems like an important first step. @nsajko already opened an [issue](https://github.com/bmad-sim/GTPSA.jl/issues/114) about this on the GTPSA Github, and I see in the [documentation](https://github.com/bmad-sim/GTPSA.jl/issues/114) some steps on how to do this. I’ll take a look and try to make sense of how to implement this interface, and will probably reach out on the issues with questions.

---

<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:** [June 10, 2024, 6:49am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/9 "2024-06-10T06:49:58Z")

</div>

> [@mattsignorelli](#):
>
> I just tried to add this as a benchmark in our [example](https://github.com/bmad-sim/GTPSA.jl/blob/main/benchmark/track.jl),

Note that for ForwardDiff.jl to be optimally efficient, you need to provide it with a precomputed `config` object. DifferentiationInterface.jl allows that easily with the backend-agnostic [“preparation” mechanism](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/operators/#Operators-Preparation).

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 10, 2024, 9:06am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/10 "2024-06-10T09:06:28Z")

</div>

> [@gdalle](#):
>
> I didn’t add any for [TaylorDiff.jl](https://juliahub.com/ui/Packages/General/TaylorDiff) because it does not support things like gradient or Hessian.

This is weird, it should. Maybe not as “interfaced” as other packages, but it should work.

---

<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:** [June 10, 2024, 9:15am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/11 "2024-06-10T09:15:16Z")

</div>

The only interface it has is `TaylorDiff.derivative`, see the [API reference](https://juliadiff.org/TaylorDiff.jl/stable/api/). Of course I could reconstruct everything else from there, I just didn’t take the time, and there might be smarter ways to do it from within TaylorDiff

---

<div class="post-metadata">

**Author:** ![oxinabox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oxinabox/32/206603_2.png) [@oxinabox](https://discourse.julialang.org/u/oxinabox)\
**Post date:** [June 10, 2024, 9:40am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/12 "2024-06-10T09:40:02Z")

</div>

One interesting case for a taylor mode AD is how it handles nesting.  
The goal is to avoid nested forward-mode and use taylor mode instead.  
But you can’t away do so easily at top level, because of encapsulation – some function is called in the middle of code you are ADing and that function happens to internally call AD.

In Diffractor, we have a flag (era-mode) that transforms nested forwards mode into taylor forward mode where it is encountered.  
I **think** it is possible to do this with operator overloading AD.

---

<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:** [June 10, 2024, 10:09am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/13 "2024-06-10T10:09:32Z")

</div>

> [@oxinabox](#):
>
> In Diffractor, we have a flag (era-mode) that transforms nested forwards mode into taylor forward mode where it is encountered.

In DifferentiationInterface.jl, we have the internal [`DifferentiationInterface.nested`](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/api/#DifferentiationInterface.nested-Tuple%7BADTypes.AbstractADType%7D) which lets the autodiff backend know it is inside another autodiff operator.

In Lux.jl there is something similar with [Nested Automatic Differentiation | Lux.jl Documentation](https://lux.csail.mit.edu/stable/manual/nested_autodiff#nested_autodiff)

---

<div class="post-metadata">

**Author:** ![oxinabox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oxinabox/32/206603_2.png) [@oxinabox](https://discourse.julialang.org/u/oxinabox)\
**Post date:** [June 11, 2024, 3:06am UTC](https://discourse.julialang.org/t/gtpsa-jl-a-julia-interface-to-the-new-gtpsa-library-for-fast-high-order-automatic-differentiation/115370/14 "2024-06-11T03:06:04Z")

</div>

Code in diffractor for reference:

> <https://github.com/JuliaDiff/Diffractor.jl/blob/main/src/stage1/forward.jl#L261-L263>