Custom non-convex optimizer (PSGD)

.gitignore

@@ -2,3 +2,6 @@ build
 cmake-build-debug
 venv
 .idea
+*.so
+*.pyc
+__pycache__/

CMakeLists.txt

@@ -18,6 +18,8 @@ file(GLOB SOURCES "src/cpp/*/*.cpp")
 # Create the main C++ library target with a unique name
 add_library(finmath_library SHARED ${SOURCES}
     "src/cpp/InterestAndAnnuities/simple_interest.cpp"
+    "src/cpp/Optimization/psgd.cpp"             # Add PSGD source
+    "src/cpp/Optimization/linalg_helpers.cpp" # Add linalg helpers
     "include/finmath/InterestAndAnnuities/simple_interest.h"
     "include/finmath/OptionPricing/options_pricing.h"
     "include/finmath/OptionPricing/options_pricing_types.h"
@@ -26,43 +28,46 @@ add_library(finmath_library SHARED ${SOURCES}
     "include/finmath/TimeSeries/rsi.h"
     "include/finmath/TimeSeries/ema.h")
 
-# Test executables
-add_executable(black_scholes_test test/OptionPricing/black_scholes_test.cpp)
-target_link_libraries(black_scholes_test finmath_library)
-
-add_executable(binomial_option_pricing_test test/OptionPricing/binomial_option_pricing_test.cpp)
-target_link_libraries(binomial_option_pricing_test finmath_library)
+# Link pybind11 headers to the main library (needed for numpy integration in C++ files)
+target_link_libraries(finmath_library PUBLIC pybind11::headers)
 
-add_executable(compound_interest_test test/InterestAndAnnuities/compound_interest_test.cpp)
-target_link_libraries(compound_interest_test finmath_library)
+# Also link Python libraries/headers (needed for Python.h)
+find_package(Python COMPONENTS Interpreter Development REQUIRED)
+target_link_libraries(finmath_library PUBLIC Python::Python)
 
-add_executable(rsi_test test/TimeSeries/rsi_test.cpp)
-target_link_libraries(rsi_test finmath_library)
-
-# Test runner
-add_executable(run_all_tests test/test_runner.cpp)
+# Test executables
+# add_executable(black_scholes_test test/OptionPricing/black_scholes_test.cpp)
+# target_link_libraries(black_scholes_test finmath_library)
+# 
+# add_executable(binomial_option_pricing_test test/OptionPricing/binomial_option_pricing_test.cpp)
+# target_link_libraries(binomial_option_pricing_test finmath_library)
+# 
+# add_executable(compound_interest_test test/InterestAndAnnuities/compound_interest_test.cpp)
+# target_link_libraries(compound_interest_test finmath_library)
+# 
+# add_executable(rsi_test test/TimeSeries/rsi_test.cpp) # This was the problematic one
+# target_link_libraries(rsi_test finmath_library)
+# 
+# # Test runner (can be removed if using ctest directly)
+# add_executable(run_all_tests test/test_runner.cpp)
 
 # Enable testing
 enable_testing()
 
-# Define individual tests
-add_test(NAME BlackScholesTest COMMAND black_scholes_test)
-add_test(NAME BinomialOptionPricingTest COMMAND binomial_option_pricing_test)
-add_test(NAME CompoundInterestTest COMMAND compound_interest_test)
-add_test(NAME RSITest COMMAND rsi_test)
-
-# Add a custom target to run all tests
-add_custom_target(build_and_test
-    COMMAND ${CMAKE_COMMAND} --build . --target black_scholes_test
-    COMMAND ${CMAKE_COMMAND} --build . --target binomial_option_pricing_test
-    COMMAND ${CMAKE_COMMAND} --build . --target compound_interest_test
-    COMMAND ${CMAKE_COMMAND} --build . --target rsi_test
-    COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
-    WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
-)
-
-# Make 'build_and_test' the default target
-add_custom_target(default ALL DEPENDS build_and_test)
+# Helper macro to add a test executable and link it
+macro(add_cpp_test test_name source_file)
+    message(STATUS "Adding C++ test: ${test_name} from ${source_file}")
+    add_executable(${test_name}_executable ${source_file})
+    target_link_libraries(${test_name}_executable PRIVATE finmath_library)
+    add_test(NAME ${test_name} COMMAND ${test_name}_executable)
+endmacro()
+
+# Add C++ tests using the macro
+add_cpp_test(RSITest test/TimeSeries/RSI/C++/rsi_test.cpp)
+add_cpp_test(RollingVolatilityTest test/TimeSeries/RollingVolatility/C++/rolling_volatility_test.cpp)
+add_cpp_test(SimpleMovingAverageTest test/TimeSeries/SimpleMovingAverage/C++/simple_moving_average_test.cpp)
+add_cpp_test(EMATest test/TimeSeries/EMA/C++/ema_test.cpp)
+# Add tests for EMA here later...
 
 # Add pybind11 for Python bindings
 include(FetchContent)
@@ -79,5 +84,8 @@ pybind11_add_module(finmath_bindings src/python_bindings.cpp ${SOURCES})
 # Set the output name of the bindings to 'finmath' to match your desired module name
 set_target_properties(finmath_bindings PROPERTIES OUTPUT_NAME "finmath")
 
+# Set the library output directory to be alongside the source bindings file
+set_target_properties(finmath_bindings PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${PROJECT_SOURCE_DIR}/src")
+
 # Link the Python bindings target with the C++ library
 target_link_libraries(finmath_bindings PRIVATE finmath_library)

include/finmath/Optimization/linalg_helpers.h

@@ -0,0 +1,45 @@
+#ifndef FINMATH_OPTIMIZATION_LINALG_HELPERS_H
+#define FINMATH_OPTIMIZATION_LINALG_HELPERS_H
+
+#include <vector>
+#include <cmath>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+
+namespace finmath
+{
+    namespace optimization
+    {
+        namespace linalg
+        {
+
+            // Calculate L2 norm (Euclidean norm)
+            double norm(const std::vector<double> &v);
+
+            // In-place vector addition: a += b
+            void add_vectors(std::vector<double> &a, const std::vector<double> &b);
+
+            // In-place vector subtraction: a -= b
+            void subtract_vectors(std::vector<double> &a, const std::vector<double> &b);
+
+            // In-place scalar multiplication: v *= scalar
+            void scale_vector(std::vector<double> &v, double scalar);
+
+            // In-place addition with scaling: a += scalar * b
+            void add_scaled_vector(std::vector<double> &a, const std::vector<double> &b, double scalar);
+
+            // In-place subtraction with scaling: a -= scalar * b
+            void subtract_scaled_vector(std::vector<double> &a, const std::vector<double> &b, double scalar);
+
+            // Clip vector v in-place if its norm exceeds max_norm
+            void clip_vector_norm(std::vector<double> &v, double max_norm);
+
+            // Sample a vector uniformly from a ball of given radius and dimension
+            std::vector<double> sample_uniform_ball(double radius, size_t dimension, std::mt19937 &rng_engine);
+
+        } // namespace linalg
+    } // namespace optimization
+} // namespace finmath
+
+#endif // FINMATH_OPTIMIZATION_LINALG_HELPERS_H

include/finmath/Optimization/psgd.h

@@ -0,0 +1,62 @@
+#ifndef FINMATH_OPTIMIZATION_PSGD_H
+#define FINMATH_OPTIMIZATION_PSGD_H
+
+#include <vector>
+#include <functional>
+#include <stdexcept> // For potential exceptions
+
+namespace finmath
+{
+    namespace optimization
+    {
+
+        /**
+         * @brief Implements the Enhanced Perturbed Stochastic Gradient Descent (PSGD-C) algorithm.
+         *
+         * Based on the algorithm described in the user-provided LaTeX source and Python implementation.
+         * Designed for non-convex optimization, incorporating EMA smoothing, gradient/parameter clipping,
+         * and noise injection to escape saddle points.
+         *
+         * @param stochastic_grad Function providing a stochastic gradient estimate for a given point x.
+         *                        Signature: std::vector<double>(const std::vector<double>& x)
+         * @param objective_f Function providing the objective function value f(x). Needed for threshold calculation.
+         *                     Signature: double(const std::vector<double>& x)
+         * @param x0 Initial starting point (vector).
+         * @param ell Smoothness parameter (Lipschitz constant of the gradient).
+         * @param rho Hessian Lipschitz parameter.
+         * @param eps Target accuracy for the norm of the smoothed gradient (termination criterion).
+         * @param sigma Standard deviation estimate of the noise in the stochastic gradient samples (||grad_f_i(x) - grad_f(x)||).
+         * @param delta Confidence parameter (probability of failure).
+         * @param batch_size Mini-batch size used by stochastic_grad (influences g_thresh).
+         * @param step_size_coeff Coefficient 'c' to calculate step size eta = c / ell.
+         * @param ema_beta Decay factor for the exponential moving average of the gradient (typically 0.8-0.95).
+         * @param max_iters Maximum number of iterations to perform.
+         * @param grad_clip_norm Maximum L2 norm for the raw stochastic gradient (g_max). Set <= 0 to disable.
+         * @param param_clip_norm Maximum L2 norm for the parameter vector x (x_max). Set <= 0 to disable.
+         * @return The final optimized parameter vector x.
+         *
+         * @throws std::invalid_argument if input parameters are inconsistent (e.g., ell <= 0, rho <= 0, eps <= 0, sigma < 0, 0 < delta < 1, batch_size <= 0, c <= 0, 0 <= ema_beta < 1).
+         */
+        std::vector<double> perturbed_sgd(
+            const std::function<std::vector<double>(const std::vector<double> &)> &stochastic_grad,
+            const std::function<double(const std::vector<double> &)> &objective_f,
+            const std::vector<double> &x0,
+            // Problem params
+            double ell,
+            double rho,
+            double eps = 1e-3,
+            double sigma = 0.1,
+            double delta = 0.1,
+            // Algo params
+            int batch_size = 32,
+            double step_size_coeff = 0.5,
+            double ema_beta = 0.9,
+            int max_iters = 100000,
+            double grad_clip_norm = 10.0,  // Default g_max (or disable if <= 0)
+            double param_clip_norm = 100.0 // Default x_max (or disable if <= 0)
+        );
+
+    } // namespace optimization
+} // namespace finmath
+
+#endif // FINMATH_OPTIMIZATION_PSGD_H

include/finmath/TimeSeries/ema.h

@@ -2,11 +2,18 @@
 #define EMA_H
 
 #include <vector>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
 
 // Function to compute the Exponential Moving Average (EMA) using window size
-std::vector<double> compute_ema(const std::vector<double>& prices, size_t window);
+std::vector<double> compute_ema(const std::vector<double> &prices, size_t window);
 
 // Function to compute the Exponential Moving Average (EMA) using a smoothing factor
-std::vector<double> compute_ema_with_smoothing(const std::vector<double>& prices, double smoothing_factor);
+std::vector<double> compute_ema_with_smoothing(const std::vector<double> &prices, double smoothing_factor);
+
+// NumPy overloads
+std::vector<double> compute_ema_np(py::array_t<double> prices_arr, size_t window);
+std::vector<double> compute_ema_with_smoothing_np(py::array_t<double> prices_arr, double smoothing_factor);
 
 #endif // EMA_H

include/finmath/TimeSeries/rolling_volatility.h

@@ -2,14 +2,20 @@
 #define ROLLING_VOLATILITY_H
 
 #include <vector>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
 
 // Function to compute the logarithmic returns from prices
-std::vector<double> compute_log_returns(const std::vector<double>& prices);
+std::vector<double> compute_log_returns(const std::vector<double> &prices);
 
 // Function to compute the standard deviation of a vector
-double compute_std(const std::vector<double>& data);
+double compute_std(const std::vector<double> &data);
+
+// Function to compute the rolling volatility from a time series of prices (vector version)
+std::vector<double> rolling_volatility(const std::vector<double> &prices, size_t window_size);
 
-// Function to compute the rolling volatility from a time series of prices
-std::vector<double> rolling_volatility(const std::vector<double>& prices, size_t window_size);
+// Overloaded function to compute rolling volatility from a NumPy array
+std::vector<double> rolling_volatility_np(py::array_t<double> prices_arr, size_t window_size);
 
 #endif // ROLLING_VOLATILITY_H

include/finmath/TimeSeries/rsi.h

@@ -1,15 +1,21 @@
 #ifndef RSI_H
 #define RSI_H
 
-#include<vector>
+#include <vector>
+#include <pybind11/numpy.h>
 
-//function to compute the average gain over a window
-double compute_avg_gain(const std::vector<double>& price_changes, size_t window_size);
+namespace py = pybind11;
+
+// function to compute the average gain over a window
+double compute_avg_gain(const std::vector<double> &price_changes, size_t window_size);
 
 // Function to compute the average loss over a window
-double compute_avg_loss(const std::vector<double>& price_changes, size_t window_size);
+double compute_avg_loss(const std::vector<double> &price_changes, size_t window_size);
 
 // Function to compute the RSI from a time series of prices
-std::vector<double> compute_smoothed_rsi(const std::vector<double>& prices, size_t window_size);
+std::vector<double> compute_smoothed_rsi(const std::vector<double> &prices, size_t window_size);
+
+// NumPy overload
+std::vector<double> compute_smoothed_rsi_np(py::array_t<double> prices_arr, size_t window_size);
 
 #endif // RSI_H

include/finmath/TimeSeries/simple_moving_average.h

@@ -2,8 +2,14 @@
 #define SIMPLE_MOVING_AVERAGE_H
 
 #include <vector>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
 
 // Function to compute the moving average from a time series
-std::vector<double> simple_moving_average(const std::vector<double>& data, size_t window_size);
+std::vector<double> simple_moving_average(const std::vector<double> &data, size_t window_size);
+
+// NumPy overload
+std::vector<double> simple_moving_average_np(py::array_t<double> data_arr, size_t window_size);
 
 #endif // MOVING_AVERAGE_H