(Hyper) Dual Numbers for automatic differentiation (CAS)

[mcasl]

Dear all,
I'm new to HP Prime and I am taking my first steps in programming it.
In the context of CAS, I'm working on a way of automating the automatic differentiation of a function using dual numbers.
For example, let \[f\left(x,y\right)=\left(x^2+7y\right)*\left(7x^2+5y\right)\] . To find the derivative at x=3, y=2 we apply the function at \(x=3+\epsilon_1\) and \(y=2+\epsilon_2\), giving, after simplification:
\[f\left(3+\epsilon_1,2+\epsilon_2\right)=7\epsilon_1^4+84\epsilon_1^3+54\epsilon_​1^2 \epsilon_2+486 \epsilon_1^2+324\epsilon_1\epsilon_2+1404\epsilon_1+35\epsilon_2+626\epsilon_2+1​679 \]
Now, that expression can be further simplified considering that, by definition, \(\epsilon_1^2=\epsilon_2^2=\epsilon_1\epsilon_2=0\), yielding:
\[f\left(3+\epsilon_1,2+\epsilon_2\right)=1679+1404\epsilon_1+626\epsilon_2\]
The result shows that \(f\left(3,2\right)=1679\) and that \(\frac{\partial f}{\partial x}=1404\) and \(\frac{\partial f}{\partial y}=626\).
I'm currently using the CAS function taylor2d that I downloaded from the forums in order to discard the terms with \(\epsilon_i^n |_{ n>2}\) as well as \(\epsilon_i\epsilon_j\), but that approach does not cover multivariate cases with more than two variables.
I'm interested in knowing if there is an alternative way of simplifying the result, especially if it involves letting xcas know that the rules for the \(\epsilon_i\) are the ones previously mentioned.
Thanks for your patience and attention

[Albert Chan]

The whole point about automatic differentiation is that it is "automatic".
There is no need for Cas, taylor series ...

We store variable's value and its derivative in a list as {f, f'}. to represent (f + f' ε)
Then we just add rules for + - × ÷ ..., example:

(06-19-2022 10:50 AM)Albert Chan Wrote:  Apply f to dual number {g, g'} produce {f(g), (f(g))' = f'(g) * g'}

Above green link also have code for some useful rules, enough to solve your example

lua> a, b = D.new(7), D.new(5)
lua> function f(x,y) local x2=x*x; return (x2+a*y)*(a*x2+b*y) end
lua>
lua> unpack( f(D.new(3,1), D.new(2)) ) -- y as constant
1679      1404
lua> unpack( f(D.new(3), D.new(2,1)) ) -- x as constant
1679      626

Since f(real, real) = real, we can represent (ε1, ε2) as (ε1 + ε2*I)

lua> unpack( f(D.new(3,1), D.new(2,I)) ) -- derivative = ∂f/∂x + ∂f/∂y*I
1679      (1404+626*I)

There is no Cas. Everything is numerical.

[Thomas Klemm]

This is an implementation for the HP-42S: Automatic differentiation using dual numbers.
I used complex number as representation, so we get addition, subtraction and scalar multiplication for free.

But you could use the following matrix representation instead:

\(
\begin{bmatrix}
u & {u}' \\
0 & u \\
\end{bmatrix}
\)

With this you get also multiplication and division for free:

\(
\begin{bmatrix}
u & {u}' \\
0 & u \\
\end{bmatrix}
\cdot
\begin{bmatrix}
v & {v}' \\
0 & v \\
\end{bmatrix}
=
\begin{bmatrix}
uv & u{v}'+{u}'v \\
0 & uv \\
\end{bmatrix}
\)

\(
\begin{bmatrix}
u & {u}' \\
0 & u \\
\end{bmatrix}
\cdot
\begin{bmatrix}
v & {v}' \\
0 & v \\
\end{bmatrix}^{-1}
=
\begin{bmatrix}
\frac{u}{v} & \frac{{u}'v - u{v}'}{v^2} \\
0 & \frac{u}{v} \\
\end{bmatrix}
\)

In case of function application it really depends on how it is implemented.
The result should implement the chain-rule:

\(
f \left(
\begin{bmatrix}
u & {u}' \\
0 & u \\
\end{bmatrix}
\right) =
\begin{bmatrix}
f(u) & {f}'(u){u}' \\
0 & f(u) \\
\end{bmatrix}
\)

---

Example

Code:

from sympy import Matrix

x = Matrix([
    [3, 1],
    [0, 3]
])

y = Matrix([
    [2, 1],
    [0, 2]
])

(x**2 + 7*y)*(7*x**2 + 5*y)

\(
\begin{bmatrix}
1679 & 2030 \\
0 & 1679
\end{bmatrix}
\)

But here we use the same \(\varepsilon = \varepsilon_{1} = \varepsilon_{2}\).
So it might not be exactly what you want.

[mcasl]

Hi! Thank you for your interest.The context is that I'm teaching automatic differentiation and I propose my students simple exercises like the one I showed for them to solve using pen and paper. Thus the convenience of having a tool to assess the results.
I also tried the matrix approach, and I was in awe watching how xcas was capable of handling to perfection trigonometric expressions. In fact, it caused me great curiosity how xcas was smart enough to compute them.
I can think of solutions to the topic involving loops and substitution, or using symb2poly and further string processing, but I was intrigued in knowing if xcas was capable of learning the rules for the \( \epsilon_i \) and apply them in its intermediate computations, not just relying on user at the end of the process for discarding terms.

[KeithB]

(02-13-2024 08:57 AM)mcasl Wrote:  In fact, it caused me great curiosity how xcas was smart enough to compute them.

Because Parisse is smart.
https://www-fourier.ujf-grenoble.fr/~parisse/giac.html

[Thomas Klemm]

(02-13-2024 08:57 AM)mcasl Wrote:  The context is that I'm teaching automatic differentiation …

These articles could be useful:

• Understanding Automatic Differentiation in 30 lines of Python
  
• Building an LLM from scratch