Arbitrary Taylor expansion

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: nn2poly
 Title: Neural Network Weights Transformation into Polynomial Coefficients
-Version: 1.0.1
+Version: 1.0.1.9000
 Authors@R: c(
     person(given = "Pablo", family = "Morala", role = c("aut", "cre"),
            email = "moralapablo@gmail.com",

NEWS.md

@@ -1,3 +1,7 @@
+# nn2poly 1.0.1.9000
+
+- Implement Taylor expansion at any given `taylor_center` (#85 addressing #82).
+
 # nn2poly 1.0.1
 
 - Fix UB in hash implementation (#84).

R/RcppExports.R

@@ -5,16 +5,16 @@ obtain_taylor_vector <- function(taylor_orders, af_string_list) {
     .Call(`_nn2poly_obtain_taylor_vector`, taylor_orders, af_string_list)
 }
 
-obtain_derivatives_list <- function(taylor_orders, af_string_list) {
-    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list)
+obtain_derivatives_list <- function(taylor_orders, af_string_list, taylor_center = 0.0) {
+    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list, taylor_center)
 }
 
 alg_non_linear <- function(coeffs_input, labels_input, labels_output, previous_order, q_layer, g) {
     .Call(`_nn2poly_alg_non_linear`, coeffs_input, labels_input, labels_output, previous_order, q_layer, g)
 }
 
-nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders) {
-    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders)
+nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders, taylor_center = 0.0) {
+    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders, taylor_center)
 }
 
 combinations_with_repetition <- function(n, k) {

R/nn2poly.R

@@ -45,6 +45,9 @@ NULL
 #  layers. If a vector is used, each value corresponds to the Taylor order used
 #' at each layer activation function. Default set to \code{8}.
 #'
+#' @param taylor_center \code{numeric} value that sets the point at which the
+#' Taylor expansion is computed. Default set to \code{0}.
+#'
 #' @param ... Ignored.
 #'
 #' @return Returns an object of class `nn2poly`.
@@ -114,6 +117,7 @@ nn2poly <- function(object,
                     max_order = 2,
                     keep_layers = FALSE,
                     taylor_orders = 8,
+                    taylor_center = 0.0,
                     ...) {
   UseMethod("nn2poly")
 }
@@ -123,9 +127,10 @@ nn2poly.list <- function(object,
                          max_order = 2,
                          keep_layers = FALSE,
                          taylor_orders = 8,
+                         taylor_center = 0.0,
                          ...) {
-  result <- nn2poly_algorithm(object, names(object),
-                              max_order, keep_layers, taylor_orders)
+  result <- nn2poly_algorithm(object, names(object), max_order, keep_layers,
+                              taylor_orders, taylor_center)
   class(result) <- "nn2poly"
   result
 }

src/RcppExports.cpp

@@ -25,14 +25,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // obtain_derivatives_list
-CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list);
-RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP) {
+CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list, double taylor_center);
+RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP, SEXP taylor_centerSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
     Rcpp::traits::input_parameter< const Functions& >::type af_string_list(af_string_listSEXP);
-    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list));
+    Rcpp::traits::input_parameter< double >::type taylor_center(taylor_centerSEXP);
+    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list, taylor_center));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -53,8 +54,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // nn2poly_algorithm
-List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders);
-RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP) {
+List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders, double taylor_center);
+RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP, SEXP taylor_centerSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -63,7 +64,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type max_order(max_orderSEXP);
     Rcpp::traits::input_parameter< bool >::type keep_layers(keep_layersSEXP);
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
-    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders));
+    Rcpp::traits::input_parameter< double >::type taylor_center(taylor_centerSEXP);
+    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders, taylor_center));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -82,9 +84,9 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_nn2poly_obtain_taylor_vector", (DL_FUNC) &_nn2poly_obtain_taylor_vector, 2},
-    {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 2},
+    {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 3},
     {"_nn2poly_alg_non_linear", (DL_FUNC) &_nn2poly_alg_non_linear, 6},
-    {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 5},
+    {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 6},
     {"_nn2poly_combinations_with_repetition", (DL_FUNC) &_nn2poly_combinations_with_repetition, 2},
     {NULL, NULL, 0}
 };

src/approx.h

@@ -0,0 +1,14 @@
+#ifndef nn2poly__approx_h
+#define nn2poly__approx_h
+
+#include "nn2poly_types.h"
+#include "approx_taylor.h"
+
+template <typename... Args>
+void coeffs(const std::string& type, const std::string& af, Coeffs& c, Args&&... args) {
+  if (type == "taylor") if constexpr (sizeof...(Args) < 2)
+    return nn2poly::detail::coeffs_taylor(af, c, std::forward<Args>(args)...);
+  throw std::invalid_argument("approximation type '" + type + "' not supported");
+}
+
+#endif

src/approx_taylor.h

@@ -0,0 +1,64 @@
+#ifndef nn2poly__approx_taylor_h
+#define nn2poly__approx_taylor_h
+
+namespace nn2poly {
+namespace detail {
+
+// Helper to handle the shared Taylor convolution loop structure
+template<typename Func>
+void compute_taylor_conv(Coeffs& c, Func step, double c0) {
+  c[0] = c0;
+  for (size_t n = 0; n < c.size() - 1; n++) {
+    double conv = 0.0;
+    for (size_t k = 0; k <= n; k++)
+      conv += c[k] * c[n - k];
+    c[n + 1] = step(n, conv, c) / static_cast<double>(n + 1);
+  }
+}
+
+inline void coeffs_taylor_sigmoid(Coeffs& c, double a) {
+  compute_taylor_conv(c, [](size_t n, double conv, const Coeffs& c) {
+    return (c[n] - conv);
+  }, 1.0 / (1.0 + std::exp(-a)));
+}
+
+inline void coeffs_taylor_tanh(Coeffs& c, double a) {
+  compute_taylor_conv(c, [](size_t n, double conv, const Coeffs&) {
+    return (n == 0 ? 1.0 : 0.0) - conv;
+  }, std::tanh(a));
+}
+
+inline void coeffs_taylor_linear(Coeffs& c, double a) {
+  c[0] = a;
+  if (c.size() > 1) {
+    c[1] = 1.0;
+    std::fill(c.begin() + 2, c.end(), 0.0);
+  }
+}
+
+inline void coeffs_taylor_softplus(Coeffs& c, double a) {
+  coeffs_taylor_sigmoid(c, a);
+  for (size_t n = c.size() - 1; n > 0; n--)
+    c[n] = c[n - 1] / static_cast<double>(n);
+  if (a > 0.0)
+    c[0] = a + std::log(1.0 + std::exp(-a));
+  else c[0] = std::log(1.0 + std::exp(a));
+}
+
+inline void coeffs_taylor(const std::string& af, Coeffs& c, double a = 0.0) {
+  if (c.empty()) return;
+  if (af == "sigmoid")
+    return coeffs_taylor_sigmoid(c, a);
+  if (af == "tanh")
+    return coeffs_taylor_tanh(c, a);
+  if (af == "softplus")
+    return coeffs_taylor_softplus(c, a);
+  if (af == "linear")
+    return coeffs_taylor_linear(c, a);
+  throw std::invalid_argument("function '" + af + "' not supported");
+}
+
+} // namespace detail
+} // namespace nn2poly
+
+#endif

src/linalg_arma.h

@@ -28,6 +28,10 @@ inline Weights alg_linear(Weights& coeffs_list, const Weights& layer) {
   return mat;
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.col(i) -= scalar;
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.col(i) += scalar * vec;
 }

src/linalg_eigen.h

@@ -28,6 +28,10 @@ inline Weights alg_linear(const Weights& coeffs_list, const Weights& layer) {
   return mat;
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.col(i) -= scalar;
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.col(i) += scalar * vec;
 }

src/linalg_torch.h

@@ -30,6 +30,10 @@ inline Weights alg_linear(Weights& coeffs_list, const Weights& layer) {
   return mat;
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.select(1, i).sub_(scalar);
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.select(1, i).add_(vec * scalar);
 }

src/nn2poly.cpp

@@ -1,6 +1,6 @@
 #include "partitions.h"
 #include "multiset.h"
-#include "taylor.h"
+#include "approx.h"
 using namespace nn2poly::linalg;
 
 // [[Rcpp::export]]
@@ -24,7 +24,8 @@ Term obtain_taylor_vector(const Term& taylor_orders,
 
 // [[Rcpp::export]]
 CoeffsList obtain_derivatives_list(const Term& taylor_orders,
-                                   const Functions& af_string_list) {
+                                   const Functions& af_string_list,
+                                   double taylor_center = 0.0) {
   if (taylor_orders.size() != af_string_list.size())
     throw std::invalid_argument(
       "`taylor_orders` length does not match provided number of layers");
@@ -34,8 +35,9 @@ CoeffsList obtain_derivatives_list(const Term& taylor_orders,
     if (taylor_orders[i] < 0)
       throw std::invalid_argument("`taylor_orders` must be non-negative");
     // Obtain the vector with the derivatives of the activation function up to
-    // the given degree centered at 0
-    out[i] = coeffs_taylor(af_string_list[i], taylor_orders[i]);
+    // the given degree centered at taylor_center
+    out[i].resize(static_cast<size_t>(taylor_orders[i] + 1));
+    coeffs("taylor", af_string_list[i], out[i], taylor_center);
   }
 
   return out;
@@ -135,7 +137,7 @@ inline void check_weights_dimensions(const Layers& layers) {
 // [[Rcpp::export]]
 List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
                        int max_order, bool keep_layers,
-                       const Term& taylor_orders) {
+                       const Term& taylor_orders, double taylor_center = 0.0) {
   if (layers.empty())
     throw std::invalid_argument("`layers` is empty");
   if (af_list.empty())
@@ -154,7 +156,7 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
   // depending on the last layer being linear or not
   WeightsLists results(static_cast<size_t>(last_linear ? 2 * L - 1 : 2 * L));
   const Term taylor = obtain_taylor_vector(taylor_orders, af_list);
-  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list);
+  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list, taylor_center);
 
   // Starting point for the algorithm: Set weights as coefficients
   // of an order 1 polynomial.
@@ -222,6 +224,9 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
 
     NN2POLY_DEBUG_LOG(2, "[layer", l, "]",
       DTAG(previous_order), DTAG(q_layer), DTAG(labels_list.size()));
+
+    // Subtract the center from the intercept
+    if (taylor_center) sub_scalar(coeffs_list, 0, taylor_center);
     // Apply non-linear algorithm
     coeffs_list = alg_non_linear_impl(
       coeffs_list,