Merge pull request #27 from HamletTanyavong/dev

Add second-order, reverse-mode automatic differentiation

README.md

@@ -13,7 +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 automatic differentiation.[^2]
+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.
 
 [^1]: Please visit the [documentation site](https://mathematics.hamlettanyavong.com) for detailed information.
-[^2]: So far, only first-order, reverse-mode automatic differentiation is supported; first-order, forward-mode, as well as second-order forward and reverse-modes are planned features.

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

@@ -2,7 +2,7 @@
 
 Support for first-order, reverse-mode automatic differentiation (autodiff) is provided by the `GradientTape` class.
 
-## Gradient tapes
+## Gradient Tapes
 
 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
@@ -98,11 +98,11 @@ graph BT
 ```
 We can then calculate the gradient of our function by using the `ReverseAccumulation` method.
 ```csharp
-tape.ReverseAccumulation(out var gradients);
+tape.ReverseAccumulation(out var gradient);
 ```
 Since this is a single variable equation, we can access the first element of `gradients` to get our result.
 ```csharp
-Console.WriteLine(gradients[0]);
+Console.WriteLine(gradient[0]);
 ```
 The correct value for the derivative should be `3.525753368769319`. The complete code looks as follows:
 ```csharp
@@ -119,12 +119,12 @@ var result = tape.Divide(
 // Optional: examine the nodes on the gradient tape
 tape.PrintNodes();
 
-tape.ReverseAccumulation(out var gradients);
+tape.ReverseAccumulation(out var gradient);
 
 // The value of the function at the point x = 1.23: 0.6675110878078776
 Console.WriteLine("Value: {0}", result);
 // The derivative of the function with respect to x at the point x = 1.23: 3.525753368769319
-Console.WriteLine("Derivative: {0}", gradients[0]);
+Console.WriteLine("Derivative: {0}", gradient[0]);
 ```
 
 ## Multivariable Equations
@@ -200,17 +200,17 @@ graph BT
     mul -- adj(w₆) = adj(w₇) ∂w₇/∂w₆ --> div
     div -- adj(f) = adj(w₇) = 1 (seed) --> function["f(x, y, z)"]
 ```
-As before, we can use the `ReverseAccumulation` to get our gradients
+As before, we can use `ReverseAccumulation` to get our gradients
 ```csharp
-tape.ReverseAccumulation(out var gradients);
+tape.ReverseAccumulation(out var gradient);
 ```
 and print them to the console with
 ```csharp
 using Mathematics.NET.LinearAlgebra;
 
 // code
 
-Console.WriteLine(gradients.ToDisplayString());
+Console.WriteLine(gradient.ToDisplayString());
 ```
 This will print the following to the console:
 ```
@@ -303,12 +303,12 @@ var result = tape.Cos(
 tape.PrintNodes(CancellationToken.None);
 Console.WriteLine();
 
-tape.ReverseAccumulation(out var gradients);
+tape.ReverseAccumulation(out var gradient);
 
 // The value of the function at the point z = 1.23 + i2.34 and w = -0.66 + i0.23
 Console.WriteLine("Value: {0}", result);
-// The gradients of the function: ∂f/∂z and ∂f/∂w, respectively
-Console.WriteLine("Gradients: {0}", gradients.ToDisplayString());
+// The gradient of the function: ∂f/∂z and ∂f/∂w, respectively
+Console.WriteLine("Gradient: {0}", gradient.ToDisplayString());
 ```
 which is almost the exact same code we would have written in the real case. (Note that some methods such as `Atan2` are not available for complex gradient tapes.) This should output the following to the console:
 ```
@@ -333,7 +333,7 @@ Node 5:
     Parents: [4, 5]
 
 Value: (27.784322505370138, 24.753716703326287)
-Gradients: [(126.28638563049401, -98.74954259806483),  (-38.801295827094066, -109.6878698782088)  ]
+Gradient: [(126.28638563049401, -98.74954259806483),  (-38.801295827094066, -109.6878698782088)  ]
 ```
 
 ## Custom Operations
@@ -348,9 +348,9 @@ var result = tape.CustomOperation(
     x => Real.Sin(x),  // The function
     x => Real.Cos(x)); // The derivative of the function
 
-tape.ReverseAccumulation(out var gradients);
+tape.ReverseAccumulation(out var gradient);
 Console.WriteLine("Value: {0}", result);
-Console.WriteLine("Gradient: {0}", gradients.ToDisplayString());
+Console.WriteLine("Gradient: {0}", gradient.ToDisplayString());
 ```
 For custom binary operations, we can write
 ```csharp

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

@@ -0,0 +1,58 @@
+# Second-Order, Reverse Mode Automatic Differentiation
+
+Support for first-order, reverse-mode automatic differentiation (autodiff) is provided by the `HessianTape` class.
+
+## Hessian Tapes
+
+The steps needed to perform second-order, reverse-mode autodiff is similar to the steps needed to perform the first-order case. This time, however, we have access to the following overloads and/or versions of `ReverseAccumulation`:
+```csharp
+HessianTape<Complex> tape = new();
+
+// Do some math...
+
+// Use when we are only interested in the gradient
+tape.ReverseAccumulation(out ReadOnlySpan<Complex> gradient);
+// Use when we are only interested in the Hessian
+tape.ReverseAccumulation(out ReadOnlySpan2D<Complex> hessian);
+// Use when we are interested in both the gradient and Hessian
+tape.ReverseAccumulation(out var gradient, out var hessian);
+```
+The last version may be useful for calculations such as finding the Laplacian of a scalar function in spherical coordinates which involves derivatives of first and second orders:
+$$
+\begin{align}
+    \nabla^2f(r,\theta,\phi)    &   =\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac{1}{r^2\sin{\theta}}\frac{\partial}{\partial\theta}\left(\sin{\theta}\frac{\partial f}{\partial\theta}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2f}{\partial\phi^2}   \\
+    &   =\frac{2}{r}\frac{\partial f}{\partial r}+\frac{\partial^2f}{\partial r^2}+\frac{1}{r^2\sin{\theta}}\left(\cos{\theta}\frac{\partial f}{\partial\theta}+\sin{\theta}\frac{\partial^2f}{\partial\theta^2}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2f}{\partial\phi^2}
+\end{align}
+$$
+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)}
+$$
+we can write
+```csharp
+using Mathematics.NET.AutoDiff;
+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))
+_ = tape.Divide(
+        tape.Cos(x.X1),
+        tape.Multiply(
+            tape.Add(x.X1, x.X2),
+            tape.Sin(x.X3)));
+
+tape.ReverseAccumulation(out var gradient, out var hessian);
+
+// Manual Laplacian computation
+var u = Real.One / (x.X1.Value * Real.Sin(x.X2.Value)); // 1 / (r * sin(θ))
+var laplacian = 2.0 * gradient[0] / x.X1.Value +
+                hessian[0, 0] +
+                u * Real.Cos(x.X2.Value) * gradient[1] / x.X1.Value +
+                hessian[1, 1] / (x.X1.Value * x.X1.Value) +
+                u * u * hessian[2, 2];
+
+Console.WriteLine(laplacian);
+```
+which should give us `48.80966092022821`.

docs/guide/toc.yml

@@ -14,3 +14,5 @@
     href: autodiff/first-order-reverse-mode.md
   - name: First-Order, Forward-Mode Automatic Differentiation
     href: autodiff/first-order-forward-mode.md
+  - name: Second-Order, Reverse-Mode Automatic Differentiation
+    href: autodiff/second-order-reverse-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 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.

src/Mathematics.NET.SourceGenerators/Mathematics.NET.SourceGenerators.csproj

@@ -15,7 +15,7 @@
       <PrivateAssets>all</PrivateAssets>
       <IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
     </PackageReference>
-    <PackageReference Include="Microsoft.CodeAnalysis.CSharp" Version="4.7.0" />
+    <PackageReference Include="Microsoft.CodeAnalysis.CSharp" Version="4.8.0" />
   </ItemGroup>
 
   <ItemGroup>