Non-linear algebra module.

non_linear_algebra/doc.txt

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de> 
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup    non_linear_algebra NON_LINEAR_ALGEBRA
+ *
+ * @brief       Non-linear algebra operations
+ *
+ * @details The non-linear algebra module contains functions to solve multi-variant nonlinear 
+ *          equations as wells algorithms solving problems of regression smoothing and curve 
+ *          fitting. This module also enables enables the optimization of an approximate solution
+ *          by using Non-linear Least Squares (NLS) methods such as modified Gauss--Newton (GN) or 
+ *          the Levenberg--Marquardt (LVM) algorithms.
+ *
+ * @author  Zakaria Kasmi
+ */

non_linear_algebra/optimization/doc.txt

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de> 
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup    optimization OPTIMIZATION
+ * @ingroup     non_linear_algebra
+ *
+ * @brief       Solving problems of regression smoothing and curve fitting. 
+ *          
+ * @details     This module also enables the optimization of an 
+ *              approximate solution by using Non-linear Least Squares 
+ *              (NLS) methods such as modified Gauss--Newton (GN) or 
+ *              the Levenberg--Marquardt (LVM) algorithms.
+ *
+ * @author      Zakaria Kasmi
+ */

non_linear_algebra/optimization/include/levenberg_marquardt.h

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup     optimization
+ * @{
+ *
+ * @file
+ * @brief       Implement the Levenberg--Marquardt (LVM) algorithm.
+ * @note        This function is generally implemented.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author      Naouar Guerchali
+ *
+ * @}
+ */
+
+#ifndef LEVENBERG_MARQUARDT_H_
+#define LEVENBERG_MARQUARDT_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief   Implements the Levenberg--Marquardt (LVM) algorithm.
+ * @details The user should provide pointers to the error and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] n                        length of the start vector.
+ * @param[in] x0_vec[]                 start vector.
+ * @param[in] data_vec[]               data vector.
+ * @param[in] eps                      accuracy bound.
+ * @param[in] tau                      \f$ \tau \f$ factor.
+ * @param[in] beta0                    \f$ \beta_0 \f$ factor.
+ * @param[in] beta1                    \f$ \beta_1 \f$ factor.
+ * @param[in] max_iter_num             maximal iteration number of the LVM algorithm.
+ * @param[out] est_x_vec[]             estimated (optimized) vector.
+ * @param[in] (*get_f_error)           pointer to the error function.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  required iteration number.
+ *
+ */
+uint8_t opt_levenberg_marquardt(uint8_t f_length, uint8_t n,
+                                vector_t x0_vec[n],
+                                vector_t data_vec[f_length],
+                                matrix_t eps, matrix_t tau, matrix_t beta0,
+                                matrix_t beta1,
+                                uint8_t max_iter_num,
+                                vector_t est_x_vec[n],
+                                void (*get_f_error)(vector_t x0_vec[],
+                                                    vector_t data_vec[],
+                                                    vector_t f_vec[]),
+                                void (*get_jacobian)(vector_t x0_vec[],
+                                                     matrix_t J[][n])
+                                );
+
+/**
+ * @brief   Compute the initial value \f$ \mu_0 \f$ of the Levenberg--Marquardt (LVM) algorithm.
+ * @details The user should provide a pointer to the matrix
+ *          \f$ J_f^{T} J_{f} \f$.
+ *
+ * @param[in] n          column number of the matrix \f$J_f^T J_f\f$.
+ * @param[in] tau        \f$ \tau \f$ factor.
+ * @param[in] JTJ[][]    pointer to the matrix \f$J_f^T J_f\f$.
+ *
+ * @return  parameter \f$ \rho_{\mu} \f$
+ */
+matrix_t opt_levenberg_marquardt_get_mu0(uint8_t n, matrix_t tau,
+                                         matrix_t JTJ[][n]);
+
+#endif /* LEVENBERG_MARQUARDT_H_ */

non_linear_algebra/optimization/include/modified_gauss_newton.h

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup     optimization
+ * @{
+ *
+ * @file
+ * @brief       Implement the Gauss--Newton algorithm.
+ * @note        This function is generally implemented.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author      Abdelmoumen Norrdine <a.norrdine@googlemail.com>
+ *
+ * @}
+ */
+
+#ifndef MODIFIED_GAUSS_NEWTON_H_
+#define MODIFIED_GAUSS_NEWTON_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief   Implements the modified Gauss--Newton algorithm.
+ * @details The user should provide pointers to the error and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] n                        length of the start vector.
+ * @param[in] x0_vec[]                 start vector.
+ * @param[in] data_vec[]               data vector.
+ * @param[in] eps                      accuracy bound.
+ * @param[in] fmin                     termination tolerance on the error function.
+ * @param[in] max_iter_num             maximal iteration number of the Gauss--Newton algorithm.
+ * @param[out] est_x_vec[]             estimated (optimized) vector.
+ * @param[in] (*get_f_error)           pointer to the error function.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  required iteration number.
+ *
+ */
+uint8_t modified_gauss_newton(uint8_t f_length, uint8_t n,
+                              vector_t x0_vec[n],
+                              vector_t data_vec[f_length],
+                              matrix_t eps, matrix_t fmin, uint8_t max_iter_num,
+                              vector_t est_x_vec[n],
+                              void (*get_f_error)(vector_t x0_vec[],
+                                                  vector_t data_vec[],
+                                                  vector_t f_vec[]),
+                              void (*get_jacobian)(vector_t x0_vec[],
+                                                   matrix_t J[][n])
+                              );
+
+#endif /* MODIFIED_GAUSS_NEWTON_H_ */

non_linear_algebra/optimization/levenberg_marquardt.c

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief       Implement the Levenberg--Marquardt (LVM) algorithm.
+ * @note        This function is generally implemented.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author      Naouar Guerchali
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+#include <inttypes.h>
+
+#include "levenberg_marquardt.h"
+#include "matrix.h"
+#include "vector.h"
+#include "solve.h"
+#include "utils.h"
+
+/**
+ * @brief    Implements the correction-function of the Levenberg--Marquardt (LVM) algorithm.
+ * @details  The user should provide pointers to the error and Jacobian functions.
+ *
+ * @note     This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] n                        length of the start vector.
+ * @param[in] x_vec[]                  start vector.
+ * @param[in] data_vec[]               data vector.
+ * @param[in] mu                       regularization parameter \f$ \mu \f$.
+ * @param[out] s[]                     correction vector.
+ * @param[in] (*get_f_error)           pointer to the error function that calculates the matrix
+ *                                     \f$ J_f^{T} J_{f} \f$.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ *
+ * @return  the parameter \f$ \rho_{\mu} \f$
+ */
+matrix_t opt_levenberg_marquardt_correction(uint8_t f_length, uint8_t n,
+                                            matrix_t x_vec[n],
+                                            matrix_t data_vec[f_length],
+                                            matrix_t mu,
+                                            matrix_t s[n],
+                                            void (*get_f_error)(
+                                                vector_t x_vec[],
+                                                vector_t data_vec[],
+                                                vector_t f_vec[]),
+                                            void (*get_jacobian)(
+                                                vector_t x_vec[],
+                                                matrix_t J[][n])
+                                            )
+{
+
+    matrix_t Fx[f_length];
+    matrix_t J[f_length][n];
+    vector_t f_vec[f_length];
+    matrix_t JTJ_mu2_I[n][n];
+    vector_t JTF[n];
+    matrix_t ro_mu = 0;
+    matrix_t Fx_square = 0;
+    matrix_t Fx_plus_s_square = 0;
+    matrix_t Fx_plus_J_mul_s_square = 0;
+    vector_t Fx_plus_J_mul_s[f_length];
+    vector_t F_x_plus_s[f_length];
+    vector_t x_plus_s[n];
+    matrix_t denom;
+
+    /*
+     * compute: -(J'J + mu^2*I)
+     */
+    //JT_J = J'*J
+    get_jacobian(x_vec, J);
+
+    //store J'J in JTJ_mu2_I
+    matrix_trans_mul_itself(f_length, n, J, JTJ_mu2_I);
+
+    //compute: J'J + mu^2*I
+    matrix_add_to_diag(n, JTJ_mu2_I, n, pow(mu, 2));
+
+    // -(J'*J+mu^2*eye(n)), is an (nxn) matrix
+    matrix_mul_scalar(n, n, JTJ_mu2_I, -1, JTJ_mu2_I);
+
+    //JT_f = J'*f
+    get_f_error(x_vec, data_vec, f_vec);
+    matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JTF);
+
+    //solve the equation: s = -(J'*J+mu^2*I)\(J'*Fx);
+    solve_householder(n, n, JTJ_mu2_I, JTF, s);
+
+    //f(x)
+    get_f_error(x_vec, data_vec, Fx);
+
+    //f(x)^2
+    Fx_square = vector_get_scalar_product(f_length, Fx, Fx);
+
+    //(x+s)
+    vector_add(n, x_vec, s, x_plus_s);
+
+    //f(x+s)
+    get_f_error(x_plus_s, data_vec, F_x_plus_s);
+
+    //f(x+s)^2
+    Fx_plus_s_square = vector_get_scalar_product(f_length, F_x_plus_s,
+                                                 F_x_plus_s);
+
+    //compute: J(x)*s
+    matrix_mul_vec(f_length, n, J, s, Fx_plus_J_mul_s);
+
+    //compute: f(x) + J(x)*s
+    vector_add(f_length, Fx, Fx_plus_J_mul_s, Fx_plus_J_mul_s);
+
+    //compute: ||f(x) + J(x)*s||^2
+    Fx_plus_J_mul_s_square = vector_get_scalar_product(f_length,
+                                                       Fx_plus_J_mul_s,
+                                                       Fx_plus_J_mul_s);
+
+    denom = Fx_square - Fx_plus_J_mul_s_square;
+    if (denom != 0) {
+        ro_mu = (Fx_square - Fx_plus_s_square) / denom;
+    }
+    else {
+        puts("ro_mu is infinite !!!");
+        ro_mu = FLT_MAX;
+    }
+
+    return ro_mu;
+}
+
+uint8_t opt_levenberg_marquardt(uint8_t f_length, uint8_t n,
+                                vector_t x0_vec[n],
+                                vector_t data_vec[f_length],
+                                matrix_t eps, matrix_t tau, matrix_t beta0,
+                                matrix_t beta1,
+                                uint8_t max_iter_num,
+                                vector_t est_x_vec[n],
+                                void (*get_f_error)(vector_t x0_vec[],
+                                                    vector_t data_vec[],
+                                                    vector_t f_vec[]),
+                                void (*get_jacobian)(vector_t x0_vec[],
+                                                     matrix_t J[][n])
+                                )
+{
+    matrix_t J[f_length][n];
+    vector_t JT_f[n];
+    matrix_t JTJ_mu2_I[n][n];
+    vector_t s[n];
+    vector_t f_vec[f_length];
+    matrix_t mu;
+    matrix_t ro_mu;
+    uint8_t it;
+
+    /*
+     * compute mu0 and -(J'J + mu0^2*I)
+     */
+    //JT_J = J'*J
+    get_jacobian(x0_vec, J);
+
+    //store J'J in JTJ_mu2_I
+    matrix_trans_mul_itself(f_length, n, J, JTJ_mu2_I);
+
+    // compute mu_0: mu = mu_0
+    mu = opt_levenberg_marquardt_get_mu0(n, tau, JTJ_mu2_I);
+
+    //compute: J'J + mu0^2*I
+    matrix_add_to_diag(n, JTJ_mu2_I, n, pow(mu, 2));
+
+    //-(J'*J+mu0^2*eye(n)) is a (n,n) matrix
+    matrix_mul_scalar(n, n, JTJ_mu2_I, -1, JTJ_mu2_I);
+
+    //compute JTF: JT_f = J'*f
+    get_f_error(x0_vec, data_vec, f_vec);
+    matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JT_f);
+
+    //solve the equation: s=-(J'*J+mu^2*eye(n))\(J'*F);
+    solve_householder(n, n, JTJ_mu2_I, JT_f, s);
+
+    vector_copy(n, x0_vec, est_x_vec);
+
+    it = 0;
+    while ((vector_get_norm2(n, s)
+            > eps * (1 + vector_get_norm2(n, est_x_vec)))
+           && (it < max_iter_num)) {   //norm(s,2)>eps*(1+norm(x,2))
+
+        ro_mu = opt_levenberg_marquardt_correction(f_length, n,
+                                                   est_x_vec, data_vec,
+                                                   mu, s, get_f_error,
+                                                   get_jacobian);
+
+        while (1) {
+            if (ro_mu <= beta0) {
+                mu = 2.0 * mu;
+                ro_mu = opt_levenberg_marquardt_correction(
+                    f_length, n,
+                    est_x_vec,
+                    data_vec, mu, s, get_f_error,
+                    get_jacobian);
+            }
+            else if (ro_mu >= beta1) {
+                mu = mu / 2.0;
+                break;
+            }
+            else {
+                break;
+            }
+        }
+        vector_add(n, est_x_vec, s, est_x_vec);
+        it = it + 1;
+    } //while
+
+    return it;
+}
+
+matrix_t opt_levenberg_marquardt_get_mu0(uint8_t n, matrix_t tau,
+                                         matrix_t JTJ[][n])
+{
+    matrix_t max_diag_JTJ = JTJ[0][0];
+
+    for (uint8_t i = 1; i < n; i++) {
+        if (JTJ[i][i] > max_diag_JTJ) {
+            max_diag_JTJ = JTJ[i][i];
+        }
+    }
+
+    return tau * max_diag_JTJ;
+}

non_linear_algebra/optimization/modified_gauss_newton.c

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief       Implement the Gauss--Newton algorithm.
+ * @note        This function is generally implemented.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author      Abdelmoumen Norrdine <a.norrdine@googlemail.com>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <math.h>
+
+#include "utils.h"
+#include "matrix.h"
+#include "vector.h"
+#include "moore_penrose_pseudo_inverse.h"
+
+//n is the size of x0 that is equal to the column number of J
+uint8_t modified_gauss_newton(uint8_t f_length, uint8_t n,
+                              vector_t x0_vec[n],
+                              vector_t data_vec[f_length],
+                              matrix_t eps, matrix_t fmin, uint8_t max_iter_num,
+                              vector_t est_x_vec[n],
+                              void (*get_f_error)(vector_t x0_vec[],
+                                                  vector_t data_vec[],
+                                                  vector_t f_vec[]),
+                              void (*get_jacobian)(vector_t x0_vec[],
+                                                   matrix_t J[][n])
+                              )
+
+{
+    matrix_t J[f_length][n];
+    matrix_t JT_J[n][n];
+    matrix_t JT_f[n];
+    vector_t f_vec[f_length];
+    matrix_t f_error;
+    matrix_t pinv_JTJ_mat[n][n];
+    vector_t correction_vec[n];
+    vector_t x_vec[n];
+    vector_t next_x_vec[n];
+    matrix_t max_error, min_error;
+    matrix_t step;
+    uint8_t iter_num;
+
+    get_f_error(x0_vec, data_vec, f_vec);
+    f_error = vector_get_norm2(f_length, f_vec);
+    max_error = f_error;
+    min_error = max_error;
+    step = eps;
+
+    vector_copy(n, x0_vec, x_vec);
+    vector_copy(n, x0_vec, est_x_vec);
+    iter_num = 0;
+
+    while ((step >= eps) && (iter_num < max_iter_num) && (f_error > fmin)) {
+        /*
+         * Compute then correction terms & next x_vec values
+         */
+        //JT_J = J'*J
+        get_jacobian(x_vec, J);
+        matrix_trans_mul_itself(f_length, n, J, JT_J);
+
+        //JT_f = J'*f
+        get_f_error(x_vec, data_vec, f_vec);
+        matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JT_f);
+
+        //solve: J'J*s = -J'*f
+        moore_penrose_get_pinv(n, n, JT_J, pinv_JTJ_mat);
+        //s = (J'J)\J'*f
+        matrix_mul_vec(n, n, pinv_JTJ_mat, JT_f, correction_vec);
+
+        // x = x - s
+        vector_sub(n, x_vec, correction_vec, next_x_vec);
+
+        //next step
+        step = vector_get_euclidean_distance(n, x_vec, next_x_vec);
+
+        // x_vec = next_x_vec
+        vector_copy(n, next_x_vec, x_vec);
+
+        //error vector
+        get_f_error(x_vec, data_vec, f_vec);
+
+        f_error = vector_get_norm2(f_length, f_vec);
+
+        if (min_error > f_error) {  //store the x_vec value with the minimum error in est_x_vec
+            vector_copy(n, x_vec, est_x_vec);
+            min_error = f_error;    // update min_error
+        }
+
+        max_error = utils_max(f_error, max_error); // update max_error
+        iter_num++;
+        if ((max_error - min_error) > 10) {
+            break;
+        }
+
+    } //while
+
+    return iter_num;
+}

non_linear_algebra/solve_non_linear_equations/damped_newton_raphson.c

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief       Implement the damped Newton--Raphson algorithm.
+ * @details     The damped Newton--Raphson algorithm enables to solve
+ *              multi-variant nonlinear equation systems.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+
+#include "moore_penrose_pseudo_inverse.h"
+#include "damped_newton_raphson.h"
+#include "vector.h"
+#include "matrix.h"
+
+//Multivariate Newtons Method
+//n is the size of x0 that is equal to the column number of J
+uint8_t damped_newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+                              double min_lamda, double eps, uint8_t max_it_num,
+                              vector_t est_x_arr[],
+                              void (*get_non_lin_sys)(vector_t x_arr[],
+                                                      vector_t f_vec[]),
+                              void (*get_jacobian)(vector_t x_arr[],
+                                                   matrix_t J[][n]))
+{
+    uint8_t iter_num;
+    double lamda;
+    double lamda_c;
+    double damped_norm_x;
+    double damped_norm_prev_x;
+    vector_t prev_x_arr[n];
+    vector_t delta_x[n];
+    vector_t lamda_mul_delta_x[n];
+
+    damped_norm_x = get_damped_norm(f_length, n, x0_arr, get_non_lin_sys,
+                                    get_jacobian);
+    iter_num = 0;
+    vector_copy(n, x0_arr, prev_x_arr);
+    while ((damped_norm_x >= eps) && (iter_num < max_it_num)) {
+
+        get_delta_x(f_length, n, prev_x_arr, get_non_lin_sys,
+                    get_jacobian, delta_x);
+        lamda = 1.0;
+
+        // x = x_k +lamda*s_k (lamda=1)
+        vector_add(n, prev_x_arr, delta_x, est_x_arr);
+
+        lamda_c = 1 - lamda / 4;
+        //|||f(x)|||_k
+        damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+                                        get_non_lin_sys, get_jacobian);
+
+        //|||f(prev_x)|||_k
+        damped_norm_prev_x = vector_get_norm2(n, delta_x);
+
+        while (damped_norm_x > (lamda_c * damped_norm_prev_x)) {
+            lamda = lamda / 2;
+            if (lamda >= min_lamda) {
+                //x = x_k +lamda*s_k
+                vector_scalar_mul(n, delta_x, lamda,
+                                  lamda_mul_delta_x);
+                vector_add(n, prev_x_arr, lamda_mul_delta_x,
+                           est_x_arr);
+            }
+            else {
+                break;
+            }
+            // |||f(x)|||_k
+            damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+                                            get_non_lin_sys, get_jacobian);
+
+        } //while
+
+        if (lamda <= min_lamda) {
+            break;
+        }
+
+        //update prev_x
+        vector_copy(n, est_x_arr, prev_x_arr);
+
+        // |||f(x)|||_k
+        damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+                                        get_non_lin_sys, get_jacobian);
+        iter_num++;
+    }
+
+    return iter_num;
+}
+
+double get_damped_norm(uint8_t m, uint8_t n, vector_t x_arr[],
+                       void (*get_non_lin_sys)(vector_t x_arr[],
+                                               vector_t f_vec[]),
+                       void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n])
+                       )
+{
+    double damp_norm = 0.0;
+    vector_t delta_x_arr[n];
+
+    get_delta_x(m, n, x_arr, get_non_lin_sys, get_jacobian, delta_x_arr);
+    damp_norm = vector_get_norm2(n, delta_x_arr);
+
+    return damp_norm;
+}
+
+//Computes the correction vector.
+void get_delta_x(uint8_t m, uint8_t n, vector_t x_arr[],
+                 void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+                 void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]),
+                 vector_t delta_x_arr[])
+{
+    matrix_t J[m][n];
+    matrix_t pinv_J[n][m];
+    vector_t f_vec[m];
+
+    // compute the Jacobian matrix for the vector x_arr.
+    get_jacobian(x_arr, J);
+    // compute the inverse of J (J^-1)
+    moore_penrose_get_pinv(m, n, J, pinv_J);
+    // compute the value of f(x)
+    get_non_lin_sys(x_arr, f_vec);
+    //delta_x = J1^-1*f(x)
+    matrix_mul_vec(n, m, pinv_J, f_vec, delta_x_arr);
+    vector_in_place_scalar_mul(n, delta_x_arr, -1);
+}

non_linear_algebra/solve_non_linear_equations/doc.txt

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de> 
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup    solve_non_linear_equations SOLVE_NON_LINEAR_EQUATIONS
+ * @ingroup     non_linear_algebra
+ *
+ * @brief       Enables to solve multi-variant nonlinear equation systems.
+ *
+ * @details     The multi-variant nonlinear equation systems are solved using 
+ *              damped or the Newton--Raphson algorithms.
+ *
+ * @author      Zakaria Kasmi
+ */

non_linear_algebra/solve_non_linear_equations/fsolve.c

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief       Solve multi-variant nonlinear equation systems.
+ * @details     The multi-variant nonlinear equation systems are solved using
+ *              damped or the Newton--Raphson algorithms.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include "vector.h"
+#include "matrix.h"
+#include "fsolve.h"
+#include "newton_raphson.h"
+#include "damped_newton_raphson.h"
+
+uint8_t fsolve(uint8_t f_length, uint8_t x0_length, vector_t x0_arr[],
+               enum NON_LIN_ALGORITHM algo, vector_t est_x_arr[],
+               void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+               void (*get_jacobian)(vector_t x_arr[], matrix_t J[][x0_length]))
+{
+    double tol = 1e-9;
+    uint8_t max_it_num = 77;
+    uint8_t iter_num = 0;
+
+    switch (algo) {
+        case Newton_Raphson:
+            iter_num = newton_raphson(f_length, x0_length, x0_arr, tol,
+                                      max_it_num,
+                                      est_x_arr, get_non_lin_sys, get_jacobian);
+            break;
+
+        case Damped_Newton_Raphson:
+        {
+            double min_lamda = 4.8828125e-04;
+            iter_num = damped_newton_raphson(f_length, x0_length, x0_arr,
+                                             min_lamda,
+                                             tol, max_it_num, est_x_arr,
+                                             get_non_lin_sys, get_jacobian);
+            break;
+        }
+
+        default:
+            iter_num = newton_raphson(f_length, x0_length, x0_arr, tol,
+                                      max_it_num,
+                                      est_x_arr, get_non_lin_sys, get_jacobian);
+    }
+
+    return iter_num;
+}

non_linear_algebra/solve_non_linear_equations/include/damped_newton_raphson.h

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup     solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief       Implement the damped Newton--Raphson algorithm.
+ * @details     The damped Newton--Raphson algorithm enables to solve
+ *              multi-variant nonlinear equation systems.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef DAMPED_NEWTON_RAPHSON_H_
+#define DAMPED_NEWTON_RAPHSON_H_
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief   Implements the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] n                        length of the start vector.
+ * @param[in] x0_arr[]                 start vector.
+ * @param[in] min_lamda                minimal damping factor.
+ * @param[in] eps                      accuracy bound.
+ * @param[in] max_it_num               maximal iteration number of the damped Newton--Raphson
+ *                                     algorithm.
+ * @param[out] est_x_arr[]             estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys)       pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  required iteration number.
+ *
+ */
+uint8_t damped_newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+                              double min_lamda, double eps, uint8_t max_it_num,
+                              vector_t est_x_arr[],
+                              void (*get_non_lin_sys)(vector_t x_arr[],
+                                                      vector_t f_vec[]),
+                              void (*get_jacobian)(vector_t x_arr[],
+                                                   matrix_t J[][n]));
+/**
+ * @brief   Compute the norm of the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] m                        number of the non-linear equations.
+ * @param[in] n                        length of the guess vector.
+ * @param[in] x_arr[]                  guess vector.
+ * @param[in] (*get_non_lin_sys)       pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  norm of the damped Newton--Raphson algorithm.
+ *
+ */
+double get_damped_norm(uint8_t m, uint8_t n, vector_t x_arr[],
+                       void (*get_non_lin_sys)(vector_t x_arr[],
+                                               vector_t f_vec[]),
+                       void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n])
+                       );
+
+/**
+ * @brief   Compute the correction vector the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] m                        number of the non-linear equations.
+ * @param[in] n                        length of the guess vector.
+ * @param[in] x_arr[]                  guess vector.
+ * @param[in] (*get_non_lin_sys)       pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ * @param[in, out] delta_x_arr[]       the correction vector (term).
+ *
+ */
+void get_delta_x(uint8_t m, uint8_t n, vector_t x_arr[],
+                 void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+                 void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]),
+                 vector_t delta_x_arr[]);
+
+#endif /* DAMPED_NEWTON_RAPHSON_H_ */

non_linear_algebra/solve_non_linear_equations/include/fsolve.h

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup     solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief       Solve multi-variant nonlinear equation systems.
+ * @details     The multi-variant nonlinear equation systems are solved using
+ *              damped or the Newton--Raphson algorithms.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef FSOLVE_H_
+#define FSOLVE_H_
+
+/**
+ * Possible algorithms to solve multi-variant nonlinear equation systems.
+ */
+enum NON_LIN_ALGORITHM {
+    Newton_Raphson,         /**< Newton--Raphson algorithm */
+    Damped_Newton_Raphson   /**< Damped Newton--Raphson algorithm */
+};
+
+/**
+ * @brief   Solve systems of multi-variant nonlinear equations.
+ *  @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ *           The user can choose between the damped or the Newton--Raphson algorithms.
+ *
+ * @note    This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] x0_length                length of the start vector.
+ * @param[in] x0_arr[]                 start vector.
+ * @param[in] algo                     damped or the Newton--Raphson algorithm.
+ * @param[in, out] est_x_arr[]         estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys)       pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  required iteration number.
+ *
+ */
+uint8_t fsolve(uint8_t f_length, uint8_t x0_length, vector_t x0_arr[],
+               enum NON_LIN_ALGORITHM algo, vector_t est_x_arr[],
+               void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+               void (*get_jacobian)(vector_t x_arr[], matrix_t J[][x0_length]));
+
+#endif /* FSOLVE_H_ */

non_linear_algebra/solve_non_linear_equations/include/newton_raphson.h

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup     solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief       Implement the Newton--Raphson algorithm.
+ * @details     The Newton--Raphson algorithm enables to solve
+ *              multi-variant nonlinear equation systems.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+#ifndef NEWTON_RAPHSON_H_
+#define NEWTON_RAPHSON_H_
+
+#include <inttypes.h>
+
+#include "vector.h"
+#include "matrix.h"
+
+/**
+ * @brief   Implements the Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note    This function is generally implemented.
+ *
+ * @param[in] f_length                 length of the error functions vector.
+ * @param[in] n                        length of the start vector.
+ * @param[in] x0_arr[]                 start vector.
+ * @param[in] eps                      accuracy bound.
+ * @param[in] max_it_num               maximal iteration number of the Newton--Raphson algorithm.
+ * @param[out] est_x_arr[]             estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys)       pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian)          pointer to the Jacobian matrix.
+ *
+ * @return  required iteration number.
+ *
+ */
+uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+                       double eps, uint8_t max_it_num, vector_t est_x_arr[],
+                       void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+                       void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]));
+
+#endif /* NEWTON_RAPHSON_H_ */

non_linear_algebra/solve_non_linear_equations/newton_raphson.c

+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *               2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief       Implement the Newton--Raphson algorithm.
+ * @details     The Newton--Raphson algorithm enables to solve
+ *              multi-variant nonlinear equation systems.
+ *
+ * @author      Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include "vector.h"
+#include "matrix.h"
+#include "moore_penrose_pseudo_inverse.h"
+
+//n is the size of x0 that is equal to the column number of J
+uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+                       double eps, uint8_t max_it_num, vector_t est_x_arr[],
+                       void (*get_non_lin_sys)(vector_t x_arr[],
+                                               vector_t f_vec[]),
+                       void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]))
+{
+    uint8_t iter_num;
+    double step;
+    vector_t prev_x_arr[n];
+    matrix_t J[f_length][n];
+    matrix_t pinv_J[n][f_length];
+    vector_t delta_x[n];
+    vector_t f_vec[f_length];
+
+    vector_copy(n, x0_arr, prev_x_arr);
+    step = eps;
+    iter_num = 0;
+
+    while ((step >= eps) && (iter_num < max_it_num)) {
+        // compute the Jacobian matrix for the prev_x_arr.
+        get_jacobian(prev_x_arr, J);
+        // compute the inverse of J (J^-1)
+        moore_penrose_get_pinv(f_length, n, J, pinv_J);
+        // compute the value of f(x)
+        get_non_lin_sys(prev_x_arr, f_vec);
+        //delta_x = J1^-1*f(x)
+        matrix_mul_vec(n, f_length, pinv_J, f_vec, delta_x);
+        //improve x: x = x - delta_x
+        vector_sub(n, prev_x_arr, delta_x, est_x_arr);
+        //next step
+        step = vector_get_euclidean_distance(n, est_x_arr, prev_x_arr);
+        //update prev_x
+        vector_copy(n, est_x_arr, prev_x_arr);
+        iter_num++;
+    }
+
+    return iter_num;
+}