Merge pull request #28 from HamletTanyavong/dev

Add second-order, forward-mode automatic differentiation

README.md

@@ -13,6 +13,6 @@
 
 ## About
 
-Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors.[^1] Mathematics.NET also supports first-order, forward and reverse-mode automatic differentiation.
+Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors.[^1] Mathematics.NET also supports first and second-order, forward and reverse-mode automatic differentiation.
 
 [^1]: Please visit the [documentation site](https://mathematics.hamlettanyavong.com) for detailed information.

docs/guide/autodiff/first-order-forward-mode.md

@@ -1,6 +1,6 @@
 # First-Order, Forward Mode Automatic Differentiation
 
-Support for first-order, forward-mode autodiff is provided by the `Dual<T>` class.
+Support for first-order, forward-mode autodiff is provided by the `Dual<T>` type.
 
 ## Dual Numbers
 
@@ -9,8 +9,13 @@ Forward-mode autodiff can be performed through the use of dual numbers which kee
 Dual<Real> x = new(1.23, 1.0);
 Dual<Real> y = new(2.34);
 ```
-The primal part represents the point at which we want to compute our derivative while the tangent part holds the information about our derivative. When we create a dual number with a tangent part, we specify a value that will be used as the seed, and it is important to know that the value of the derivative changes proportionally with this value.
+The primal part represents the point at which we want to compute our derivative while the tangent part holds the information about our derivative. When we create a dual number with a tangent part, we specify a value that will be used as the seed, and it is important to know that the value of the derivative changes proportionally with this value. We could also write, equivalently,
+```csharp
+using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
+var x = CreateVariable(1.23, 1.0);
+var y = CreateVariable(2.34);
+```
 Suppose we want to compute the partial derivative of the function
 $$
     f(x,y) = \frac{\sin(x + y)e^{-y}}{x^2+y^2+1}
@@ -19,32 +24,42 @@ at the points $ x=1.23 $ and $ y=2.34 $.
 with respect to $ x $. We must write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 // Add the following line to avoid having to type Dual<Real> every time we want to use a function.
 using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
-Dual<Real> x = new(1.23, 1.0);
-Dual<Real> y = new(2.34);
+Dual<Real> x = CreateVariable(1.23, 1.0);
+Dual<Real> y = CreateVariable(2.34);
 
 var result = Sin(x + y) * Exp(-y) / (x * x + y * y + 1);
 
 Console.WriteLine("∂f/∂x: {0}", result);
 ```
 Notice that we set the seed for the variable of interest, $ x $, to 1 while the seed for the variable we do not care about, $ y $, was set to 0. If we had set both to 1.0, then we would have computed the total derivative of the function instead. To compute the partial derivative of our function with respect to $ y $, we write
 ```csharp
-Dual<Real> x = new(1.23);
-Dual<Real> y = new(2.34, 1.0);
+Dual<Real> x = CreateVariable(1.23);
+Dual<Real> y = CreateVariable(2.34, 1.0);
 ```
 with the tangent part of the variable of interest set to 1 and the other to 0. Doing this will print the following to the console:
 ```
 ∂f/∂x: (-0.005009285670379789, -0.009425990481108835)
 ∂f/∂y: (-0.005009285670379789, -0.003024626925238263)
 ```
 
+### Total Derivative
+
+To get the total derivate of a function, set the seeds for each variable to `1.0`;
+```csharp
+Dual<Real> x = CreateVariable(1.23, 1.0);
+Dual<Real> y = CreateVariable(2.34, 1.0);
+// Repeat for each variable present
+```
+
 ## Dual Vectors
 
 We can create dual vectors to help us keep track of multiple dual numbers.
 ```csharp
-DualVector3<Real> x = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+DualVector3<Real> x = new(CreateVariable(1.23), CreateVariable(0.66), CreateVariable(2.34));
 ```
 We can use this to compute the vector-Jacobian product of the vector functions
 $$
@@ -57,9 +72,10 @@ $$
 with the vector $ v=(0.23,1.57,-1.71) $ at our points $ x_1=1.23 $, $ x_2=0.66 $, and $ x_3=2.34 $ by writing
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 using static Mathematics.NET.AutoDiff.Dual<Mathematics.NET.Core.Real>;
 
-DualVector3<Real> x = new((Dual<Real>)1.23, (Dual<Real>)0.66, (Dual<Real>)2.34);
+DualVector3<Real> x = new(CreateVariable(1.23), CreateVariable(0.66), CreateVariable(2.34));
 Vector3<Real> v = new(0.23, 1.57, -1.71);
 
 var result = DualVector3<Real>.VJP(

docs/guide/autodiff/first-order-reverse-mode.md

@@ -7,6 +7,7 @@ Support for first-order, reverse-mode automatic differentiation (autodiff) is pr
 Gradient tapes keep track of operations for autodiff; unlike forward-mode autodiff, tracking is required since gradients have to be calculated in reverse order. To begin using reverse-mode autodiff, we must create a gradient tape and assign it variables to track. These variables will be passed into and returned from methods that will compute the local gradients for us and record them on the tape.
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -33,6 +34,7 @@ $$
 at the point $ x=1.23 $. We can write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -107,6 +109,7 @@ Console.WriteLine(gradient[0]);
 The correct value for the derivative should be `3.525753368769319`. The complete code looks as follows:
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Real> tape = new();
 var x = tape.CreateVariable(1.23);
@@ -289,6 +292,7 @@ $$
 at the points $ z=1.23+i2.34 $ and $ w=-0.66+i0.23 $. We can write
 ```csharp
 using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
 
 GradientTape<Complex> tape = new();
 var z = tape.CreateVariable(new(1.23, 2.34));

docs/guide/autodiff/second-order-forward-mode.md

@@ -0,0 +1,40 @@
+# Second-Order, Forward Mode Automatic Differentiation
+
+Support for first-order, forward-mode autodiff is provided by the `HyperDual<T>` type. Because this type is used in a similar manner to `Dual<T>`, please refer to that section for help.
+
+### Second-Order Derivatives
+
+Suppose we wanted to find the second derivative of the complex function
+$$
+    f(z,w) = \sin(\tan{z}*\log{w})
+$$
+with respect to $ z $, at the points $ z=1.23+i0.66 $ and $ w=2.34-i0.25 $. We can do so by writing
+```csharp
+using Mathematics.NET.AutoDiff;
+using Mathematics.NET.Core;
+using static Mathematics.NET.AutoDiff.HyperDual<Mathematics.NET.Core.Complex>;
+
+var z = CreateVariable(new(1.23, 0.66), 1.0, 1.0);
+var w = CreateVariable(new(2.34, -0.25));
+
+var result = Sin(Tan(z) * Ln(w));
+
+Console.WriteLine(result.D3);
+```
+which will give us `(-6.158582087985498, 6.391603674636932)`. Notice that we now have to provide two seed values. Each one, very loosely speaking, "represents" a first-order derivative with respect to the variable in which it appears. Since there are two seeds present with the variable $ z $, it means we want to take its derivative twice. Similarly, we can write the following if we wanted the second derivative of our function with respect to $ w $:
+```csharp
+var z = CreateVariable(new(1.23, 0.66));
+var w = CreateVariable(new(2.34, -0.25), 1.0, 1.0);
+```
+This will give us `(0.30998196902728725, -0.11498565892578178)`. To get our mixed derivative, $ \partial f/\partial{z}\partial{w} $, we must indicate with we want one of each derivative
+```csharp
+var z = CreateVariable(new(1.23, 0.66), 1.0, 0.0);
+var w = CreateVariable(new(2.34, -0.25), 0.0, 1.0);
+```
+keeping in mind that the seeds must not occupy the same "slot." This, for example, will not give us the correct answer:
+```csharp
+// Incorrect, seeds must not occupy the same "slot"
+var z = CreateVariable(new(1.23, 0.66), 1.0, 0.0);
+var w = CreateVariable(new(2.34, -0.25), 1.0, 0.0);
+```
+This will print `(0.6670456012622978, 2.2955143408553718)` to the console.

docs/guide/autodiff/second-order-reverse-mode.md

@@ -1,6 +1,6 @@
 # Second-Order, Reverse Mode Automatic Differentiation
 
-Support for first-order, reverse-mode automatic differentiation (autodiff) is provided by the `HessianTape` class.
+Support for first-order, reverse-mode autodiff is provided by the `HessianTape` class.
 
 ## Hessian Tapes
 
@@ -26,7 +26,7 @@ $$
 $$
 Note that, in the future, we will not have to do this manually since there will be a method made specifically to compute Laplacians in spherical coordinates. For now, if we wanted to compute the Laplacian of the function
 $$
-    f(x,y,z) = \frac{\cos(x)}{(x+y)\sin(z)}
+    f(r,\theta,\phi) = \frac{\cos(r)}{(r+\theta)\sin(\phi)}
 $$
 we can write
 ```csharp
@@ -36,7 +36,7 @@ using Mathematics.NET.Core;
 HessianTape<Real> tape = new();
 var x = tape.CreateVariableVector(1.23, 0.66, 0.23);
 
-// f(x, y, z) = cos(x) / ((x + y) * sin(z))
+// f(r, θ, ϕ) = cos(r) / ((r + θ) * sin(ϕ))
 _ = tape.Divide(
         tape.Cos(x.X1),
         tape.Multiply(

docs/guide/toc.yml

@@ -16,3 +16,5 @@
     href: autodiff/first-order-forward-mode.md
   - name: Second-Order, Reverse-Mode Automatic Differentiation
     href: autodiff/second-order-reverse-mode.md
+  - name: Second-Order, Forward-Mode Automatic Differentiation
+    href: autodiff/second-order-forward-mode.md

docs/index.md

@@ -54,4 +54,4 @@
 
 ## About
 
-Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors. Mathematics.NET also supports forward and reverse-mode automatic differentiation.
+Mathematics.NET provides custom types for complex, real, and rational numbers as well as other mathematical objects such as vectors, matrices, and tensors. Mathematics.NET also supports first and second-order, forward and reverse-mode automatic differentiation.

src/Mathematics.NET/AutoDiff/AutoDiffExtensions.cs

@@ -33,19 +33,19 @@ namespace Mathematics.NET.AutoDiff;
 public static class AutoDiffExtensions
 {
     //
-    // Variable creation
+    // Variable vector creation
     //
 
-    /// <summary>Create a three-element vector from a seed vector of length three.</summary>
+    /// <summary>Create a dual vector from a seed vector.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A type that implements <see cref="ITape{T}"/></param>
-    /// <param name="x">A three-element vector of seed values</param>
+    /// <param name="x">A vector of seed values</param>
     /// <returns>A variable vector of length three</returns>
     public static VariableVector3<T> CreateVariableVector<T>(this ITape<T> tape, Vector3<T> x)
         where T : IComplex<T>, IDifferentiableFunctions<T>
         => new(tape.CreateVariable(x.X1), tape.CreateVariable(x.X2), tape.CreateVariable(x.X3));
 
-    /// <summary>Create a three-element vector from seed values.</summary>
+    /// <summary>Create a vector from seed values.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A type that implements <see cref="ITape{T}"/></param>
     /// <param name="x1Seed">The first seed value</param>
@@ -159,6 +159,25 @@ public static Vector3<T> Gradient<T>(this ITape<T> tape, Func<ITape<T>, Variable
         return new(gradients[0], gradients[1], gradients[2]);
     }
 
+    /// <summary>Compute the Hessian of a scaler function using reverse-mode automatic differentiation.</summary>
+    /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
+    /// <param name="tape">A gradient or Hessian tape</param>
+    /// <param name="f">A scalar function</param>
+    /// <param name="x">The point at which to compute the gradient</param>
+    /// <returns>The Hessian of the scalar function</returns>
+    public static Matrix3x3<T> Hessian<T>(this HessianTape<T> tape, Func<HessianTape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
+        where T : IComplex<T>, IDifferentiableFunctions<T>
+    {
+        _ = f(tape, x);
+
+        tape.ReverseAccumulation(out ReadOnlySpan2D<T> hessian);
+
+        return new(
+            hessian[0, 0], hessian[0, 1], hessian[0, 2],
+            hessian[1, 0], hessian[1, 1], hessian[1, 2],
+            hessian[2, 0], hessian[2, 1], hessian[2, 2]);
+    }
+
     /// <summary>Compute the Jacobian of a vector function using reverse-mode automatic differentiation: $ \nabla^\text{T}f_i(\textbf{x}) $ for $ i=\left\{1,2,3\right\} $.</summary>
     /// <typeparam name="T">A type that implements <see cref="IComplex{T}"/> and <see cref="IDifferentiableFunctions{T}"/></typeparam>
     /// <param name="tape">A gradient or Hessian tape</param>
@@ -233,7 +252,7 @@ public static Vector3<T> JVP<T>(
     /// <param name="f">A scalar function</param>
     /// <param name="x">The point at which to compute the Laplacian</param>
     /// <returns>The Laplacian of the scalar function</returns>
-    public static T Laplacian<T>(this HessianTape<T> tape, Func<ITape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
+    public static T Laplacian<T>(this HessianTape<T> tape, Func<HessianTape<T>, VariableVector3<T>, Variable<T>> f, VariableVector3<T> x)
         where T : IComplex<T>, IDifferentiableFunctions<T>
     {
         _ = f(tape, x);