Hessian and Jacobian

src/CMakeLists.txt

@@ -41,7 +41,9 @@ set(public_headers
     XAD/Complex.hpp
     XAD/Exceptions.hpp
     XAD/Expression.hpp
+    XAD/Hessian.hpp	
     XAD/Interface.hpp	
+    XAD/Jacobian.hpp
     XAD/Literals.hpp
     XAD/Macros.hpp
     XAD/MathFunctions.hpp

src/XAD/Hessian.hpp

@@ -0,0 +1,74 @@
+/*******************************************************************************
+
+   Implementation of Hessian computing methods.
+
+   This file is part of XAD, a comprehensive C++ library for
+   automatic differentiation.
+
+   Copyright (C) 2010-2024 Xcelerit Computing Ltd.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Affero General Public License as published
+   by the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Affero General Public License for more details.
+
+   You should have received a copy of the GNU Affero General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+******************************************************************************/
+
+#include <functional>
+#include <vector>
+
+namespace xad
+{
+
+template <class T>
+class Hessian
+{
+  public:
+    Hessian(std::function<T(std::vector<T> &)> func, std::vector<T> &v) : foo(func), v(v) {}
+
+    std::vector<std::vector<T>> compute()
+    {
+        xad::Tape<xad::FReal<double>> tape;
+        domain = static_cast<unsigned int>(v.size());
+        std::vector<std::vector<T>> matrix(domain, std::vector<T>(domain, 0.0));
+
+        tape.registerInputs(v);
+
+        for (unsigned int i = 0; i < domain; i++)
+        {
+            derivative(value(v[i])) = 1.0;
+            tape.newRecording();
+
+            T y = foo(v);
+            tape.registerOutput(y);
+            value(derivative(y)) = 1.0;
+
+            tape.computeAdjoints();
+
+            for (unsigned int j = 0; j < domain; j++)
+            {
+                // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j]))
+                //           << "\n";
+                matrix[i][j] = derivative(derivative(v[j]));
+            }
+
+            derivative(value(v[i])) = 0.0;
+        }
+
+        return matrix;
+    }
+
+  private:
+    std::function<T(std::vector<T> &)> foo;
+    std::vector<T> v;
+    unsigned int domain;
+};
+}  // namespace xad

src/XAD/Jacobian.hpp

@@ -0,0 +1,79 @@
+/*******************************************************************************
+
+   Implementation of Jacobian computing methods.
+
+   This file is part of XAD, a comprehensive C++ library for
+   automatic differentiation.
+
+   Copyright (C) 2010-2024 Xcelerit Computing Ltd.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Affero General Public License as published
+   by the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Affero General Public License for more details.
+
+   You should have received a copy of the GNU Affero General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+******************************************************************************/
+
+#include <functional>
+#include <vector>
+
+namespace xad
+{
+
+template <class T>
+class Jacobian
+{
+  public:
+    Jacobian(std::vector<std::function<T(std::vector<T> &)>> foos, std::vector<T> &v,
+             xad::Tape<double> *tape)
+        : foos(foos),
+          v(v),
+          tape(tape),
+          domain(static_cast<unsigned int>(v.size())),
+          codomain(static_cast<unsigned int>(foos.size()))
+    {
+    }
+
+    std::vector<std::vector<T>> compute()
+    {
+        std::vector<std::vector<T>> matrix(codomain, std::vector<T>(domain, 0.0));
+
+        tape->registerInputs(v);
+
+        for (unsigned int i = 0; i < domain; i++)
+        {
+            for (unsigned int j = 0; j < codomain; j++)
+            {
+                derivative(v[i]) = 1.0;
+                tape->newRecording();
+
+                T y = foos[j](v);
+                tape->registerOutput(y);
+                derivative(y) = 1.0;
+
+                tape->computeAdjoints();
+
+                // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl;
+                matrix[i][j] = derivative(v[i]);
+                derivative(v[i]) = 0.0;
+            }
+        }
+
+        return matrix;
+    }
+
+  private:
+    std::vector<std::function<T(std::vector<T> &)>> foos;
+    std::vector<T> v;
+    Tape<double> *tape;
+    unsigned int domain, codomain;
+};
+}  // namespace xad

src/XAD/XAD.hpp

@@ -35,7 +35,9 @@
 #include <XAD/Config.hpp>
 #include <XAD/Exceptions.hpp>
 #include <XAD/Expression.hpp>
+#include <XAD/Hessian.hpp>
 #include <XAD/Interface.hpp>
+#include <XAD/Jacobian.hpp>
 #include <XAD/Literals.hpp>
 #include <XAD/MathFunctions.hpp>
 #include <XAD/Tape.hpp>

test/CMakeLists.txt

@@ -44,6 +44,8 @@ set(testfiles
     StdCompatibility_test.cpp 
     ReusableRange_test.cpp
     PartialRollback_test.cpp
+    Hessian_test.cpp
+    Jacobian_test.cpp
 )
 
 if (CMAKE_COMPILER_IS_GNUCC)

test/Hessian_test.cpp

@@ -0,0 +1,92 @@
+/*******************************************************************************
+
+   Tests for xad::Hessian class.
+
+   This file is part of XAD, a comprehensive C++ library for
+   automatic differentiation.
+
+   Copyright (C) 2010-2024 Xcelerit Computing Ltd.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Affero General Public License as published
+   by the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Affero General Public License for more details.
+
+   You should have received a copy of the GNU Affero General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+******************************************************************************/
+
+#include <XAD/XAD.hpp>
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+#include <complex>
+#include <type_traits>
+#include <typeinfo>
+#include <vector>
+
+using namespace ::testing;
+
+template <class T>
+T quad(std::vector<T> &x)
+{
+    T c = x[0] * x[0];
+    T d = x[1] * x[1];
+    return c + d;
+}
+
+template <class T>
+T tquad(std::vector<T> &x)
+{
+    T c = x[0] * x[0];
+    T d = x[1] * x[1];
+    T e = x[2] * x[2];
+    return c + d + e;
+}
+
+template <class T>
+T foo(std::vector<T> &x)
+{
+    T c = exp(x[0]);
+    T d = sin(x[1]);
+    T e = cos(x[2]);
+    return c + d + e;
+}
+
+TEST(HessianTest, QuadraticForwardAdjoint)
+{
+    typedef xad::AReal<xad::FReal<double>> AD;
+
+    std::vector<AD> x = {3, 2};
+    xad::Hessian<AD> h(quad<AD>, x);
+
+    std::vector<std::vector<AD>> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}};
+    std::vector<std::vector<AD>> computed_hessian = h.compute();
+
+    for (unsigned int i = 0; i < cross_hessian.size(); i++)
+        for (unsigned int j = 0; j < cross_hessian.size(); j++)
+            ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]);
+}
+
+TEST(HessianTest, ThreeVarQuadratic)
+{
+    typedef xad::fwd_adj<double> mode;
+    typedef mode::active_type AD;
+
+    std::vector<AD> x = {3, 2, 4};
+    xad::Hessian<AD> hes(tquad<AD>, x);
+
+    std::vector<std::vector<AD>> cross_hessian = {
+        {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}};
+    std::vector<std::vector<AD>> computed_hessian = hes.compute();
+
+    for (unsigned int i = 0; i < cross_hessian.size(); i++)
+        for (unsigned int j = 0; j < cross_hessian.size(); j++)
+            ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]);
+}

test/Jacobian_test.cpp

@@ -0,0 +1,65 @@
+/*******************************************************************************
+
+   Tests for xad::Jacobian class.
+
+   This file is part of XAD, a comprehensive C++ library for
+   automatic differentiation.
+
+   Copyright (C) 2010-2024 Xcelerit Computing Ltd.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Affero General Public License as published
+   by the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Affero General Public License for more details.
+
+   You should have received a copy of the GNU Affero General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+******************************************************************************/
+
+#include <XAD/XAD.hpp>
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+#include <functional>
+#include <vector>
+
+template <typename T>
+T foo1(std::vector<T> &x)
+{
+    return x[0] + sin(x[1]);
+}
+
+template <typename T>
+T foo2(std::vector<T> &x)
+{
+    return x[1] + sin(x[0]);
+}
+
+TEST(JacobianTest, SimpleJacobianAdjoint)
+{
+    typedef xad::adj<double> mode;
+    typedef mode::tape_type tape_type;
+    typedef mode::active_type AD;
+
+    tape_type tape;
+
+    std::vector<std::function<AD(std::vector<AD> &)>> funcs(2);
+    funcs[0] = foo1<AD>;
+    funcs[1] = foo2<AD>;
+
+    std::vector<AD> x = {-2, 1};
+    xad::Jacobian<AD> jac(funcs, x, &tape);
+
+    std::vector<std::vector<AD>> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}};
+    std::vector<std::vector<AD>> computed_jacobian = jac.compute();
+
+    for (unsigned int i = 0; i < cross_jacobian.size(); i++)
+        for (unsigned int j = 0; j < cross_jacobian.size(); j++)
+            ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]);
+}