diff --git a/include/math/sqrt.hpp b/include/math/sqrt.hpp new file mode 100644 index 0000000..572bb51 --- /dev/null +++ b/include/math/sqrt.hpp @@ -0,0 +1,163 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: math + * File Name: sqrt.cpp + * Revision: 0.1.0 + * Date: 06-08-2026 + * Author: Michelle Bausager + * + * Description: + * Functions to calculate the sqrt of numbers + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once + +#include // for panic::vector +#include // for panic::matrix + +namespace panic{ +namespace math{ + + + +/** + * @brief calculates the sqrt of a value. + * + * Computes: + * @code + * sqrt(x, y) + * @endcode + * + * @tparam T Numeric element type. + * @param x Value to take the sqrt of. + * @param y Result. + * + * + * @note This function is omp-friendly. + */ +template +bool sqrt(const T x, T& y); + + +/** + * @brief calculates the sqrt of a value. + * + * Computes: + * @code + * result = sqrt(k) + * @endcode + * + * @tparam T Numeric element type. + * @param x Value to take the sqrt of. + * + * @return The calculated value + * + * @note This function is omp-friendly. + */ +template +T sqrt(const T x); + + +/** + * @brief Calculates the sqrt elementwise in a vector + * + * Computes: + * @code + * c[i] = sqrt(a[i]) + * @endcode + * + * @tparam T Numeric element type. + * @param a Input vector. + * @param c Output vector. Resized to match @p a. + * + * @return true if @p c was resized and filled successfully. + * @return false if resizing @p c failed. + * + * @note This overload writes the result into an existing vector to avoid + * unnecessary temporary allocations. + */ +template +bool sqrt(const panic::tensor::vector& a, panic::tensor::vector& c); + +/** + * @brief Calculates the sqrt elementwise in a vector + * + * Computes: + * @code + * result[i] = sqrt(a[i]) + * @endcode + * + * @tparam T Numeric element type. + * @param a Input vector. + * + * @return A new vector containing the result. + * @return An empty vector if the operation fails. + * + * @note This overload is convenient, but may allocate a new vector. + */ +template +panic::tensor::vector sqrt(const panic::tensor::vector& a); + +/** + * @brief Calculates the sqrt elementwise of a matrix + * + * Computes: + * @code + * C(i,j) = sqrt(A(i,j)) + * @endcode + * + * @tparam T Numeric element type. + * @param A Input matrix. + * @param C Output Matrix. Resized to match @p A. + * + * @return true if @p C was resized and filled successfully. + * @return false if resizing @p C failed. + * + * @note This overload writes the result into an existing vector to avoid + * unnecessary temporary allocations. + */ +template +bool sqrt(const panic::tensor::matrix& A, panic::tensor::matrix& C); + +/** + * @brief Returns the calculated sqrt elementwise of the matrix + * + * Computes: + * @code + * result(i,j) = sqrt(A(i,j)) + * @endcode + * + * @tparam T Numeric element type. + * @param A Input matrix. + * + * @return A new matrix containing the result. + * @return An empty matrix if the operation fails. + * + * @note This overload is convenient, but may allocate a new vector. + */ +template +panic::tensor::matrix sqrt(const panic::tensor::matrix& A); + +} // namespace math +} // namespace panic \ No newline at end of file diff --git a/include/neural_network/layer/trainable_layer.hpp b/include/neural_network/layer/trainable_layer.hpp index 1566b16..c37ce7f 100644 --- a/include/neural_network/layer/trainable_layer.hpp +++ b/include/neural_network/layer/trainable_layer.hpp @@ -63,6 +63,24 @@ struct trainable_layer:layer{ panic::tensor::real_matrix dweights; panic::tensor::real_vector dbiases; + /** + * @brief Previous parameter updates used by momentum SGD. + * + * These remain empty unless an optimizer using momentum + * initializes them. + */ + panic::tensor::real_matrix weight_momentums; + panic::tensor::real_vector bias_momentums; + + /** + * @brief Previous parameter updates used by AdaGrad. + * + * These remain empty unless an optimizer using cache + * initializes them. + */ + panic::tensor::real_matrix weight_cache; + panic::tensor::real_vector bias_cache; + virtual ~trainable_layer() = default; diff --git a/include/neural_network/model/model.hpp b/include/neural_network/model/model.hpp index fefa82a..f7df230 100644 --- a/include/neural_network/model/model.hpp +++ b/include/neural_network/model/model.hpp @@ -323,9 +323,70 @@ struct model{ * * */ - bool add_optimizer_sgd(const panic::types::real_t learning_rate = static_cast(1)); + bool add_optimizer_sgd(const panic::types::real_t learning_rate = 1, + const panic::types::real_t decay = 0, + const panic::types::real_t momentum = 0); + /** + * @brief Adds optimizer_adagrad to the model + * + * Computes: + * @code + * model.optimizer_adagrad(1e-4) + * @endcode + * + * @param learning_rate Learning rate for update_param (default 1e-3). + * + * @return true If looped and optimized every trainable layer. + * + * + */ + bool add_optimizer_adagrad(const panic::types::real_t learning_rate = 1, + const panic::types::real_t decay = 0, + const panic::types::real_t epsilon = 1e-7); + + + /** + * @brief Adds optimizer_rmsprop to the model + * + * Computes: + * @code + * model.optimizer_adagrad(1e-4) + * @endcode + * + * @param learning_rate Learning rate for update_param (default 1e-3). + * + * @return true If looped and optimized every trainable layer. + * + * + */ + bool add_optimizer_rmsprop(const panic::types::real_t learning_rate = 0.001, + const panic::types::real_t decay = 0, + const panic::types::real_t epsilon = 1e-7, + const panic::types::real_t rho = 0.9); + + + + /** + * @brief Adds optimizer_adam to the model + * + * Computes: + * @code + * model.optimizer_adam(1e-4) + * @endcode + * + * @param learning_rate Learning rate for update_param (default 1e-3). + * + * @return true If looped and optimized every trainable layer. + * + * + */ + bool add_optimizer_adam(const panic::types::real_t learning_rate = 0.001, + const panic::types::real_t decay = 0, + const panic::types::real_t epsilon = 1e-7, + const panic::types::real_t beta_1 = 0.9, + const panic::types::real_t beta_2 = 0.999); /** * @brief Finalizes the model configuration. diff --git a/include/neural_network/optimizers/optimizer.hpp b/include/neural_network/optimizers/optimizer.hpp index ff08e61..f548aca 100644 --- a/include/neural_network/optimizers/optimizer.hpp +++ b/include/neural_network/optimizers/optimizer.hpp @@ -56,6 +56,8 @@ namespace panic{ */ struct optimizer{ + panic::types::real_t current_learning_rate; + /** * @brief Default de-constructor * @@ -64,15 +66,33 @@ struct optimizer{ /** - * @brief Virtual forward function for derivative layers + * @brief Virtual update parameters function for derivative optimizers * - * @param inputs Data matrix input for forward function. + * @param layer trainable layer to have their parameters updated * * @Note It's equal to 0 because it make the derivative * object NEEDS to have these function to work. */ virtual bool update_params(trainable_layer& layer) = 0; + /** + * @brief Virtual function to update internal parameters before update_params() + * + * + */ + virtual bool pre_update_params(){ + return true; + } + + /** + * @brief Virtual function to update internal parameters after update_params() + * + * + */ + virtual bool post_update_params(){ + return true; + } + }; diff --git a/include/neural_network/optimizers/optimizer_adagrad.hpp b/include/neural_network/optimizers/optimizer_adagrad.hpp new file mode 100644 index 0000000..64d6ce1 --- /dev/null +++ b/include/neural_network/optimizers/optimizer_adagrad.hpp @@ -0,0 +1,112 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_adagrad.hpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_adagrad struct used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + + +namespace panic{ + namespace neural_network{ + + + +/** + * @brief optimizer_adagrad for the rest of the neural network library to use + * + */ +struct optimizer_adagrad: optimizer{ + + panic::types::real_t learning_rate; + + panic::types::real_t decay; + + panic::types::real_t epsilon; + + panic::types::uint_t iterations; + + + + /** + * @brief Constructor + * + * @param learning_rate The learning rate for the optimization. + * @param decay The decay for the learning rate over the interations. + * + */ + optimizer_adagrad(const panic::types::real_t learning_rate = static_cast(1), + const panic::types::real_t decay = static_cast(0), + const panic::types::real_t epsilon = static_cast(1e-7)); + + + /** + * @brief Default de-constructor + * + */ + ~optimizer_adagrad() = default; + + + /** + * @brief Updates weights and biases in trainable layers + * + * @param layer Trianable layer to update. + * + */ + bool update_params(trainable_layer& layer) override; + + + + /** + * @brief function to update internal parameters before update_params() + * + * @param layer Trianable layer to update. + * + */ + bool pre_update_params() override; + + /** + * @brief function to update internal parameters after update_params() + * + */ + bool post_update_params() override; + +}; + + } // namespace tensor +} // namespace panic + + + diff --git a/include/neural_network/optimizers/optimizer_adam.hpp b/include/neural_network/optimizers/optimizer_adam.hpp new file mode 100644 index 0000000..0e8043d --- /dev/null +++ b/include/neural_network/optimizers/optimizer_adam.hpp @@ -0,0 +1,120 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_rmsprop.hpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_adam struct used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + + +namespace panic{ + namespace neural_network{ + + + +/** + * @brief optimizer_adam for the rest of the neural network library to use + * + */ +struct optimizer_adam: optimizer{ + + panic::types::real_t learning_rate; + + panic::types::real_t decay; + + panic::types::real_t epsilon; + + panic::types::real_t beta_1; + panic::types::real_t beta_2; + + panic::types::real_t beta_1_power; + panic::types::real_t beta_2_power; + + panic::types::uint_t iterations; + + + + /** + * @brief Constructor + * + * @param learning_rate The learning rate for the optimization. + * @param decay The decay for the learning rate over the interations. + * + */ + optimizer_adam(const panic::types::real_t learning_rate = static_cast(0.001), + const panic::types::real_t decay = static_cast(0), + const panic::types::real_t epsilon = static_cast(1e-7), + const panic::types::real_t beta_1 = static_cast(0.9), + const panic::types::real_t beta_2 = static_cast(0.999)); + + + /** + * @brief Default de-constructor + * + */ + ~optimizer_adam() = default; + + + /** + * @brief Updates weights and biases in trainable layers + * + * @param layer Trianable layer to update. + * + */ + bool update_params(trainable_layer& layer) override; + + + + /** + * @brief function to update internal parameters before update_params() + * + * @param layer Trianable layer to update. + * + */ + bool pre_update_params() override; + + /** + * @brief function to update internal parameters after update_params() + * + */ + bool post_update_params() override; + +}; + + } // namespace tensor +} // namespace panic + + + diff --git a/include/neural_network/optimizers/optimizer_rmsprop.hpp b/include/neural_network/optimizers/optimizer_rmsprop.hpp new file mode 100644 index 0000000..fd673f1 --- /dev/null +++ b/include/neural_network/optimizers/optimizer_rmsprop.hpp @@ -0,0 +1,115 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_rmsprop.hpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_rmsprop struct used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + + +namespace panic{ + namespace neural_network{ + + + +/** + * @brief optimizer_rmsprop for the rest of the neural network library to use + * + */ +struct optimizer_rmsprop: optimizer{ + + panic::types::real_t learning_rate; + + panic::types::real_t decay; + + panic::types::real_t epsilon; + + panic::types::real_t rho; + + panic::types::uint_t iterations; + + + + /** + * @brief Constructor + * + * @param learning_rate The learning rate for the optimization. + * @param decay The decay for the learning rate over the interations. + * + */ + optimizer_rmsprop(const panic::types::real_t learning_rate = static_cast(0.001), + const panic::types::real_t decay = static_cast(0), + const panic::types::real_t epsilon = static_cast(1e-7), + const panic::types::real_t rho = static_cast(0.9)); + + + /** + * @brief Default de-constructor + * + */ + ~optimizer_rmsprop() = default; + + + /** + * @brief Updates weights and biases in trainable layers + * + * @param layer Trianable layer to update. + * + */ + bool update_params(trainable_layer& layer) override; + + + + /** + * @brief function to update internal parameters before update_params() + * + * @param layer Trianable layer to update. + * + */ + bool pre_update_params() override; + + /** + * @brief function to update internal parameters after update_params() + * + */ + bool post_update_params() override; + +}; + + } // namespace tensor +} // namespace panic + + + diff --git a/include/neural_network/optimizers/optimizer_sgd.hpp b/include/neural_network/optimizers/optimizer_sgd.hpp index 3459195..2e61af6 100644 --- a/include/neural_network/optimizers/optimizer_sgd.hpp +++ b/include/neural_network/optimizers/optimizer_sgd.hpp @@ -52,13 +52,24 @@ struct optimizer_sgd: optimizer{ panic::types::real_t learning_rate; + panic::types::real_t decay; + + panic::types::real_t momentum; + + panic::types::uint_t iterations; + + + /** * @brief Constructor * * @param learning_rate The learning rate for the optimization. + * @param decay The decay for the learning rate over the interations. * */ - optimizer_sgd(const panic::types::real_t learning_rate = static_cast(1e-3)); + optimizer_sgd(const panic::types::real_t learning_rate = static_cast(1), + const panic::types::real_t decay = static_cast(0), + const panic::types::real_t momentum = static_cast(0)); /** @@ -76,10 +87,24 @@ struct optimizer_sgd: optimizer{ */ bool update_params(trainable_layer& layer) override; + + + /** + * @brief function to update internal parameters before update_params() + * + * @param layer Trianable layer to update. + * + */ + bool pre_update_params() override; + + /** + * @brief function to update internal parameters after update_params() + * + */ + bool post_update_params() override; + }; - - } // namespace tensor } // namespace panic diff --git a/main.cpp b/main.cpp index b6de104..d32e5af 100644 --- a/main.cpp +++ b/main.cpp @@ -64,6 +64,7 @@ #include #include #include +#include #include @@ -709,8 +710,9 @@ int main(void) { panic::types::uint_t classes = 3; - // create spiral data + // create spiral data panic::neural_network::spiral_data(samples, classes, X, y); + //panic::neural_network::vertical_data(samples, classes, X, y); // Initilise my model panic::neural_network::model mymodel; @@ -741,10 +743,42 @@ int main(void) { return false; } - if (! mymodel.add_optimizer_sgd()){ + /* + panic::types::real_t learning_rate = 1; + panic::types::real_t decay = 1e-3; + panic::types::real_t momentum = 0.9; + if (! mymodel.add_optimizer_sgd(learning_rate, decay, momentum)){ return false; } - + */ + /* + panic::types::real_t learning_rate = 1; + panic::types::real_t decay = 1e-4; + panic::types::real_t epsilon = 1e-7; + if (! mymodel.add_optimizer_adagrad(learning_rate, decay, epsilon)){ + return false; + } + */ + /* + panic::types::real_t learning_rate = 0.001; + panic::types::real_t decay = 1e-4; + panic::types::real_t epsilon = 1e-7; + panic::types::real_t rho = 0.999; + if (! mymodel.add_optimizer_rmsprop(learning_rate, decay, epsilon, rho)){ + return false; + } + */ + + panic::types::real_t learning_rate = 0.02; + panic::types::real_t decay = 1e-5; + panic::types::real_t epsilon = 1e-7; + panic::types::real_t beta_1 = 0.9; + panic::types::real_t beta_2 = 0.999; + + if (! mymodel.add_optimizer_adam(learning_rate, decay, epsilon, beta_1, beta_2)){ + return false; + } + if (!mymodel.finalize()){ return false; } @@ -752,15 +786,13 @@ int main(void) { panic::types::uint_t epochs = 10000; - panic::types::uint_t print_every = 100; + panic::types::uint_t print_every = 250; if (!mymodel.train(X, y, epochs, print_every)){ std::cout << "Training failed" << std::endl; return false; } - - return 0; } \ No newline at end of file diff --git a/src/math/sqrt.cpp b/src/math/sqrt.cpp new file mode 100644 index 0000000..e381ecb --- /dev/null +++ b/src/math/sqrt.cpp @@ -0,0 +1,322 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: math + * File Name: sqrt.cpp + * Revision: 0.1.0 + * Date: 06-08-2026 + * Author: Michelle Bausager + * + * Description: + * Functions to calculate the sqrt of numbers + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//----------------------------------------------------------------------------------------------------- + +#include +#include +#include + +#include // for panic::vector +#include // for panic::matrix + +//--------------------------------------------------------------------------------------------------------------------------- +// PRIVATE CONSTANTS +//--------------------------------------------------------------------------------------------------------------------------- +/** + * @brief Minimum number of element operations before using the OpenMP-enabled loop. + * + * Small vectors and matrices are kept serial because the overhead of starting + * worker threads can be larger than the work itself. + */ +static const panic::types::uint_t sqrt_omp_min_work = 250; +//--------------------------------------------------------------------------------------------------------------------------- +// INPLEMENTATION +//--------------------------------------------------------------------------------------------------------------------------- + + +namespace panic { + namespace math { +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// Calculates the square root of a non-negative real number +// using Heron's method. +// +// Returns 0 when: +// - input is negative +// - input is NaN +// - input is positive infinity +// - the method does not converge within the iteration limit +//-------------------------------------------------------------------------------------------------------------------------- +template +bool sqrt(const T x, T& y){ + + const T zero = static_cast(0); + + const T one = static_cast(1); + + const T half = static_cast(0.5); + + const T tolerance = static_cast(1e-6); + + const panic::types::uint_t max_iterations = static_cast(128); + + // Reject negative numbers and NaN. + // NaN >= 0 evaluates to false. + if (!(x >= zero)){ + return false; + } + + // sqrt(0) is exactly 0. + // This special case also avoids division by zero. + if (x == zero){ + y = T{0}; + return true; + } + + // Choose an initial estimate that works reasonably for + // values both above and below 1. + T current; + + if (x >= one){ + current = x; + } + else { + current = one; + } + + for (panic::types::uint_t iteration = static_cast(0); iteration < max_iterations; ++iteration){ + + const T next = half *(current + x / current); + + T difference = next - current; + + if (difference < zero){ + difference = - difference; + } + + /* + * next == current means floating-point rounding has + * prevented any further improvement. + */ + if (next == current || difference <= tolerance * next){ + y = next; + return true; + } + y = next; + current = next; + } + + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template bool sqrt(const panic::types::real_t x, + panic::types::real_t& y +); + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// Calculates the square root of a non-negative real number +// using Heron's method. +// +// Returns false when: +// - input is negative +// - input is NaN +// - input is positive infinity +// - the method does not converge within the iteration limit +//-------------------------------------------------------------------------------------------------------------------------- +template +T sqrt(const T x){ + T y; + if (! sqrt(x,y)){ + return T(); + } + + return y; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template panic::types::real_t sqrt(const panic::types::real_t x +); + + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// Calculates the sqrt elementwise of a vector +//-------------------------------------------------------------------------------------------------------------------------- +template +bool sqrt(const panic::tensor::vector& a, panic::tensor::vector& c){ + + if (!c.resize(a.size())){ + return false; + } + + bool valid = true; + + // valid remains true only if every sqrt is valid. + PANIC_OMP_PARALLEL_FOR_REDUCTION_IF( a.size() > sqrt_omp_min_work, &&, valid ) + for (panic::types::uint_t i = 0; i < a.size(); ++i){ + const bool nonzero = sqrt(a[i], c[i]); + valid = valid && nonzero; + } + + + return valid; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template bool sqrt(const panic::tensor::vector& a, + panic::tensor::vector& c +); + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// Calculates the sqrt elementwise for a vector +//-------------------------------------------------------------------------------------------------------------------------- +template +panic::tensor::vector sqrt(const panic::tensor::vector& a){ + panic::tensor::vector c(a.size()); + + if (!sqrt(a, c)){ + return panic::tensor::vector(); + } + + return c; +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template panic::tensor::vector + sqrt(const panic::tensor::vector& a +); + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// calculates the natrual sqrt elementwise of a matrix +//-------------------------------------------------------------------------------------------------------------------------- +template +bool sqrt(const panic::tensor::matrix& A, panic::tensor::matrix& C){ + + panic::types::uint_t rows = A.rows(); + panic::types::uint_t cols = A.cols(); + panic::types::uint_t work = rows*cols; + + if ( !C.resize(rows, cols) ){ + return false; + } + + bool valid = true; + + // valid remains true only if every sqrt is valid. + PANIC_OMP_PARALLEL_FOR_IF(work > sqrt_omp_min_work) + for (panic::types::uint_t i = 0; i < rows; ++i){ + for (panic::types::uint_t j = 0; j < cols; ++j){ + const bool nonzero = sqrt(A(i,j), C(i,j)); + valid = valid && nonzero; + } + } + + return valid; +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template bool sqrt(const panic::tensor::matrix& A, + panic::tensor::matrix& C +); + + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::math::sqrt +// +// Description: +// Calculates the natrual sqrt element-wise of a matrix +//-------------------------------------------------------------------------------------------------------------------------- +template +panic::tensor::matrix sqrt(const panic::tensor::matrix& A){ + panic::tensor::matrix C; + + if (!sqrt(A, C)){ + return panic::tensor::matrix(); + } + + return C; +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template panic::tensor::matrix + sqrt(const panic::tensor::matrix& A +); + + + } // namespace math +} // namespace panic diff --git a/src/math/sum.cpp b/src/math/sum.cpp index 9718517..d552e00 100644 --- a/src/math/sum.cpp +++ b/src/math/sum.cpp @@ -264,7 +264,7 @@ bool sum_colwise(const panic::tensor::matrix& A, panic::tensor::vector& b) for (panic::types::uint_t i = 0; i < cols; ++i){ b[i] = T{0}; for (panic::types::uint_t j = 0; j < rows; ++j){ - b[i] += A(i,j); + b[i] += A(j,i); } } diff --git a/src/neural_network/activation/activation_relu.cpp b/src/neural_network/activation/activation_relu.cpp index c986bfb..da9be88 100644 --- a/src/neural_network/activation/activation_relu.cpp +++ b/src/neural_network/activation/activation_relu.cpp @@ -51,9 +51,10 @@ * Small vectors and matrices are kept serial because the overhead of starting * worker threads can be larger than the work itself. */ -static const panic::types::uint_t activation_ReLU_omp_min_size = 500; +static const panic::types::uint_t activation_ReLU_omp_min_size = 250; static const panic::types::real_t almost_zero = static_cast (1e-7); +static const panic::types::real_t zero = static_cast (0); //--------------------------------------------------------------------------------------------------------------------------- // INPLEMENTATION //--------------------------------------------------------------------------------------------------------------------------- @@ -82,7 +83,7 @@ bool activation_relu::forward(const panic::tensor::real_matrix& input_data){ inputs = input_data; - panic::math::clip_lower(inputs, almost_zero, outputs); + panic::math::clip_lower(inputs, zero, outputs); return true; @@ -98,8 +99,28 @@ bool activation_relu::forward(const panic::tensor::real_matrix& input_data){ //-------------------------------------------------------------------------------------------------------------------------- bool activation_relu::backward(const panic::tensor::real_matrix& dvalues){ - // Zero gradients where input values were negative - dinputs = panic::math::clip_lower(dvalues, almost_zero); + if (dvalues.rows() != inputs.rows() || dvalues.cols() != inputs.cols()){ + return false; + } + + dinputs = dvalues; + + const panic::types::uint_t rows = dinputs.rows(); + const panic::types::uint_t cols = dinputs.cols(); + const panic::types::uint_t work = rows * cols; + + PANIC_OMP_PARALLEL_FOR_IF(work > activation_ReLU_omp_min_size) + for (panic::types::uint_t i = 0; i < rows; ++i){ + for (panic::types::uint_t j = 0; j < cols; ++j){ + + // ReLU output did not depend on this input, + // so no gradient passes backward. + if (inputs(i, j) <= static_cast(0)){ + + dinputs(i, j) = static_cast(0); + } + } + } return true; diff --git a/src/neural_network/datasets/spiral_data.cpp b/src/neural_network/datasets/spiral_data.cpp index 542fc8a..d142a40 100644 --- a/src/neural_network/datasets/spiral_data.cpp +++ b/src/neural_network/datasets/spiral_data.cpp @@ -96,7 +96,7 @@ bool spiral_data(const panic::types::uint_t samples, const panic::types::uint_t } panic::tensor::matrix random_matrix(samples*classes, 2); - panic::random::uniform(random_matrix, T{-0.15}, T{0.15}); + panic::random::uniform(random_matrix, T{-0.0015}, T{0.0015}); if (!panic::math::add(X, random_matrix, X)){ return false; diff --git a/src/neural_network/layer/layer_dense.cpp b/src/neural_network/layer/layer_dense.cpp index f2f3bb4..0c30afa 100644 --- a/src/neural_network/layer/layer_dense.cpp +++ b/src/neural_network/layer/layer_dense.cpp @@ -85,7 +85,7 @@ layer_dense::layer_dense() { //-------------------------------------------------------------------------------------------------------------------------- layer_dense::layer_dense(panic::types::uint_t input_size, panic::types::uint_t neurons) { weights.resize(input_size, neurons); - panic::random::uniform(weights); + panic::random::uniform(weights, static_cast(-1), static_cast(1)); panic::math::mul(weights, static_cast(0.01), weights); biases.resize(neurons); diff --git a/src/neural_network/model/model.cpp b/src/neural_network/model/model.cpp index c57cc33..31dc973 100644 --- a/src/neural_network/model/model.cpp +++ b/src/neural_network/model/model.cpp @@ -49,7 +49,10 @@ #include -#include +#include +#include +#include +#include #include #include @@ -545,9 +548,86 @@ bool model::add_loss_categorical_crossentropy(){ // Example: // model.add_optimizer_sgd(1e-4); //-------------------------------------------------------------------------------------------------------------------------- -bool model::add_optimizer_sgd(const panic::types::real_t learning_rate){ +bool model::add_optimizer_sgd(const panic::types::real_t learning_rate, const panic::types::real_t decay, const panic::types::real_t momentum){ - optimizer_sgd* new_optimizer = new optimizer_sgd(learning_rate); + optimizer_sgd* new_optimizer = new optimizer_sgd(learning_rate, decay, momentum); + + if (new_optimizer == 0){ + return false; + } + + delete optimizer_function; + optimizer_function = new_optimizer; + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::neural_network::model::add_optimizer_adagrad +// +// Description: +// Adds Stocastient Gradient Decent as an optimizer to the model. +// +// Example: +// model.add_optimizer_adagrad(1e-4); +//-------------------------------------------------------------------------------------------------------------------------- +bool model::add_optimizer_adagrad(const panic::types::real_t learning_rate, const panic::types::real_t decay, const panic::types::real_t epsilon){ + + optimizer_adagrad* new_optimizer = new optimizer_adagrad(learning_rate, decay, epsilon); + + if (new_optimizer == 0){ + return false; + } + + delete optimizer_function; + optimizer_function = new_optimizer; + + return true; +} + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::neural_network::model::add_optimizer_rmsprop +// +// Description: +// Adds Stocastient Gradient Decent as an optimizer to the model. +// +// Example: +// model.add_optimizer_rmsprop(1e-4); +//-------------------------------------------------------------------------------------------------------------------------- +bool model::add_optimizer_rmsprop(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t epsilon, + const panic::types::real_t rho){ + + optimizer_rmsprop* new_optimizer = new optimizer_rmsprop(learning_rate, decay, epsilon, rho); + + if (new_optimizer == 0){ + return false; + } + + delete optimizer_function; + optimizer_function = new_optimizer; + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::neural_network::model::add_optimizer_adam +// +// Description: +// Adds Stocastient Gradient Decent as an optimizer to the model. +// +// Example: +// model.add_optimizer_adam(1e-4); +//-------------------------------------------------------------------------------------------------------------------------- +bool model::add_optimizer_adam(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t epsilon, + const panic::types::real_t beta_1, + const panic::types::real_t beta_2){ + + optimizer_adam* new_optimizer = new optimizer_adam(learning_rate, decay, epsilon, beta_1, beta_2); if (new_optimizer == 0){ return false; @@ -618,7 +698,9 @@ bool model::optimize(){ return false; } - + if (! optimizer_function->pre_update_params()){ + return false; + } for(panic::types::uint_t i = 0; i < trainable_layer_count; ++i){ @@ -627,10 +709,16 @@ bool model::optimize(){ return false; } + + if (! optimizer_function->update_params(*trainable_layers[i]) ){ return false; } + } + + if (! optimizer_function->post_update_params()){ + return false; } @@ -654,7 +742,7 @@ bool model::train(const panic::tensor::real_matrix& X_train, panic::types::real_t accuracy; panic::tensor::uint_vector comparisons; - for (panic::types::uint_t epoch = 0; epoch < epochs; ++epoch){ + for (panic::types::uint_t epoch = 0; epoch < epochs+1; ++epoch){ forward(X_train); @@ -688,8 +776,9 @@ bool model::train(const panic::tensor::real_matrix& X_train, // Its mean is therefore the classification accuracy. accuracy = panic::math::mean(comparisons); - if (epoch % print_every == static_cast(0)){ + if ((epoch) % print_every == static_cast(0) ){ std::cout << "Epoch: " << epoch; + std::cout << " lr: " << optimizer_function->current_learning_rate; std::cout << " data loss: " << loss_function->data_loss; std::cout << " acc: " << accuracy << std::endl; } diff --git a/src/neural_network/optimizers/optimizer_adagrad.cpp b/src/neural_network/optimizers/optimizer_adagrad.cpp new file mode 100644 index 0000000..2406ea2 --- /dev/null +++ b/src/neural_network/optimizers/optimizer_adagrad.cpp @@ -0,0 +1,217 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_adagrad.cpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_adagrad used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + +#include +#include +#include +#include +#include + + +//--------------------------------------------------------------------------------------------------------------------------- +// PRIVATE CONSTANTS +//--------------------------------------------------------------------------------------------------------------------------- +/** + * @brief Minimum number of element operations before using the OpenMP-enabled loop. + * + * Small vectors and matrices are kept serial because the overhead of starting + * worker threads can be larger than the work itself. + */ +static const panic::types::uint_t optimizer_adagrad_omp_min_size = 250; + +static const panic::types::real_t real_t_1 = static_cast(1); +//--------------------------------------------------------------------------------------------------------------------------- +// INPLEMENTATION +//--------------------------------------------------------------------------------------------------------------------------- + +namespace panic{ + namespace neural_network{ + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::optimizer_adagrad +// +// Description: +// Constructor for optimizer_adagrad. +//-------------------------------------------------------------------------------------------------------------------------- +optimizer_adagrad::optimizer_adagrad(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t epsilon){ + + this->learning_rate = learning_rate; + this->decay = decay; + this->epsilon = epsilon; + + iterations = static_cast(0); + current_learning_rate = learning_rate; + + +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::update_params +// +// Description: +// Updates weights and biases in layer. +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adagrad::update_params(trainable_layer& layer) { + + + // Gradients must match their corresponding parameters. + if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){ + return false; + } + + if (layer.weights.rows() == static_cast(0) || layer.weights.cols() == static_cast(0) || layer.biases.size() == static_cast(0) ){ + return false; + } + + + // If layer doesn't contain cache arrays, create them + // and fill them with zeros + if(layer.weights.rows() != layer.weight_cache.rows() || + layer.weights.cols() != layer.weight_cache.cols() || + layer.biases.size() != layer.bias_cache.size()){ + + if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){ + return false; + } + layer.weight_cache.fill(0); + + if (! layer.bias_cache.resize(layer.biases.size())){ + return false; + } + layer.bias_cache.fill(0); + } + + // Update cache with squared current gradients + if(!panic::math::add( + layer.bias_cache, + panic::math::mul(layer.dbiases, layer.dbiases), + layer.bias_cache)){ + return false; + } + + + if(!panic::math::add( + layer.weight_cache, + panic::math::mul(layer.dweights, layer.dweights), + layer.weight_cache)){ + return false; + } + + + + panic::tensor::real_matrix weight_updates; + panic::tensor::real_vector bias_updates; + + // Vanilla SGD parameter update + normalization + // with square rooted cache + + // weight_updates = -current_learning_rate * dweights + if (!panic::math::mul(layer.dweights, -current_learning_rate, weight_updates)){ + return false; + } + // weight_updates = weight_updates / (sqrt(weight_cache) + epsilon) + if (!panic::math::div(weight_updates, panic::math::add(panic::math::sqrt(layer.weight_cache), epsilon), weight_updates)){ + return false; + } + // Add weight updates to layer + if (! panic::math::add(layer.weights, weight_updates, layer.weights)){ + return false; + } + + + // bias_updates = -current_learning_rate * dbiases + if (!panic::math::mul(layer.dbiases, -current_learning_rate, bias_updates)){ + return false; + } + // bias_updates = bias_updates / (sqrt(bias_cache) + epsilon) + if (!panic::math::div(bias_updates, panic::math::add(panic::math::sqrt(layer.bias_cache), epsilon), bias_updates)){ + return false; + } + // Add bias updates to layer + if (! panic::math::add(layer.biases, bias_updates, layer.biases)){ + return false; + } + + + + return true; + +} + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::pre_update_params +// +// Description: +// function to update internal parameters before update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adagrad::pre_update_params(){ + + if (decay){ + current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations))); + } + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::post_update_params +// +// Description: +// function to update internal parameters after update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adagrad::post_update_params(){ + + iterations += static_cast(1); + + return true; +} + + + + + + } // namespace tensor +} // namespace panic diff --git a/src/neural_network/optimizers/optimizer_adam.cpp b/src/neural_network/optimizers/optimizer_adam.cpp new file mode 100644 index 0000000..10e2d25 --- /dev/null +++ b/src/neural_network/optimizers/optimizer_adam.cpp @@ -0,0 +1,380 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_adam.cpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_adam used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + +#include +#include +#include +#include +#include + + +//--------------------------------------------------------------------------------------------------------------------------- +// PRIVATE CONSTANTS +//--------------------------------------------------------------------------------------------------------------------------- +/** + * @brief Minimum number of element operations before using the OpenMP-enabled loop. + * + * Small vectors and matrices are kept serial because the overhead of starting + * worker threads can be larger than the work itself. + */ +static const panic::types::uint_t optimizer_adam_omp_min_size = 250; + +static const panic::types::real_t real_t_1 = static_cast(1); +static const panic::types::real_t zero = static_cast(0); +//--------------------------------------------------------------------------------------------------------------------------- +// INPLEMENTATION +//--------------------------------------------------------------------------------------------------------------------------- + +namespace panic{ + namespace neural_network{ + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::optimizer_adam +// +// Description: +// Constructor for optimizer_adam. +//-------------------------------------------------------------------------------------------------------------------------- +optimizer_adam::optimizer_adam(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t epsilon, + const panic::types::real_t beta_1, + const panic::types::real_t beta_2){ + + this->learning_rate = learning_rate; + this->decay = decay; + this->epsilon = epsilon; + this->beta_1 = beta_1; + this->beta_2 = beta_2; + + this->beta_1_power = beta_1; + this->beta_2_power = beta_2; + + iterations = static_cast(0); + current_learning_rate = learning_rate; + + +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::update_params +// +// Description: +// Updates weights and biases in layer. +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adam::update_params(trainable_layer& layer) { + + + // Gradients must match their corresponding parameters. + if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){ + return false; + } + + if (layer.weights.rows() == static_cast(0) || layer.weights.cols() == static_cast(0) || layer.biases.size() == static_cast(0) ){ + return false; + } + + + // If layer doesn't contain cache arrays, create them + // and fill them with zeros + if(layer.weights.rows() != layer.weight_cache.rows() || + layer.weights.cols() != layer.weight_cache.cols() || + layer.biases.size() != layer.bias_cache.size()){ + + if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){ + return false; + } + layer.weight_cache.fill(0); + + if (! layer.bias_cache.resize(layer.biases.size())){ + return false; + } + layer.bias_cache.fill(0); + } + + // If layer doesn't contain momentum arrays, create them + // and fill them with zeros + if(layer.weights.rows() != layer.weight_momentums.rows() || + layer.weights.cols() != layer.weight_momentums.cols() || + layer.biases.size() != layer.bias_momentums.size()){ + + if (! layer.weight_momentums.resize(layer.weights.rows(), layer.weights.cols())){ + return false; + } + layer.weight_momentums.fill(0); + + if (! layer.bias_momentums.resize(layer.biases.size())){ + return false; + } + layer.bias_momentums.fill(0); + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Update momentum calculations + //-------------------------------------------------------------------------------------------------------------------------- + const panic::types::real_t one = static_cast(1); + const panic::types::real_t one_minus_beta_1 = one - beta_1; + + panic::tensor::real_matrix retained_weight_momentum; + panic::tensor::real_matrix new_weight_gradient; + + + if (!panic::math::mul(layer.weight_momentums, beta_1, retained_weight_momentum)){ + return false; + } + + if (!panic::math::mul(layer.dweights, one_minus_beta_1, new_weight_gradient)){ + return false; + } + + if (!panic::math::add(retained_weight_momentum, new_weight_gradient, layer.weight_momentums)){ + return false; + } + + + panic::tensor::real_vector retained_bias_momentum; + panic::tensor::real_vector new_bias_gradient; + + if (!panic::math::mul(layer.bias_momentums, beta_1, retained_bias_momentum)){ + return false; + } + + if (!panic::math::mul(layer.dbiases, one_minus_beta_1, new_bias_gradient)){ + return false; + } + + if (!panic::math::add(retained_bias_momentum, new_bias_gradient, layer.bias_momentums)){ + return false; + } + + + //-------------------------------------------------------------------------------------------------------------------------- + // Update momentum corrected + //-------------------------------------------------------------------------------------------------------------------------- + + const panic::types::real_t momentum_correction = one - beta_1_power; + + if (momentum_correction == zero){ + return false; + } + + panic::tensor::real_matrix corrected_weight_momentums; + panic::tensor::real_vector corrected_bias_momentums; + + if (!panic::math::div(layer.weight_momentums, momentum_correction, corrected_weight_momentums)){ + return false; + } + + if (!panic::math::div(layer.bias_momentums, momentum_correction, corrected_bias_momentums)){ + return false; + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Update cache calculations + //-------------------------------------------------------------------------------------------------------------------------- + const panic::types::real_t one_minus_beta_2 = one - beta_2; + + panic::tensor::real_matrix retained_weight_cache; + panic::tensor::real_matrix squared_dweights; + panic::tensor::real_matrix new_weight_cache_values; + + if (!panic::math::mul(layer.weight_cache, beta_2, retained_weight_cache)){ + return false; + } + + if (!panic::math::mul(layer.dweights, layer.dweights, squared_dweights)){ + return false; + } + + if (!panic::math::mul(squared_dweights, one_minus_beta_2, new_weight_cache_values)){ + return false; + } + + if (!panic::math::add(retained_weight_cache, new_weight_cache_values, layer.weight_cache)){ + return false; + } + + + panic::tensor::real_vector retained_bias_cache; + panic::tensor::real_vector squared_dbiases; + panic::tensor::real_vector new_bias_cache_values; + + if (!panic::math::mul(layer.bias_cache, beta_2, retained_bias_cache)){ + return false; + } + + if (!panic::math::mul(layer.dbiases, layer.dbiases, squared_dbiases)){ + return false; + } + + if (!panic::math::mul(squared_dbiases, one_minus_beta_2, new_bias_cache_values)){ + return false; + } + + if (!panic::math::add(retained_bias_cache, new_bias_cache_values, layer.bias_cache)){ + return false; + } + + + + //-------------------------------------------------------------------------------------------------------------------------- + // Update cache corrected + //-------------------------------------------------------------------------------------------------------------------------- + + const panic::types::real_t cache_correction = one - beta_2_power; + + if (cache_correction == zero){ + return false; + } + + panic::tensor::real_matrix corrected_weight_cache; + panic::tensor::real_vector corrected_bias_cache; + + if (!panic::math::div(layer.weight_cache, cache_correction, corrected_weight_cache)){ + return false; + } + + if (!panic::math::div(layer.bias_cache, cache_correction, corrected_bias_cache)){ + return false; + } + + + //-------------------------------------------------------------------------------------------------------------------------- + // Calculate Adam update + //-------------------------------------------------------------------------------------------------------------------------- + + + + panic::tensor::real_matrix sqrt_corrected_weight_cache; + panic::tensor::real_matrix weight_denominator; + panic::tensor::real_matrix scaled_weight_momentums; + panic::tensor::real_matrix weight_updates; + + if (!panic::math::sqrt(corrected_weight_cache, sqrt_corrected_weight_cache)){ + return false; + } + + if (!panic::math::add(sqrt_corrected_weight_cache, epsilon, weight_denominator)){ + return false; + } + + if (!panic::math::mul(corrected_weight_momentums, -current_learning_rate, scaled_weight_momentums)){ + return false; + } + + if (!panic::math::div(scaled_weight_momentums, weight_denominator, weight_updates)){ + return false; + } + + if (!panic::math::add(layer.weights, weight_updates, layer.weights)){ + return false; + } + + + panic::tensor::real_vector sqrt_corrected_bias_cache; + panic::tensor::real_vector bias_denominator; + panic::tensor::real_vector scaled_bias_momentums; + panic::tensor::real_vector bias_updates; + + if (!panic::math::sqrt(corrected_bias_cache, sqrt_corrected_bias_cache)){ + return false; + } + + if (!panic::math::add(sqrt_corrected_bias_cache, epsilon, bias_denominator)){ + return false; + } + + if (!panic::math::mul(corrected_bias_momentums, -current_learning_rate, scaled_bias_momentums )){ + return false; + } + + if (!panic::math::div(scaled_bias_momentums, bias_denominator, bias_updates)){ + return false; + } + + if (!panic::math::add(layer.biases, bias_updates, layer.biases)){ + return false; + } + + + return true; + +} + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::pre_update_params +// +// Description: +// function to update internal parameters before update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adam::pre_update_params(){ + + if (decay){ + current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations))); + } + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::post_update_params +// +// Description: +// function to update internal parameters after update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_adam::post_update_params(){ + + iterations += static_cast(1); + + beta_1_power *= beta_1; + beta_2_power *= beta_2; + + return true; +} + + + + + + } // namespace tensor +} // namespace panic diff --git a/src/neural_network/optimizers/optimizer_rmsprop.cpp b/src/neural_network/optimizers/optimizer_rmsprop.cpp new file mode 100644 index 0000000..25093f7 --- /dev/null +++ b/src/neural_network/optimizers/optimizer_rmsprop.cpp @@ -0,0 +1,349 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC 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 LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: neural_network + * File Name: optimizer_rmsprop.cpp + * Revision: 0.1.0 + * Date: 04-08-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the optimizer_rmsprop used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + +#include +#include +#include +#include +#include + + +//--------------------------------------------------------------------------------------------------------------------------- +// PRIVATE CONSTANTS +//--------------------------------------------------------------------------------------------------------------------------- +/** + * @brief Minimum number of element operations before using the OpenMP-enabled loop. + * + * Small vectors and matrices are kept serial because the overhead of starting + * worker threads can be larger than the work itself. + */ +static const panic::types::uint_t optimizer_rmsprop_omp_min_size = 250; + +static const panic::types::real_t real_t_1 = static_cast(1); +//--------------------------------------------------------------------------------------------------------------------------- +// INPLEMENTATION +//--------------------------------------------------------------------------------------------------------------------------- + +namespace panic{ + namespace neural_network{ + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::optimizer_rmsprop +// +// Description: +// Constructor for optimizer_rmsprop. +//-------------------------------------------------------------------------------------------------------------------------- +optimizer_rmsprop::optimizer_rmsprop(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t epsilon, + const panic::types::real_t rho){ + + this->learning_rate = learning_rate; + this->decay = decay; + this->epsilon = epsilon; + this->rho = rho; + + iterations = static_cast(0); + current_learning_rate = learning_rate; + + +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::update_params +// +// Description: +// Updates weights and biases in layer. +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_rmsprop::update_params(trainable_layer& layer) { + + + // Gradients must match their corresponding parameters. + if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){ + return false; + } + + if (layer.weights.rows() == static_cast(0) || layer.weights.cols() == static_cast(0) || layer.biases.size() == static_cast(0) ){ + return false; + } + + + // If layer doesn't contain cache arrays, create them + // and fill them with zeros + if(layer.weights.rows() != layer.weight_cache.rows() || + layer.weights.cols() != layer.weight_cache.cols() || + layer.biases.size() != layer.bias_cache.size()){ + + if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){ + return false; + } + layer.weight_cache.fill(0); + + if (! layer.bias_cache.resize(layer.biases.size())){ + return false; + } + layer.bias_cache.fill(0); + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Update caches with squared current gradients + //-------------------------------------------------------------------------------------------------------------------------- + + const panic::types::real_t one_minus_rho = + static_cast(1) - rho; + + //-------------------------------------------------------------------------------------------------------------------------- + // Weight cache + //-------------------------------------------------------------------------------------------------------------------------- + + panic::tensor::real_matrix retained_weight_cache; + panic::tensor::real_matrix squared_dweights; + panic::tensor::real_matrix weighted_dweights; + + if (!panic::math::mul( + layer.weight_cache, + rho, + retained_weight_cache + )){ + return false; + } + + if (!panic::math::mul( + layer.dweights, + layer.dweights, + squared_dweights + )){ + return false; + } + + if (!panic::math::mul( + squared_dweights, + one_minus_rho, + weighted_dweights + )){ + return false; + } + + if (!panic::math::add( + retained_weight_cache, + weighted_dweights, + layer.weight_cache + )){ + return false; + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Bias cache + //-------------------------------------------------------------------------------------------------------------------------- + + panic::tensor::real_vector retained_bias_cache; + panic::tensor::real_vector squared_dbiases; + panic::tensor::real_vector weighted_dbiases; + + if (!panic::math::mul( + layer.bias_cache, + rho, + retained_bias_cache + )){ + return false; + } + + if (!panic::math::mul( + layer.dbiases, + layer.dbiases, + squared_dbiases + )){ + return false; + } + + if (!panic::math::mul( + squared_dbiases, + one_minus_rho, + weighted_dbiases + )){ + return false; + } + + if (!panic::math::add( + retained_bias_cache, + weighted_dbiases, + layer.bias_cache + )){ + return false; + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Calculate weight updates + //-------------------------------------------------------------------------------------------------------------------------- + + panic::tensor::real_matrix weight_updates; + + if (!panic::math::mul( + layer.dweights, + -current_learning_rate, + weight_updates + )){ + return false; + } + + panic::tensor::real_matrix sqrt_weight_cache; + panic::tensor::real_matrix weight_denominator; + + if (!panic::math::sqrt( + layer.weight_cache, + sqrt_weight_cache + )){ + return false; + } + + if (!panic::math::add( + sqrt_weight_cache, + epsilon, + weight_denominator + )){ + return false; + } + + if (!panic::math::div( + weight_updates, + weight_denominator, + weight_updates + )){ + return false; + } + + if (!panic::math::add( + layer.weights, + weight_updates, + layer.weights + )){ + return false; + } + + //-------------------------------------------------------------------------------------------------------------------------- + // Calculate bias updates + //-------------------------------------------------------------------------------------------------------------------------- + + panic::tensor::real_vector bias_updates; + + if (!panic::math::mul( + layer.dbiases, + -current_learning_rate, + bias_updates + )){ + return false; + } + + panic::tensor::real_vector sqrt_bias_cache; + panic::tensor::real_vector bias_denominator; + + if (!panic::math::sqrt( + layer.bias_cache, + sqrt_bias_cache + )){ + return false; + } + + if (!panic::math::add( + sqrt_bias_cache, + epsilon, + bias_denominator + )){ + return false; + } + + if (!panic::math::div( + bias_updates, + bias_denominator, + bias_updates + )){ + return false; + } + + if (!panic::math::add( + layer.biases, + bias_updates, + layer.biases + )){ + return false; + } + + return true; + +} + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::pre_update_params +// +// Description: +// function to update internal parameters before update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_rmsprop::pre_update_params(){ + + if (decay){ + current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations))); + } + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::post_update_params +// +// Description: +// function to update internal parameters after update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_rmsprop::post_update_params(){ + + iterations += static_cast(1); + + return true; +} + + + + + + } // namespace tensor +} // namespace panic diff --git a/src/neural_network/optimizers/optimizer_sgd.cpp b/src/neural_network/optimizers/optimizer_sgd.cpp index 4a14981..471a924 100644 --- a/src/neural_network/optimizers/optimizer_sgd.cpp +++ b/src/neural_network/optimizers/optimizer_sgd.cpp @@ -40,6 +40,7 @@ #include #include +#include //--------------------------------------------------------------------------------------------------------------------------- // PRIVATE CONSTANTS @@ -51,6 +52,8 @@ * worker threads can be larger than the work itself. */ static const panic::types::uint_t optimizer_sgd_omp_min_size = 250; + +static const panic::types::real_t real_t_1 = static_cast(1); //--------------------------------------------------------------------------------------------------------------------------- // INPLEMENTATION //--------------------------------------------------------------------------------------------------------------------------- @@ -65,8 +68,18 @@ namespace panic{ // Description: // Constructor for optimizer_sgd. //-------------------------------------------------------------------------------------------------------------------------- -optimizer_sgd::optimizer_sgd(const panic::types::real_t learning_rate) { +optimizer_sgd::optimizer_sgd(const panic::types::real_t learning_rate, + const panic::types::real_t decay, + const panic::types::real_t momentum) { + this->learning_rate = learning_rate; + this->decay = decay; + this->momentum = momentum; + + iterations = static_cast(0); + current_learning_rate = learning_rate; + + } //-------------------------------------------------------------------------------------------------------------------------- @@ -87,28 +100,79 @@ bool optimizer_sgd::update_params(trainable_layer& layer) { return false; } - panic::tensor::real_matrix weight_updates; panic::tensor::real_vector bias_updates; - // weight_updates = -learning_rate * dweights - if (!panic::math::mul(layer.dweights, -learning_rate, weight_updates)){ - return false; + + if(momentum){ + + // If layer doesn't contain momentum arrays, create them + // and fill them with zeros + if(layer.weights.rows() != layer.weight_momentums.rows() || + layer.weights.cols() != layer.weight_momentums.cols() || + layer.biases.size() != layer.bias_momentums.size()){ + + if (! layer.weight_momentums.resize(layer.weights.rows(), layer.weights.cols())){ + return false; + } + layer.weight_momentums.fill(0); + + if (! layer.bias_momentums.resize(layer.biases.size())){ + return false; + } + layer.bias_momentums.fill(0); + } + + // Build weight updates with momentum - take previous + // updates multiplied by retain factor and update with + // current gradient + if (!panic::math::sub( + panic::math::mul(layer.weight_momentums, momentum), + panic::math::mul(layer.dweights, current_learning_rate), + weight_updates)){ + + return false; + } + + // Build bias updates + if (!panic::math::sub( + panic::math::mul(layer.bias_momentums, momentum), + panic::math::mul(layer.dbiases, current_learning_rate), + bias_updates)){ + + return false; + } + + layer.weight_momentums = weight_updates; + layer.bias_momentums = bias_updates; + + + } + else{ + + // weight_updates = -current_learning_rate * dweights + if (!panic::math::mul(layer.dweights, -current_learning_rate, weight_updates)){ + return false; + } + + + + // bias_updates = -current_learning_rate * dbiases + if (!panic::math::mul(layer.dbiases, -current_learning_rate, bias_updates)){ + return false; + } + + } // weights += weight_updates if (!panic::math::add(layer.weights, weight_updates, layer.weights)){ - return false; - } - - // bias_updates = -learning_rate * dbiases - if (!panic::math::mul(layer.dbiases, -learning_rate, bias_updates)){ - return false; + return false; } // biases += bias_updates if (!panic::math::add(layer.biases, bias_updates, layer.biases)){ - return false; + return false; } return true; @@ -116,5 +180,39 @@ bool optimizer_sgd::update_params(trainable_layer& layer) { } + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::pre_update_params +// +// Description: +// function to update internal parameters before update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_sgd::pre_update_params(){ + + if (decay){ + current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations))); + } + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::neural_network::post_update_params +// +// Description: +// function to update internal parameters after update_params() +//-------------------------------------------------------------------------------------------------------------------------- +bool optimizer_sgd::post_update_params(){ + + iterations += static_cast(1); + + return true; +} + + + + + } // namespace tensor } // namespace panic