Post Reply 
(Hyper) Dual Numbers for automatic differentiation (CAS)
02-13-2024, 12:05 AM
Post: #1
(Hyper) Dual Numbers for automatic differentiation (CAS)
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
Find all posts by this user
Quote this message in a reply
02-13-2024, 03:51 AM (This post was last modified: 02-13-2024 10:23 PM by Albert Chan.)
Post: #2
RE: (Hyper) Dual Numbers for automatic differentiation (CAS)
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.
Find all posts by this user
Quote this message in a reply
02-13-2024, 07:18 AM (This post was last modified: 02-13-2024 07:39 AM by Thomas Klemm.)
Post: #3
RE: (Hyper) Dual Numbers for automatic differentiation (CAS)
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.
Find all posts by this user
Quote this message in a reply
02-13-2024, 08:57 AM (This post was last modified: 02-13-2024 02:12 PM by mcasl.)
Post: #4
RE: (Hyper) Dual Numbers for automatic differentiation (CAS)
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.
Find all posts by this user
Quote this message in a reply
02-13-2024, 02:22 PM
Post: #5
RE: (Hyper) Dual Numbers for automatic differentiation (CAS)
(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
Find all posts by this user
Quote this message in a reply
02-21-2024, 06:27 AM (This post was last modified: 02-21-2024 06:37 AM by Thomas Klemm.)
Post: #6
RE: (Hyper) Dual Numbers for automatic differentiation (CAS)
(02-13-2024 08:57 AM)mcasl Wrote:  The context is that I'm teaching automatic differentiation …

These articles could be useful:
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: