RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
damped_newton_raphson.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
23 #include <stdio.h>
24 
26 #include "damped_newton_raphson.h"
27 #include "vector.h"
28 #include "matrix.h"
29 
30 //Multivariate Newton’s Method
31 //n is the size of x0 that is equal to the column number of J
32 uint8_t damped_newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
33  double min_lamda, double eps, uint8_t max_it_num,
34  vector_t est_x_arr[],
35  void (*get_non_lin_sys)(vector_t x_arr[],
36  vector_t f_vec[]),
37  void (*get_jacobian)(vector_t x_arr[],
38  matrix_t J[][n]))
39 {
40  uint8_t iter_num;
41  double lamda;
42  double lamda_c;
43  double damped_norm_x;
44  double damped_norm_prev_x;
45  vector_t prev_x_arr[n];
46  vector_t delta_x[n];
47  vector_t lamda_mul_delta_x[n];
48 
49  damped_norm_x = get_damped_norm(f_length, n, x0_arr, get_non_lin_sys,
50  get_jacobian);
51  iter_num = 0;
52  vector_copy(n, x0_arr, prev_x_arr);
53  while ((damped_norm_x >= eps) && (iter_num < max_it_num)) {
54 
55  get_delta_x(f_length, n, prev_x_arr, get_non_lin_sys,
56  get_jacobian, delta_x);
57  lamda = 1.0;
58 
59  // x = x_k +lamda*s_k (lamda=1)
60  vector_add(n, prev_x_arr, delta_x, est_x_arr);
61 
62  lamda_c = 1 - lamda / 4;
63  //|||f(x)|||_k
64  damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
65  get_non_lin_sys, get_jacobian);
66 
67  //|||f(prev_x)|||_k
68  damped_norm_prev_x = vector_get_norm2(n, delta_x);
69 
70  while (damped_norm_x > (lamda_c * damped_norm_prev_x)) {
71  lamda = lamda / 2;
72  if (lamda >= min_lamda) {
73  //x = x_k +lamda*s_k
74  vector_scalar_mul(n, delta_x, lamda,
75  lamda_mul_delta_x);
76  vector_add(n, prev_x_arr, lamda_mul_delta_x,
77  est_x_arr);
78  }
79  else {
80  break;
81  }
82  // |||f(x)|||_k
83  damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
84  get_non_lin_sys, get_jacobian);
85 
86  } //while
87 
88  if (lamda <= min_lamda) {
89  break;
90  }
91 
92  //update prev_x
93  vector_copy(n, est_x_arr, prev_x_arr);
94 
95  // |||f(x)|||_k
96  damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
97  get_non_lin_sys, get_jacobian);
98  iter_num++;
99  }
100 
101  return iter_num;
102 }
103 
104 double get_damped_norm(uint8_t m, uint8_t n, vector_t x_arr[],
105  void (*get_non_lin_sys)(vector_t x_arr[],
106  vector_t f_vec[]),
107  void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n])
108  )
109 {
110  double damp_norm = 0.0;
111  vector_t delta_x_arr[n];
112 
113  get_delta_x(m, n, x_arr, get_non_lin_sys, get_jacobian, delta_x_arr);
114  damp_norm = vector_get_norm2(n, delta_x_arr);
115 
116  return damp_norm;
117 }
118 
119 //Computes the correction vector.
120 void get_delta_x(uint8_t m, uint8_t n, vector_t x_arr[],
121  void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
122  void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]),
123  vector_t delta_x_arr[])
124 {
125  matrix_t J[m][n];
126  matrix_t pinv_J[n][m];
127  vector_t f_vec[m];
128 
129  // compute the Jacobian matrix for the vector x_arr.
130  get_jacobian(x_arr, J);
131  // compute the inverse of J (J^-1)
132  moore_penrose_get_pinv(m, n, J, pinv_J);
133  // compute the value of f(x)
134  get_non_lin_sys(x_arr, f_vec);
135  //delta_x = J1^-1*f(x)
136  matrix_mul_vec(n, m, pinv_J, f_vec, delta_x_arr);
137  vector_in_place_scalar_mul(n, delta_x_arr, -1);
138 }
moore_penrose_get_pinv
int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t pinv_A[n][m])
Calculate the Moore–Penrose inverse of a rectangular matrix.
Definition: moore_penrose_pseudo_inverse.c:38
vector_get_norm2
vector_t vector_get_norm2(uint8_t length, vector_t arr[])
Compute the 2-norm norm of a vector.
Definition: vector.c:42
matrix_mul_vec
void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t vec[n], matrix_t dst_arr[m])
Compute the multiplication of a matrix with a column vector.
Definition: matrix.c:434
get_delta_x
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[])
Compute the correction vector the damped Newton–Raphson algorithm.
Definition: damped_newton_raphson.c:120
vector_t
#define vector_t
Define the data type of the vector elements.
Definition: vector.h:33
damped_newton_raphson
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]))
Implements the damped Newton–Raphson algorithm.
Definition: damped_newton_raphson.c:32
vector_in_place_scalar_mul
void vector_in_place_scalar_mul(uint8_t size, vector_t a_vec[size], vector_t scl)
Compute the product of a vector with a real number.
Definition: vector.c:131
matrix.h
Matrix computations.
damped_newton_raphson.h
Implement the damped Newton–Raphson algorithm.
moore_penrose_pseudo_inverse.h
Moore–Penrose algorithm to compute the pseudo-inverse of a matrix.
get_damped_norm
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]))
Compute the norm of the damped Newton–Raphson algorithm.
Definition: damped_newton_raphson.c:104
vector_copy
void vector_copy(uint8_t size, vector_t src_arr[], vector_t dest_arr[])
Copy the elements of the source vector to the destination vector.
Definition: vector.c:37
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
vector_add
void vector_add(uint8_t size, vector_t a_vec[size], vector_t b_vec[size], vector_t a_plus_b_vec[size])
Compute the addition of two vectors.
Definition: vector.c:104
vector.h
Vector computations.
vector_scalar_mul
void vector_scalar_mul(uint8_t size, vector_t src_vec[size], vector_t scl, vector_t dest_vec[])
Compute the product of a vector with a real number.
Definition: vector.c:141