Enhance finite_jacobian functionality with tensor order support

- Follows the convention that even-order tensors (e.g., matrices) are stored in column-blocks, while odd-order tensors are stored in row-blocks. This is based on the tensor vectorization presented in "Dynamic Deformables" by Kim and Eberle [2022].
- Update CMake configuration to use C++17
- Add tests for size mismatch and tensor layouts

Author: zfergus

CMakeLists.txt

@@ -27,7 +27,7 @@ option(FINITE_DIFF_BUILD_UNIT_TESTS  "Build unit-tests"  ${FINITE_DIFF_TOPLEVEL_
 
 # Set default minimum C++ standard
 if(FINITE_DIFF_TOPLEVEL_PROJECT)
-    set(CMAKE_CXX_STANDARD 11)
+    set(CMAKE_CXX_STANDARD 17)
     set(CMAKE_CXX_STANDARD_REQUIRED ON)
     set(CMAKE_CXX_EXTENSIONS OFF)
 endif()
@@ -76,8 +76,8 @@ target_link_libraries(finitediff_finitediff PUBLIC spdlog::spdlog)
 # Compiler options
 ################################################################################
 
-# Use C++11
-target_compile_features(finitediff_finitediff PUBLIC cxx_std_11)
+# Use C++17
+target_compile_features(finitediff_finitediff PUBLIC cxx_std_17)
 
 ################################################################################
 # Tests

src/finitediff.cpp

@@ -92,6 +92,7 @@ void finite_gradient(
     }
 }
 
+template <bool IsTensorOrderEven>
 void finite_jacobian(
     const Eigen::Ref<const Eigen::VectorXd>& x,
     const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
@@ -113,18 +114,31 @@ void finite_jacobian(
         const Eigen::MatrixXd tmp = f(x);
         f_rows = tmp.rows(), f_cols = tmp.cols();
     }
-    jac.setZero(f_rows, f_cols * x.size());
+    if constexpr (IsTensorOrderEven) {
+        jac.setZero(f_rows, f_cols * x.size());
+    } else {
+        jac.setZero(f_rows * x.size(), f_cols);
+    }
 
     // f: ℝ^n ↦ ℝ^{p×q} ⟹ ∇f: ℝ^n ↦ ℝ^{p×(qn)}
     Eigen::VectorXd x_mutable = x;
     for (size_t i = 0; i < x.size(); i++) {
         for (size_t ci = 0; ci < inner_steps; ci++) {
             x_mutable[i] += internal_coeffs[ci] * eps;
-            jac.middleCols(f_cols * i, f_cols) +=
-                external_coeffs[ci] * f(x_mutable);
+            if constexpr (IsTensorOrderEven) {
+                jac.middleCols(f_cols * i, f_cols) +=
+                    external_coeffs[ci] * f(x_mutable);
+            } else {
+                jac.middleRows(f_rows * i, f_rows) +=
+                    external_coeffs[ci] * f(x_mutable);
+            }
             x_mutable[i] = x[i];
         }
-        jac.middleCols(f_cols * i, f_cols) /= denom;
+        if constexpr (IsTensorOrderEven) {
+            jac.middleCols(f_cols * i, f_cols) /= denom;
+        } else {
+            jac.middleRows(f_rows * i, f_rows) /= denom;
+        }
     }
 }
 
@@ -172,7 +186,12 @@ bool compare_gradient(
     const double test_eps,
     const std::string& msg)
 {
-    assert(x.rows() == y.rows());
+    if (x.rows() != y.rows()) {
+        spdlog::debug(
+            "{} gradient size mismatch: x.rows()={} y.rows()={}", msg, x.rows(),
+            y.rows());
+        return false;
+    }
 
     bool same = true;
     for (long i = 0; i < x.rows(); i++) {
@@ -199,8 +218,13 @@ bool compare_jacobian(
     const double test_eps,
     const std::string& msg)
 {
-    assert(x.rows() == y.rows());
-    assert(x.cols() == y.cols());
+    if (x.rows() != y.rows() || x.cols() != y.cols()) {
+        spdlog::debug(
+            "{} jacobian size mismatch: x.rows()={} x.cols()={} "
+            "y.rows()={} y.cols()={}",
+            msg, x.rows(), x.cols(), y.rows(), y.cols());
+        return false;
+    }
 
     bool same = true;
     for (long i = 0; i < x.rows(); i++) {
@@ -256,4 +280,9 @@ Eigen::MatrixXd unflatten(const Eigen::Ref<const Eigen::VectorXd>& x, int dim)
     return X;
 }
 
+// clang-format off
+template void finite_jacobian<false>(const Eigen::Ref<const Eigen::VectorXd>& x, const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f, Eigen::MatrixXd& jac, const AccuracyOrder accuracy, const double eps);
+template void finite_jacobian<true>(const Eigen::Ref<const Eigen::VectorXd>& x, const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f, Eigen::MatrixXd& jac, const AccuracyOrder accuracy, const double eps);
+// clang-format on
+
 } // namespace fd

src/finitediff.hpp

@@ -41,18 +41,48 @@ void finite_gradient(
 /**
  * @brief Compute the Jacobian of a function using finite differences.
  *
+ * @param[in]  x              Point at which to compute the Jacobian.
+ * @param[in]  f              Compute the Jacobian of this function.
+ * @param[out] jac            Computed Jacobian.
+ * @param[in]  accuracy       Accuracy of the finite differences.
+ * @param[in]  eps            Value of the finite difference step.
+ * @tparam IsTensorOrderEven  If true, the Jacobian is stored in column-blocks.
+ *                            Otherwise, in row-blocks.
+ */
+template <bool IsTensorOrderEven = true>
+void finite_jacobian(
+    const Eigen::Ref<const Eigen::VectorXd>& x,
+    const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
+    Eigen::MatrixXd& jac,
+    const AccuracyOrder accuracy = SECOND,
+    const double eps = 1.0e-8);
+
+/**
+ * @brief Compute the Jacobian of a function using finite differences with
+ *        explicit tensor order.
+ *
+ * Follows the convention that even-order tensors (e.g., matrices) are stored in
+ * column-blocks, while odd-order tensors are stored in row-blocks. This is
+ * based on the tensor vectorization presented in "Dynamic Deformables" by Kim
+ * and Eberle [2022].
+ *
  * @param[in]  x         Point at which to compute the Jacobian.
  * @param[in]  f         Compute the Jacobian of this function.
  * @param[out] jac       Computed Jacobian.
  * @param[in]  accuracy  Accuracy of the finite differences.
  * @param[in]  eps       Value of the finite difference step.
+ * @tparam TensorOrder   Order of the output tensor of the function f.
  */
-void finite_jacobian(
+template <int TensorOrder>
+inline void finite_jacobian_tensor(
     const Eigen::Ref<const Eigen::VectorXd>& x,
     const std::function<Eigen::MatrixXd(const Eigen::VectorXd&)>& f,
     Eigen::MatrixXd& jac,
     const AccuracyOrder accuracy = SECOND,
-    const double eps = 1.0e-8);
+    const double eps = 1.0e-8)
+{
+    return finite_jacobian<TensorOrder % 2 == 0>(x, f, jac, accuracy, eps);
+}
 
 /**
  * @brief Compute the Hessian of a function using finite differences.

tests/CMakeLists.txt

@@ -5,6 +5,8 @@
 add_executable(finitediff_tests
   test_gradient.cpp
   test_jacobian.cpp
+  test_compare_size_mismatch.cpp
+  test_jacobian_tensor.cpp
   test_hessian.cpp
   test_flatten.cpp
 )

tests/test_compare_size_mismatch.cpp

@@ -0,0 +1,15 @@
+#include <catch2/catch_test_macros.hpp>
+
+#include <Eigen/Core>
+
+#include <finitediff.hpp>
+
+using namespace fd;
+
+TEST_CASE("compare_jacobian returns false for different sizes", "[compare]")
+{
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 3);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 2);
+
+    CHECK_FALSE(compare_jacobian(A, B));
+}

tests/test_gradient.cpp

@@ -9,7 +9,7 @@
 
 using namespace fd;
 
-TEST_CASE("Test finite difference gradient of quadratic", "[gradient]")
+TEST_CASE("Finite difference gradient of quadratic", "[gradient]")
 {
     int n = GENERATE(1, 2, 4, 10, 100);
 
@@ -33,7 +33,7 @@ TEST_CASE("Test finite difference gradient of quadratic", "[gradient]")
     CHECK(compare_gradient(grad, fgrad));
 }
 
-TEST_CASE("Test finite difference gradient of Rosenbrock", "[gradient]")
+TEST_CASE("Finite difference gradient of Rosenbrock", "[gradient]")
 {
     const auto f = [](const Eigen::VectorXd& x) {
         double t1 = 1 - x[0];
@@ -59,7 +59,7 @@ TEST_CASE("Test finite difference gradient of Rosenbrock", "[gradient]")
     CHECK(compare_gradient(grad, fgrad));
 }
 
-TEST_CASE("Test finite difference gradient of trig", "[gradient]")
+TEST_CASE("Finite difference gradient of trig", "[gradient]")
 {
     int n = GENERATE(1, 2, 4, 10, 100);
 

tests/test_hessian.cpp

@@ -10,7 +10,7 @@
 
 using namespace fd;
 
-TEST_CASE("Test finite difference hessian of quadratic", "[hessian]")
+TEST_CASE("Finite difference hessian of quadratic", "[hessian]")
 {
     AccuracyOrder accuracy = GENERATE(SECOND, FOURTH, SIXTH, EIGHTH);
 
@@ -35,7 +35,7 @@ TEST_CASE("Test finite difference hessian of quadratic", "[hessian]")
     CHECK(compare_hessian(hess, fhess));
 }
 
-TEST_CASE("Test finite difference hessian of Rosenbrock", "[hessian]")
+TEST_CASE("Finite difference hessian of Rosenbrock", "[hessian]")
 {
     AccuracyOrder accuracy = GENERATE(SECOND, FOURTH, SIXTH, EIGHTH);
     const auto f = [](const Eigen::VectorXd& x) {
@@ -58,7 +58,7 @@ TEST_CASE("Test finite difference hessian of Rosenbrock", "[hessian]")
     CHECK(compare_hessian(hess, fhess));
 }
 
-TEST_CASE("Test finite difference hessian of trig", "[hessian]")
+TEST_CASE("Finite difference hessian of trig", "[hessian]")
 {
     AccuracyOrder accuracy = GENERATE(SECOND, FOURTH, SIXTH, EIGHTH);
     int n = GENERATE(1, 2, 4, 10, 25);

tests/test_jacobian.cpp

@@ -9,7 +9,7 @@
 
 using namespace fd;
 
-TEST_CASE("Test finite difference jacobian of linear", "[jacobian]")
+TEST_CASE("Finite difference jacobian of linear", "[jacobian]")
 {
     int n = GENERATE(1, 2, 4, 10, 100);
 
@@ -32,7 +32,7 @@ TEST_CASE("Test finite difference jacobian of linear", "[jacobian]")
     CHECK(compare_jacobian(jac, fjac));
 }
 
-TEST_CASE("Test finite difference jacobian of trig", "[jacobian]")
+TEST_CASE("Finite difference jacobian of trig", "[jacobian]")
 {
     int n = GENERATE(1, 2, 4, 10, 100);
 

tests/test_jacobian_tensor.cpp

@@ -0,0 +1,52 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators_all.hpp>
+
+#include <Eigen/Core>
+
+#include <finitediff.hpp>
+
+using namespace fd;
+
+TEST_CASE("finite_jacobian_tensor layouts", "[jacobian_tensor]")
+{
+    int n = GENERATE(1, 2, 4, 10);
+
+    const int p = 2;
+    const int q = 3;
+
+    // Construct a tensor T_k (p x q) for each input component k such that
+    // f(x) = sum_k x_k * T_k. The analytic derivative w.r.t x_k is T_k.
+    std::vector<Eigen::MatrixXd> T;
+    T.reserve(n);
+    for (int k = 0; k < n; ++k) {
+        T.emplace_back(Eigen::MatrixXd::Random(p, q));
+    }
+
+    const auto f = [&](const Eigen::VectorXd& x) -> Eigen::MatrixXd {
+        Eigen::MatrixXd R = Eigen::MatrixXd::Zero(p, q);
+        for (int k = 0; k < n; ++k) {
+            R += x[k] * T[k];
+        }
+        return R;
+    };
+
+    Eigen::VectorXd x = Eigen::VectorXd::Random(n);
+
+    // Analytic Jacobians for the two storage conventions
+    Eigen::MatrixXd jac_even = Eigen::MatrixXd::Zero(p, q * n);
+    Eigen::MatrixXd jac_odd = Eigen::MatrixXd::Zero(p * n, q);
+    for (int k = 0; k < n; ++k) {
+        jac_even.block(0, q * k, p, q) = T[k];
+        jac_odd.block(p * k, 0, p, q) = T[k];
+    }
+
+    AccuracyOrder accuracy = GENERATE(SECOND, FOURTH, SIXTH, EIGHTH);
+
+    Eigen::MatrixXd fjac_even;
+    finite_jacobian_tensor<2>(x, f, fjac_even, accuracy);
+    CHECK(compare_jacobian(jac_even, fjac_even));
+
+    Eigen::MatrixXd fjac_odd;
+    finite_jacobian_tensor<1>(x, f, fjac_odd, accuracy);
+    CHECK(compare_jacobian(jac_odd, fjac_odd));
+}