From ad47ce2dbb4f6d4a5728624f3bf3f529887c4d71 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 21 Nov 2023 04:55:25 -0600
Subject: [PATCH 01/29] Update README.md
---
README.md | 3 +--
1 file changed, 1 insertion(+), 2 deletions(-)
diff --git a/README.md b/README.md
index 6c737b4b..50f9be0e 100644
--- a/README.md
+++ b/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.
From c7e345307deac8db82e44e78395e05478816b9a4 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 21 Nov 2023 06:00:31 -0600
Subject: [PATCH 02/29] Create HessianNode.cs
---
src/Mathematics.NET/AutoDiff/HessianNode.cs | 89 +++++++++++++++++++++
1 file changed, 89 insertions(+)
create mode 100644 src/Mathematics.NET/AutoDiff/HessianNode.cs
diff --git a/src/Mathematics.NET/AutoDiff/HessianNode.cs b/src/Mathematics.NET/AutoDiff/HessianNode.cs
new file mode 100644
index 00000000..6c18d064
--- /dev/null
+++ b/src/Mathematics.NET/AutoDiff/HessianNode.cs
@@ -0,0 +1,89 @@
+//
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//
+
+using System.Runtime.InteropServices;
+
+namespace Mathematics.NET.AutoDiff;
+
+/// Represents a node on a Hessian tape
+/// A type that implements
+[StructLayout(LayoutKind.Sequential)]
+internal readonly record struct HessianNode
+ where T : IComplex
+{
+ /// The first derivative of the left component of the binary operation
+ public readonly T DX;
+ /// The second derivative of the left component of the binary operation
+ public readonly T DXX;
+ /// The derivative of the left or right component of the binary operation with respect to the left and right variables
+ public readonly T DXY;
+ /// The first derivative of the right component of the binary operation
+ public readonly T DY;
+ /// The second derivative of the right component of the binary operation
+ public readonly T DYY;
+
+ /// The parent index of the left node
+ public readonly int PX;
+ /// The parent index of the right node
+ public readonly int PY;
+
+ public HessianNode(int index)
+ {
+ DX = T.Zero;
+ DXX = T.Zero;
+ DXY = T.Zero;
+ DY = T.Zero;
+ DYY = T.Zero;
+
+ PX = index;
+ PY = index;
+ }
+
+ public HessianNode(T dfx, T dfxx, int px, int py)
+ {
+ DX = dfx;
+ DXX = dfxx;
+ DXY = T.Zero;
+ DY = T.Zero;
+ DYY = T.Zero;
+
+ PX = px;
+ PY = py;
+ }
+
+ public HessianNode(T dfx, T dfxx, T dfxy, T dfy, T dfyy, int px, int py)
+ {
+ DX = dfx;
+ DXX = dfxx;
+ DXY = dfxy;
+ DY = dfy;
+ DYY = dfyy;
+
+ PX = px;
+ PY = py;
+ }
+}
From fefc20de8a80eb3af3eddc83b67bbf8d424db1e7 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Tue, 21 Nov 2023 09:44:43 -0600
Subject: [PATCH 03/29] Create HessianTape.cs
---
src/Mathematics.NET/AutoDiff/HessianTape.cs | 629 ++++++++++++++++++++
1 file changed, 629 insertions(+)
create mode 100644 src/Mathematics.NET/AutoDiff/HessianTape.cs
diff --git a/src/Mathematics.NET/AutoDiff/HessianTape.cs b/src/Mathematics.NET/AutoDiff/HessianTape.cs
new file mode 100644
index 00000000..232ccf49
--- /dev/null
+++ b/src/Mathematics.NET/AutoDiff/HessianTape.cs
@@ -0,0 +1,629 @@
+//
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//
+
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Mathematics.NET.AutoDiff;
+
+/// Represents a Hessian tape
+/// A type that implements and
+public record class HessianTape
+ where T : IComplex, IDifferentiableFunctions
+{
+ private List> _nodes;
+ private int _variableCount;
+
+ public HessianTape()
+ {
+ _nodes = [];
+ }
+
+ /// Get the number of nodes on the gradient tape.
+ public int NodeCount => _nodes.Count;
+
+ /// Get the number of variables that are being tracked.
+ public int VariableCount => _variableCount;
+
+ //
+ // Methods
+ //
+
+ /// Create a variable for the gradient tape to track.
+ /// A seed value
+ /// A variable
+ public Variable CreateVariable(T seed)
+ {
+ _nodes.Add(new(_variableCount));
+ Variable variable = new(_variableCount++, seed);
+ return variable;
+ }
+
+ /// Print the nodes of the gradient tape to the console.
+ /// A cancellation token
+ /// The total number of nodes to print
+ public void PrintNodes(CancellationToken cancellationToken, int limit = 100)
+ {
+ const string tab = " ";
+
+ ReadOnlySpan> nodeSpan = CollectionsMarshal.AsSpan(_nodes);
+ HessianNode node;
+
+ int i = 0;
+ while (i < Math.Min(_variableCount, limit))
+ {
+ CheckForCancellation(cancellationToken);
+ node = nodeSpan[i];
+ Console.WriteLine($"Root Node {i}:");
+ Console.WriteLine($"{tab}Weights: [[{node.DX}, {node.DXX}, {node.DXY}],");
+ Console.WriteLine($"{tab} [{node.DY}, {node.DYY}, {node.DXY}]]");
+ Console.WriteLine($"{tab}Parents: [{node.PX}, {node.PY}]");
+ i++;
+ }
+
+ CheckForCancellation(cancellationToken);
+ Console.WriteLine();
+
+ while (i < Math.Min(nodeSpan.Length, limit))
+ {
+ CheckForCancellation(cancellationToken);
+ node = nodeSpan[i];
+ Console.WriteLine($"Node {i}:");
+ Console.WriteLine($"{tab}Weights: [[{node.DX}, {node.DXX}, {node.DXY}],");
+ Console.WriteLine($"{tab} [{node.DY}, {node.DYY}, {node.DXY}]]");
+ Console.WriteLine($"{tab}Parents: [{node.PX}, {node.PY}]");
+ i++;
+ }
+
+ static void CheckForCancellation(CancellationToken cancellationToken)
+ {
+ if (cancellationToken.IsCancellationRequested)
+ {
+ Console.WriteLine("Print node operation cancelled");
+ cancellationToken.ThrowIfCancellationRequested();
+ }
+ }
+ }
+
+ /// Perform reverse accumulation on the Hessian tape and output the resulting Hessian.
+ /// The gradient
+ /// The Hessian
+ [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
+ public void ReverseAccumulation(out ReadOnlySpan gradient, out ReadOnlySpan2D hessian)
+ => ReverseAccumulation(out gradient, out hessian, T.One);
+
+ // The following method uses the edge-pushing algorithm outlined by Gower and Mello: https://arxiv.org/pdf/2007.15040.pdf.
+ // TODO: use newer variations/versions of this algorithm since they are more performant
+
+ /// Perform reverse accumulation on the Hessian tape and output the resulting Hessian.
+ /// The gradient
+ /// The Hessian
+ /// A seed value
+ /// The Hessian tape does not have any tracked variables.
+ [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
+ public void ReverseAccumulation(out ReadOnlySpan gradient, out ReadOnlySpan2D hessian, T seed)
+ {
+ if (_variableCount == 0)
+ {
+ throw new Exception("Hessian tape contains no root nodes");
+ }
+
+ ReadOnlySpan> nodes = CollectionsMarshal.AsSpan(_nodes);
+ ref var start = ref MemoryMarshal.GetReference(nodes);
+ var length = nodes.Length;
+
+ Span gradientSpan = new T[length];
+ gradientSpan[length - 1] = seed;
+
+ Span2D hessianSpan = new T[length, length];
+
+ for (int i = length - 1; i >= _variableCount; i--)
+ {
+ var node = Unsafe.Add(ref start, i);
+ var gradientElement = gradientSpan[i];
+
+ EdgePush(hessianSpan, ref node, i);
+ Accumulate(hessianSpan, ref node, gradientElement);
+
+ gradientSpan[node.PX] += gradientElement * node.DX;
+ gradientSpan[node.PY] += gradientElement * node.DY;
+ }
+
+ gradient = gradientSpan[.._variableCount];
+ hessian = hessianSpan.Slice(0, 0, _variableCount, _variableCount);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
+ private static void EdgePush(Span2D weight, ref HessianNode node, int i)
+ {
+ for (int p = 0; p <= i; p++)
+ {
+ if (weight[i, p] == 0 || weight[p, i] == 0)
+ {
+ continue;
+ }
+ if (p != i)
+ {
+ if (node.PX != p)
+ {
+ var x = node.DX * weight[i, p];
+ weight[node.PX, p] += x;
+ weight[p, node.PX] += x;
+ }
+ else
+ {
+ weight[p, p] += 2 * node.DX * weight[i, p];
+ }
+
+ if (node.PY != p)
+ {
+ var x = node.DY * weight[i, p];
+ weight[node.PY, p] += x;
+ weight[p, node.PY] += x;
+ }
+ else
+ {
+ weight[p, p] += 2 * node.DY * weight[i, p];
+ }
+ }
+ else
+ {
+ var x = weight[i, i];
+ weight[node.PX, node.PX] += node.DX * node.DX * x;
+ weight[node.PX, node.PY] += node.DX * node.DY * x;
+ weight[node.PY, node.PX] += node.DY * node.DX * x;
+ weight[node.PY, node.PY] += node.DY * node.DY * x;
+ }
+ }
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
+ private static void Accumulate(Span2D weight, ref HessianNode node, T v)
+ {
+ weight[node.PX, node.PX] += v * node.DXX;
+ weight[node.PX, node.PY] += v * node.DXY;
+ weight[node.PY, node.PX] += v * node.DXY;
+ weight[node.PY, node.PY] += v * node.DYY;
+ }
+
+ //
+ // Basic operations
+ //
+
+ ///
+ public Variable Add(Variable x, Variable y)
+ {
+ _nodes.Add(new(T.One, T.Zero, T.Zero, T.One, T.Zero, x._index, y._index));
+ return new(_nodes.Count - 1, x.Value + y.Value);
+ }
+
+ ///
+ public Variable Add(T c, Variable x)
+ {
+ _nodes.Add(new(T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, c + x.Value);
+ }
+
+ ///
+ public Variable Add(Variable x, T c)
+ {
+ _nodes.Add(new(T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, x.Value + c);
+ }
+
+ ///
+ public Variable Divide(Variable x, Variable y)
+ {
+ var n = T.One / y.Value;
+ var dfxy = -n * n;
+ _nodes.Add(new(n, T.Zero, dfxy, x.Value * dfxy, -2.0 * n * x.Value * dfxy, x._index, y._index));
+ return new(_nodes.Count - 1, x.Value * n);
+ }
+
+ ///
+ public Variable Divide(T c, Variable x)
+ {
+ var n = T.One / x.Value;
+ var dfxy = -n * n;
+ _nodes.Add(new(c * dfxy, -2.0 * n * c * dfxy, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, c * n);
+ }
+
+ ///
+ public Variable Divide(Variable x, T c)
+ {
+ var n = T.One / c;
+ _nodes.Add(new(n, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, x.Value * n);
+ }
+
+ ///
+ public Variable Modulo(Variable x, Variable y)
+ {
+ _nodes.Add(new(T.One, T.Zero, T.Zero, x.Value * Real.Floor(x.Value / y.Value), T.Zero, x._index, y._index));
+ return new(_nodes.Count - 1, x.Value % y.Value);
+ }
+
+ ///
+ public Variable Modulo(Real c, Variable x)
+ {
+ _nodes.Add(new(c.Value * Real.Floor(c.Value / x.Value), T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, c % x.Value);
+ }
+
+ ///
+ public Variable Modulo(Variable x, Real c)
+ {
+ _nodes.Add(new(T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, x.Value % c);
+ }
+
+ ///
+ public Variable Multiply(Variable x, Variable y)
+ {
+ _nodes.Add(new(y.Value, T.Zero, T.One, x.Value, T.Zero, x._index, y._index));
+ return new(_nodes.Count - 1, x.Value * y.Value);
+ }
+
+ ///
+ public Variable Multiply(T c, Variable x)
+ {
+ _nodes.Add(new(c, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, c * x.Value);
+ }
+
+ ///
+ public Variable Multiply(Variable x, T c)
+ {
+ _nodes.Add(new(c, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, x.Value * c);
+ }
+
+ ///
+ public Variable Subtract(Variable x, Variable y)
+ {
+ _nodes.Add(new(T.One, T.Zero, T.Zero, -T.One, T.Zero, x._index, y._index));
+ return new(_nodes.Count - 1, x.Value - y.Value);
+ }
+
+ ///
+ public Variable Subtract(T c, Variable x)
+ {
+ _nodes.Add(new(-T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, c - x.Value);
+ }
+
+ ///
+ public Variable Subtract(Variable x, T c)
+ {
+ _nodes.Add(new(T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, x.Value - c);
+ }
+
+ //
+ // Other operations
+ //
+
+ ///
+ public Variable Negate(Variable x)
+ {
+ _nodes.Add(new(-T.One, T.Zero, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, -x.Value);
+ }
+
+ // Exponential functions
+
+ ///
+ public Variable Exp(Variable x)
+ {
+ var exp = T.Exp(x.Value);
+ _nodes.Add(new(exp, exp, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, exp);
+ }
+
+ ///
+ public Variable Exp2(Variable x)
+ {
+ var exp2 = T.Exp(x.Value);
+ var df = Real.Ln2 * exp2;
+ _nodes.Add(new(df, Real.Ln2 * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, exp2);
+ }
+
+ ///
+ public Variable Exp10(Variable x)
+ {
+ var exp10 = T.Exp(x.Value);
+ var df = Real.Ln10 * exp10;
+ _nodes.Add(new(df, Real.Ln10 * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, exp10);
+ }
+
+ // Hyperbolic functions
+
+ ///
+ public Variable Acosh(Variable x)
+ {
+ var u = x.Value - T.One;
+ var v = x.Value + T.One;
+ _nodes.Add(new(T.One / (T.Sqrt(u) * T.Sqrt(v)), -x.Value * T.Pow(u, -1.5) * T.Pow(v, -1.5), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Acosh(x.Value));
+ }
+
+ ///
+ public Variable Asinh(Variable x)
+ {
+ var u = T.One + x.Value * x.Value;
+ _nodes.Add(new(T.One / T.Sqrt(u), -x.Value * T.Pow(u, -1.5), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Asinh(x.Value));
+ }
+
+ ///
+ public Variable Atanh(Variable x)
+ {
+ var df = T.One / (T.One - x.Value * x.Value);
+ _nodes.Add(new(df, 2.0 * df * x.Value * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Atanh(x.Value));
+ }
+
+ ///
+ public Variable Cosh(Variable x)
+ {
+ var cosh = T.Cosh(x.Value);
+ _nodes.Add(new(T.Sinh(x.Value), cosh, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, cosh);
+ }
+
+ ///
+ public Variable Sinh(Variable x)
+ {
+ var sinh = T.Sinh(x.Value);
+ _nodes.Add(new(T.Cosh(x.Value), sinh, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, sinh);
+ }
+
+ ///
+ public Variable Tanh(Variable x)
+ {
+ var tanh = T.Tanh(x.Value);
+ var u = T.One / T.Cosh(x.Value);
+ var df = u * u;
+ _nodes.Add(new(df, -2.0 * df * tanh, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, tanh);
+ }
+
+ // Logarithmic functions
+
+ ///
+ public Variable Ln(Variable x)
+ {
+ var df = T.One / x.Value;
+ _nodes.Add(new(df, -df * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Ln(x.Value));
+ }
+
+ ///
+ public Variable Log(Variable x, Variable b)
+ {
+ var lnx = T.Ln(x.Value);
+ var lnb = T.Ln(b.Value);
+ var dfx = T.One / (lnb * x.Value);
+ var dfb = -lnx / (lnb * lnb * b.Value);
+ _nodes.Add(new(dfx, -dfx / x.Value, -dfx / (lnb * b.Value), dfb, -dfb * (2.0 / lnb + T.One) / b.Value, x._index, b._index));
+ return new(_nodes.Count - 1, T.Log(x.Value, b.Value));
+ }
+
+ ///
+ public Variable Log2(Variable x)
+ {
+ var u = T.One / x.Value;
+ var df = u / Real.Ln2;
+ _nodes.Add(new(df, -u * u / Real.Ln2, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Log2(x.Value));
+ }
+
+ ///
+ public Variable Log10(Variable x)
+ {
+ var u = T.One / x.Value;
+ var df = u / Real.Ln10;
+ _nodes.Add(new(df, -u * u / Real.Ln10, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Log10(x.Value));
+ }
+
+ // Power functions
+
+ ///
+ public Variable Pow(Variable x, Variable n)
+ {
+ var pow = T.Pow(x.Value, n.Value);
+ var lnx = T.Ln(x.Value);
+ var pownmo = T.Pow(x.Value, n.Value - T.One);
+ var dfn = lnx * pow;
+ _nodes.Add(new(
+ n.Value * pownmo,
+ (n.Value - T.One) * n.Value * T.Pow(x.Value, n.Value - 2.0),
+ (T.One + lnx * n.Value) * pownmo,
+ dfn,
+ lnx * dfn,
+ x._index,
+ n._index));
+ return new(_nodes.Count - 1, pow);
+ }
+
+ // Root functions
+
+ ///
+ public Variable Cbrt(Variable x)
+ {
+ var cbrt = T.Cbrt(x.Value);
+ var df = T.One / (3.0 * cbrt * cbrt);
+ _nodes.Add(new(df, -2.0 * df / (3.0 * x.Value), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, cbrt);
+ }
+
+ ///
+ public Variable Root(Variable x, Variable n)
+ {
+ var root = T.Root(x.Value, n.Value);
+ var lnx = T.Ln(x.Value);
+ var u = T.One / n.Value;
+ var v = T.One / x.Value;
+ var w = u * u;
+ var dfx = u * v * root;
+ var dfn = -lnx * root * w;
+ _nodes.Add(new(
+ dfx,
+ (n.Value - T.One) * v * root,
+ -(dfx * u + lnx * w),
+ dfn,
+ -(2.0 * u + lnx * w) * dfn,
+ x._index,
+ n._index));
+ return new(_nodes.Count - 1, root);
+ }
+
+ ///
+ public Variable Sqrt(Variable x)
+ {
+ var sqrt = T.Sqrt(x.Value);
+ var df = 0.5 / sqrt;
+ _nodes.Add(new(df, -0.5 / x.Value * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, sqrt);
+ }
+
+ // Trigonometric functions
+
+ ///
+ public Variable Cos(Variable x)
+ {
+ var cos = T.Cos(x.Value);
+ _nodes.Add(new(-T.Sin(x.Value), -cos, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, cos);
+ }
+
+ ///
+ public Variable Acos(Variable x)
+ {
+ var u = T.One - x.Value * x.Value;
+ _nodes.Add(new(-T.One / T.Sqrt(u), -x.Value * T.Pow(u, -1.5), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Acos(x.Value));
+ }
+
+ ///
+ public Variable Asin(Variable x)
+ {
+ var u = T.One - x.Value * x.Value;
+ _nodes.Add(new(T.One / T.Sqrt(u), x.Value * T.Pow(u, -1.5), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Asin(x.Value));
+ }
+
+ ///
+ public Variable Atan(Variable x)
+ {
+ var df = T.One / (T.One + x.Value * x.Value);
+ _nodes.Add(new(df, -2.0 * df * x.Value * df, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, T.Asin(x.Value));
+ }
+
+ ///
+ public Variable Atan2(Variable y, Variable x)
+ {
+ var u = y.Value * y.Value;
+ var v = x.Value * x.Value;
+ var a = T.One / (u + v);
+ var b = a * a;
+ var dfyy = -2.0 * x.Value * b * y.Value;
+ _nodes.Add(new(
+ x.Value * a,
+ dfyy,
+ (u - v) * b,
+ -y.Value * a,
+ -dfyy,
+ y._index,
+ x._index));
+ return new(_nodes.Count - 1, Real.Atan2(y.Value, x.Value));
+ }
+
+ ///
+ public Variable Sin(Variable x)
+ {
+ var sin = T.Sin(x.Value);
+ _nodes.Add(new(T.Cos(x.Value), -sin, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, sin);
+ }
+
+ ///
+ public Variable Tan(Variable x)
+ {
+ var tan = T.Tan(x.Value);
+ var sec = T.One / T.Cos(x.Value);
+ var df = sec * sec;
+ _nodes.Add(new(df, 2.0 * df * tan, x._index, _nodes.Count));
+ return new(_nodes.Count - 1, tan);
+ }
+
+ //
+ // Custom operations
+ //
+
+ /// Add a node to the Hessian tape using a custom unary operation.
+ /// A variable
+ /// A function
+ /// The derivative of the function
+ /// The second derivative of the function
+ /// A variable
+ public Variable CustomOperation(Variable x, Func f, Func dfx, Func dfxx)
+ {
+ _nodes.Add(new(dfx(x.Value), dfxx(x.Value), x._index, _nodes.Count));
+ return new(_nodes.Count - 1, f(x.Value));
+ }
+
+ /// Add a node to the Hessian tape using a custom binary operation.
+ /// The first variable
+ /// The second variable
+ /// A function
+ /// The first derivative of the function with respect to the first variable
+ /// The second derivative of the function with respect to the first variable
+ /// The second derivative of the function with respect to both variables
+ /// The first derivative of the function with respect to the second variable
+ /// The second derivative of the function with respect to the second variable
+ /// A variable
+ public Variable CustomOperation(
+ Variable x,
+ Variable y,
+ Func f,
+ Func dfx,
+ Func dfxx,
+ Func dfxy,
+ Func dfy,
+ Func dfyy)
+ {
+ _nodes.Add(new(dfx(x.Value, y.Value), dfxx(x.Value, y.Value), dfxy(x.Value, y.Value), dfy(x.Value, y.Value), dfyy(x.Value, y.Value), x._index, y._index));
+ return new(_nodes.Count - 1, f(x.Value, y.Value));
+ }
+}
From fbfe94f6581ab1dd8f8dc79c1a8b42854ad1b500 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:14:56 -0600
Subject: [PATCH 04/29] Update GradientTape.cs
- Update documentation comments
---
src/Mathematics.NET/AutoDiff/GradientTape.cs | 9 +++++----
1 file changed, 5 insertions(+), 4 deletions(-)
diff --git a/src/Mathematics.NET/AutoDiff/GradientTape.cs b/src/Mathematics.NET/AutoDiff/GradientTape.cs
index 27dd5b60..eabff864 100644
--- a/src/Mathematics.NET/AutoDiff/GradientTape.cs
+++ b/src/Mathematics.NET/AutoDiff/GradientTape.cs
@@ -142,15 +142,16 @@ static void CheckForCancellation(CancellationToken cancellationToken)
}
}
- /// Perform reverse accumulation on the gradient tape and output the resulting gradients.
- /// The gradients
+ /// Perform reverse accumulation on the gradient tape and output the resulting gradient.
+ /// The gradient
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public void ReverseAccumulation(out ReadOnlySpan gradients)
=> ReverseAccumulation(out gradients, T.One);
- /// Perform reverse accumulation on the gradient tape and output the resulting gradients.
- /// The gradients
+ /// Perform reverse accumulation on the gradient tape and output the resulting gradient.
+ /// The gradient
/// A seed value
+ /// The gradient tape does not have any tracked variables.
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public void ReverseAccumulation(out ReadOnlySpan gradients, T seed)
{
From f2a3fd0da6ec30631fc17455f621d1200a7f5ed8 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:15:25 -0600
Subject: [PATCH 05/29] Add method for formatting 2D read-only spans
---
.../LinearAlgebra/LinAlgExtensions.cs | 14 ++++++++++++++
1 file changed, 14 insertions(+)
diff --git a/src/Mathematics.NET/LinearAlgebra/LinAlgExtensions.cs b/src/Mathematics.NET/LinearAlgebra/LinAlgExtensions.cs
index acc8e214..c0686c9e 100644
--- a/src/Mathematics.NET/LinearAlgebra/LinAlgExtensions.cs
+++ b/src/Mathematics.NET/LinearAlgebra/LinAlgExtensions.cs
@@ -93,6 +93,20 @@ public static string ToDisplayString(this Span span, string? format = null
return string.Format(provider, builder.ToString());
}
+ /// Get the string representation of this
+ /// A type that implements
+ /// A 2D read-only span to format
+ /// The format to use
+ /// The provider to use to format the value
+ /// A string representation of this object
+ public static string ToDisplayString(this ReadOnlySpan2D readOnlySpan2D, string? format = null, IFormatProvider? provider = null)
+ where T : IComplex
+ {
+ Span2D span = new T[readOnlySpan2D.Width, readOnlySpan2D.Height];
+ readOnlySpan2D.CopyTo(span);
+ return span.ToDisplayString(format, provider);
+ }
+
/// Get the string representation of this object
/// A type that implements
/// The span to format
From d89d21c560605e20bcc96a4ae7c7ecfec50a0d32 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:15:55 -0600
Subject: [PATCH 06/29] Make helpers private
---
.../AutoDiff/DualVector3OfRealTests.cs | 8 ++++----
.../AutoDiff/GradientTapeExtensionsOfRealTests.cs | 8 ++++----
2 files changed, 8 insertions(+), 8 deletions(-)
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs
index 04ea9378..7ab90b59 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/DualVector3OfRealTests.cs
@@ -133,15 +133,15 @@ public void VJP_VectorAndR3VectorFunction_ReturnsVJP(double vx, double vy, doubl
// Helpers
//
- public static Dual F(DualVector3 x)
+ private static Dual F(DualVector3 x)
=> Cos(x.X1) / ((x.X1 + x.X2) * Sin(x.X3));
- public static Dual FX(DualVector3 x)
+ private static Dual FX(DualVector3 x)
=> Sin(x.X1) * (Cos(x.X2) + Sqrt(x.X3));
- public static Dual FY(DualVector3 x)
+ private static Dual FY(DualVector3 x)
=> Sqrt(x.X1 + x.X2 + x.X3);
- public static Dual FZ(DualVector3 x)
+ private static Dual FZ(DualVector3 x)
=> Sinh(Exp(x.X1) * x.X2 / x.X3);
}
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs
index 1e2c05de..4dfdfd4e 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeExtensionsOfRealTests.cs
@@ -140,7 +140,7 @@ public void VJP_VectorAndR3VectorFunction_ReturnsVJP(double vx, double vy, doubl
//
// f(x, y, z) = Cos(x) / ((x + y) * Sin(z))
- public static Variable F(GradientTape tape, VariableVector3 x)
+ private static Variable F(GradientTape tape, VariableVector3 x)
{
return tape.Divide(
tape.Cos(x.X1),
@@ -150,7 +150,7 @@ public static Variable F(GradientTape tape, VariableVector3 x)
}
// f(x, y, z) = Sin(x) * (Cos(y) + Sqrt(z))
- public static Variable FX(GradientTape tape, VariableVector3 x)
+ private static Variable FX(GradientTape tape, VariableVector3 x)
{
return tape.Multiply(
tape.Sin(x.X1),
@@ -160,7 +160,7 @@ public static Variable FX(GradientTape tape, VariableVector3 x
}
// f(x, y, z) = Sqrt(x + y + z)
- public static Variable FY(GradientTape tape, VariableVector3 x)
+ private static Variable FY(GradientTape tape, VariableVector3 x)
{
return tape.Sqrt(
tape.Add(
@@ -171,7 +171,7 @@ public static Variable FY(GradientTape tape, VariableVector3 x
}
// f(x, y, z) = Sinh(Exp(x) * y / z)
- public static Variable FZ(GradientTape tape, VariableVector3 x)
+ private static Variable FZ(GradientTape tape, VariableVector3 x)
{
return tape.Sinh(
tape.Multiply(
From 6e8b4d177bfd2a9fafbde2052a8ccb1909c5ff82 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:15:59 -0600
Subject: [PATCH 07/29] Update Mathematics.NET.SourceGenerators.csproj
---
.../Mathematics.NET.SourceGenerators.csproj | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/Mathematics.NET.SourceGenerators/Mathematics.NET.SourceGenerators.csproj b/src/Mathematics.NET.SourceGenerators/Mathematics.NET.SourceGenerators.csproj
index ea490301..31e273d9 100644
--- a/src/Mathematics.NET.SourceGenerators/Mathematics.NET.SourceGenerators.csproj
+++ b/src/Mathematics.NET.SourceGenerators/Mathematics.NET.SourceGenerators.csproj
@@ -15,7 +15,7 @@
all
runtime; build; native; contentfiles; analyzers; buildtransitive
-
+
From 5d0d3669c9665aa06ae674ba63671cedbe21df87 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:16:05 -0600
Subject: [PATCH 08/29] Update first-order-reverse-mode.md
---
docs/guide/autodiff/first-order-reverse-mode.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/docs/guide/autodiff/first-order-reverse-mode.md b/docs/guide/autodiff/first-order-reverse-mode.md
index bb5c0c2b..5b886561 100644
--- a/docs/guide/autodiff/first-order-reverse-mode.md
+++ b/docs/guide/autodiff/first-order-reverse-mode.md
@@ -200,7 +200,7 @@ 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);
```
From c5927deb6956ee8b2efb6e7840a18e05e57993fa Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:44:19 -0600
Subject: [PATCH 09/29] Update GradientTapeOfRealTests.cs
- Rename methods and some parameters
---
.../AutoDiff/GradientTapeOfRealTests.cs | 28 +++++++++----------
1 file changed, 14 insertions(+), 14 deletions(-)
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeOfRealTests.cs
index 5f91f194..06d484de 100644
--- a/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeOfRealTests.cs
+++ b/tests/Mathematics.NET.Tests/AutoDiff/GradientTapeOfRealTests.cs
@@ -69,7 +69,7 @@ public void Add_TwoVariables_ReturnsGradients(double left, double right, double
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Add, left, right);
+ var actual = ComputeGradient(_tape.Add, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-16);
}
@@ -133,7 +133,7 @@ public void Atan2_TwoVariables_ReturnsGradients(double left, double right, doubl
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Atan2, left, right);
+ var actual = ComputeGradient(_tape.Atan2, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -233,7 +233,7 @@ public void Divide_TwoVariables_ReturnsGradients(double left, double right, doub
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Divide, left, right);
+ var actual = ComputeGradient(_tape.Divide, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -306,7 +306,7 @@ public void Log_TwoVariables_ReturnsGradients(double left, double right, double
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Log, left, right);
+ var actual = ComputeGradient(_tape.Log, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -335,7 +335,7 @@ public void Modulo_TwoVariables_ReturnsGradients(double left, double right, doub
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Modulo, left, right);
+ var actual = ComputeGradient(_tape.Modulo, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -372,7 +372,7 @@ public void Multiply_TwoVariables_ReturnsGradients(double left, double right, do
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Multiply, left, right);
+ var actual = ComputeGradient(_tape.Multiply, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-16);
}
@@ -422,7 +422,7 @@ public void Pow_TwoVariables_ReturnsGradients(double left, double right, double
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Pow, left, right);
+ var actual = ComputeGradient(_tape.Pow, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -433,7 +433,7 @@ public void Root_TwoVariables_ReturnsGradients(double left, double right, double
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Root, left, right);
+ var actual = ComputeGradient(_tape.Root, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-15);
}
@@ -471,7 +471,7 @@ public void Subtract_TwoVariables_ReturnsGradients(double left, double right, do
{
Real[] expected = [expectedLeft, expectedRight];
- var actual = ComputeGradients(_tape.Subtract, left, right);
+ var actual = ComputeGradient(_tape.Subtract, left, right);
Assert.AreApproximatelyEqual(expected, actual, 1e-16);
}
@@ -528,16 +528,16 @@ private Real ComputeGradient(Func, Variable> function, Real
{
var x = _tape.CreateVariable(input);
_ = function(x);
- _tape.ReverseAccumulation(out var gradients);
- return gradients[0];
+ _tape.ReverseAccumulation(out var gradient);
+ return gradient[0];
}
- private Real[] ComputeGradients(Func, Variable, Variable> function, Real left, Real right)
+ private Real[] ComputeGradient(Func, Variable, Variable> function, Real left, Real right)
{
var x = _tape.CreateVariable(left);
var y = _tape.CreateVariable(right);
_ = function(x, y);
- _tape.ReverseAccumulation(out var gradients);
- return gradients.ToArray();
+ _tape.ReverseAccumulation(out var gradient);
+ return gradient.ToArray();
}
}
From 8182ef428f9da562dd54d42ac2a1eddc423be26a Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 04:52:19 -0600
Subject: [PATCH 10/29] Create HessianTapeOfRealTests.cs
---
.../AutoDiff/HessianTapeOfRealTests.cs | 557 ++++++++++++++++++
1 file changed, 557 insertions(+)
create mode 100644 tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
diff --git a/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs b/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
new file mode 100644
index 00000000..5f29faaa
--- /dev/null
+++ b/tests/Mathematics.NET.Tests/AutoDiff/HessianTapeOfRealTests.cs
@@ -0,0 +1,557 @@
+//
+// Mathematics.NET
+// https://github.com/HamletTanyavong/Mathematics.NET
+//
+// MIT License
+//
+// Copyright (c) 2023 Hamlet Tanyavong
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//
+
+using Mathematics.NET.AutoDiff;
+
+namespace Mathematics.NET.Tests.AutoDiff;
+
+[TestClass]
+[TestCategory("AutoDiff"), TestCategory("Hessian Tape")]
+public sealed class HessianTapeOfRealTests
+{
+ private HessianTape _tape;
+
+ public HessianTapeOfRealTests()
+ {
+ _tape = new();
+ }
+
+ //
+ // Tests
+ //
+
+ [TestMethod]
+ [DataRow(0.123, -0.125845035324435)]
+ public void Acos_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Acos, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -3.348544407755665)]
+ public void Acosh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Acosh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0, 0, 0)]
+ public void Add_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Add, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Add_ConstantAndVariable_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(right);
+ _ = _tape.Add(left, x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Add_VariableAndConstant_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(left);
+ _ = _tape.Add(x, right);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(0.123, 0.125845035324435)]
+ public void Asin_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Asin, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.3087751219667227)]
+ public void Asinh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Asinh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.3895692725912341)]
+ public void Atan_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Atan, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, -0.1178645019844717, -0.081137805227897, 0.1178645019844717)]
+ public void Atan2_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Atan2, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-14);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 9.35125088756106)]
+ public void Atanh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Atanh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.1573785841306649)]
+ public void Cbrt_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Cbrt, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.3342377271245026)]
+ public void Cos_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Cos, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 1.856761056985266)]
+ public void Cosh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Cosh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(2)]
+ [DataRow(3)]
+ [DataRow(4)]
+ public void CreateVariable_Multiple_TracksCorrectNumberOfVariables(int amount)
+ {
+ for (int i = 0; i < amount; i++)
+ {
+ _ = _tape.CreateVariable(0);
+ }
+
+ var actual = _tape.VariableCount;
+
+ Assert.AreEqual(amount, actual);
+ }
+
+ [TestMethod]
+ [DataRow(1.23)]
+ public void CreateVariable_WithSeedValue_ReturnsVariable(double value)
+ {
+ var actual = _tape.CreateVariable(value).Value;
+
+ Assert.AreEqual(value, actual);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, -0.1178645019844717, -0.081137805227897, 0.1178645019844717)]
+ public void CustomOperation_Binary_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ var y = _tape.CreateVariable(left);
+ var x = _tape.CreateVariable(right);
+
+ var u = y.Value * y.Value;
+ var v = x.Value * x.Value;
+ var a = Real.One / (u + v);
+ var b = a * a;
+ var dfyy = -2.0 * x.Value * b * y.Value;
+
+ _ = _tape.CustomOperation(
+ y,
+ x,
+ Real.Atan2,
+ (y, x) => x.Value * a,
+ (y, x) => dfyy,
+ (y, x) => (u - v) * b,
+ (y, x) => -y.Value * a,
+ (y, x) => -dfyy);
+
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian.ToArray();
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-14);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.942488801931697)]
+ public void CustomOperation_Unary_ReturnsHessian(double input, double expected)
+ {
+ var x = _tape.CreateVariable(input);
+ _ = _tape.CustomOperation(x, Real.Sin, Real.Cos, x => -Real.Sin(x));
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0, -0.1826283877565929, 0.1919939461030848)]
+ public void Divide_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Divide, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0.1919939461030848)]
+ public void Divide_ConstantAndVariable_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(right);
+ _ = _tape.Divide(left, x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Divide_VariableAndConstant_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(left);
+ _ = _tape.Divide(x, right);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 3.421229536289673)]
+ public void Exp_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Exp, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 1.126984172374114)]
+ public void Exp2_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Exp2, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 90.0391481211886)]
+ public void Exp10_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Exp10, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.6609822195782934)]
+ public void Ln_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Ln, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, -0.7774880868134957, -0.4807142135095495, 0.1753670976636016)]
+ public void Log_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Log, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.953595770301384)]
+ public void Log2_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Log2, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.2870609305990163)]
+ public void Log10_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Log10, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0, 0, 0)]
+ public void Modulo_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Modulo, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Modulo_ConstantAndVariable_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(right);
+ _ = _tape.Modulo(left, x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Modulo_VariableAndConstant_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(left);
+ _ = _tape.Modulo(x, right);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0, 1, 0)]
+ public void Multiply_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Multiply, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Multiply_ConstantAndVariable_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(right);
+ _ = _tape.Multiply(left, x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Multiply_VariableAndConstant_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(left);
+ _ = _tape.Multiply(x, right);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 0)]
+ public void Negate_Variable_ReturnsNegation(double input, double expected)
+ {
+ var x = _tape.CreateVariable(input);
+
+ var actual = ComputeHessian(_tape.Negate, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 3.364251027050465, 1.958969363947686, 0.06956296832768787)]
+ public void Pow_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Pow, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, -0.1767192454212087, -0.1765629860861052, 0.03686389570982799)]
+ public void Root_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Root, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.942488801931697)]
+ public void Sin_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Sin, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 1.564468479304407)]
+ public void Sinh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Sinh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.1832661859080147)]
+ public void Sqrt_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Sqrt, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0, 0, 0)]
+ public void Subtract_TwoVariables_ReturnsHessian(double left, double right, double expectedXX, double expectedXY, double expectedYY)
+ {
+ Real[,] expected = new Real[2, 2] { { expectedXX, expectedXY }, { expectedXY, expectedYY } };
+
+ var actual = ComputeHessian(_tape.Subtract, left, right);
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Subtract_ConstantAndVariable_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(right);
+ _ = _tape.Subtract(left, x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 2.34, 0)]
+ public void Subtract_VariableAndConstant_ReturnsHessian(double left, double right, double expected)
+ {
+ var x = _tape.CreateVariable(left);
+ _ = _tape.Subtract(x, right);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+
+ var actual = hessian[0, 0];
+
+ Assert.AreApproximatelyEqual(expected, actual, Real.Zero);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, 50.4823759141874)]
+ public void Tan_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Tan, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ [TestMethod]
+ [DataRow(1.23, -0.4887972531670078)]
+ public void Tanh_Variable_ReturnsHessian(double input, double expected)
+ {
+ var actual = ComputeHessian(_tape.Tanh, input);
+
+ Assert.AreApproximatelyEqual(expected, actual, 1e-15);
+ }
+
+ //
+ // Helpers
+ //
+
+ private Real ComputeHessian(Func, Variable> function, Real input)
+ {
+ var x = _tape.CreateVariable(input);
+ _ = function(x);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+ return hessian[0, 0];
+ }
+
+ private Real[,] ComputeHessian(Func, Variable, Variable> function, Real left, Real right)
+ {
+ var x = _tape.CreateVariable(left);
+ var y = _tape.CreateVariable(right);
+ _ = function(x, y);
+ _tape.ReverseAccumulation(out var _, out var hessian);
+ return hessian.ToArray();
+ }
+}
From 6330f3bf46190d11b9ab552270edd595fbe889f8 Mon Sep 17 00:00:00 2001
From: Hamlet Tanyavong <34531738+HamletTanyavong@users.noreply.github.com>
Date: Wed, 22 Nov 2023 23:53:41 -0600
Subject: [PATCH 11/29] Update HessianTape.cs
- Fix Exp2, Exp10, and Root methods
- Use Real.One instead of T.One in Atan2
---
src/Mathematics.NET/AutoDiff/HessianTape.cs | 16 ++++++++--------
1 file changed, 8 insertions(+), 8 deletions(-)
diff --git a/src/Mathematics.NET/AutoDiff/HessianTape.cs b/src/Mathematics.NET/AutoDiff/HessianTape.cs
index 232ccf49..2e77dd1e 100644
--- a/src/Mathematics.NET/AutoDiff/HessianTape.cs
+++ b/src/Mathematics.NET/AutoDiff/HessianTape.cs
@@ -348,7 +348,7 @@ public Variable Exp(Variable x)
///
public Variable Exp2(Variable x)
{
- var exp2 = T.Exp(x.Value);
+ var exp2 = T.Exp2(x.Value);
var df = Real.Ln2 * exp2;
_nodes.Add(new(df, Real.Ln2 * df, x._index, _nodes.Count));
return new(_nodes.Count - 1, exp2);
@@ -357,7 +357,7 @@ public Variable Exp2(Variable x)
///