Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
qr_decompose.h
Go to the documentation of this file.
1 /**
2  * @file
3  * \brief Library functions to compute [QR
4  * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
5  * matrix.
6  * \author [Krishna Vedala](https://github.com/kvedala)
7  */
8 
9 #ifndef NUMERICAL_METHODS_QR_DECOMPOSE_H_
10 #define NUMERICAL_METHODS_QR_DECOMPOSE_H_
11 
12 #include <cmath>
13 #include <cstdlib>
14 #include <iomanip>
15 #include <iostream>
16 #include <limits>
17 #include <numeric>
18 #include <valarray>
19 #ifdef _OPENMP
20 #include <omp.h>
21 #endif
22 
23 /** \namespace qr_algorithm
24  * \brief Functions to compute [QR
25  * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of any
26  * rectangular matrix
27  */
28 namespace qr_algorithm {
29 /**
30  * operator to print a matrix
31  */
32 template <typename T>
34  std::valarray<std::valarray<T>> const &v) {
35  const int width = 12;
36  const char separator = ' ';
37 
38  out.precision(4);
39  for (size_t row = 0; row < v.size(); row++) {
40  for (size_t col = 0; col < v[row].size(); col++)
41  out << std::right << std::setw(width) << std::setfill(separator)
42  << v[row][col];
43  out << std::endl;
44  }
45 
46  return out;
47 }
48 
49 /**
50  * operator to print a vector
51  */
52 template <typename T>
54  const int width = 10;
55  const char separator = ' ';
56 
57  out.precision(4);
58  for (size_t row = 0; row < v.size(); row++) {
59  out << std::right << std::setw(width) << std::setfill(separator)
60  << v[row];
61  }
62 
63  return out;
64 }
65 
66 /**
67  * Compute dot product of two vectors of equal lengths
68  *
69  * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ and
70  * \f$\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\f$ then
71  * \f$\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\f$
72  *
73  * \returns \f$\vec{a}\cdot\vec{b}\f$
74  */
75 template <typename T>
76 inline double vector_dot(const std::valarray<T> &a, const std::valarray<T> &b) {
77  return (a * b).sum();
78  // could also use following
79  // return std::inner_product(std::begin(a), std::end(a), std::begin(b),
80  // 0.f);
81 }
82 
83 /**
84  * Compute magnitude of vector.
85  *
86  * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ then
87  * \f$\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\f$
88  *
89  * \returns \f$\left|\vec{a}\right|\f$
90  */
91 template <typename T>
92 inline double vector_mag(const std::valarray<T> &a) {
93  double dot = vector_dot(a, a);
94  return std::sqrt(dot);
95 }
96 
97 /**
98  * Compute projection of vector \f$\vec{a}\f$ on \f$\vec{b}\f$ defined as
99  * \f[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\f]
100  *
101  * \returns NULL if error, otherwise pointer to output
102  */
103 template <typename T>
105  const std::valarray<T> &b) {
106  double num = vector_dot(a, b);
107  double deno = vector_dot(b, b);
108 
109  /*! check for division by zero using machine epsilon */
110  if (deno <= std::numeric_limits<double>::epsilon()) {
111  std::cerr << "[" << __func__ << "] Possible division by zero\n";
112  return a; // return vector a back
113  }
114 
115  double scalar = num / deno;
116 
117  return b * scalar;
118 }
119 
120 /**
121  * Decompose matrix \f$A\f$ using [Gram-Schmidt
122  *process](https://en.wikipedia.org/wiki/QR_decomposition).
123  *
124  * \f{eqnarray*}{
125  * \text{given that}\quad A &=&
126  *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\
127  * \text{where}\quad\mathbf{a}_i &=&
128  * \left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column
129  * vectors)}\\
130  * \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i
131  *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\
132  * \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\
133  * Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots &
134  * \mathbf{e}_{N-1}\end{bmatrix}\\
135  * R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle &
136  * \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
137  * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\
138  * 0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
139  * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\
140  * 0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle &
141  * \dots\\ \vdots & \vdots & \vdots & \ddots
142  * \end{bmatrix}\\
143  * \f}
144  */
145 template <typename T>
147  const std::valarray<std::valarray<T>> &A, /**< input matrix to decompose */
148  std::valarray<std::valarray<T>> *Q, /**< output decomposed matrix */
149  std::valarray<std::valarray<T>> *R /**< output decomposed matrix */
150 ) {
151  std::size_t ROWS = A.size(); // number of rows of A
152  std::size_t COLUMNS = A[0].size(); // number of columns of A
153  std::valarray<T> col_vector(ROWS);
154  std::valarray<T> col_vector2(ROWS);
155  std::valarray<T> tmp_vector(ROWS);
156 
157  for (int i = 0; i < COLUMNS; i++) {
158  /* for each column => R is a square matrix of NxN */
159  int j;
160  R[0][i] = 0.; /* make R upper triangular */
161 
162  /* get corresponding Q vector */
163 #ifdef _OPENMP
164 // parallelize on threads
165 #pragma omp for
166 #endif
167  for (j = 0; j < ROWS; j++) {
168  tmp_vector[j] = A[j][i]; /* accumulator for uk */
169  col_vector[j] = A[j][i];
170  }
171  for (j = 0; j < i; j++) {
172  for (int k = 0; k < ROWS; k++) {
173  col_vector2[k] = Q[0][k][j];
174  }
175  col_vector2 = vector_proj(col_vector, col_vector2);
176  tmp_vector -= col_vector2;
177  }
178 
179  double mag = vector_mag(tmp_vector);
180 
181 #ifdef _OPENMP
182 // parallelize on threads
183 #pragma omp for
184 #endif
185  for (j = 0; j < ROWS; j++) Q[0][j][i] = tmp_vector[j] / mag;
186 
187  /* compute upper triangular values of R */
188 #ifdef _OPENMP
189 // parallelize on threads
190 #pragma omp for
191 #endif
192  for (int kk = 0; kk < ROWS; kk++) {
193  col_vector[kk] = Q[0][kk][i];
194  }
195 
196 #ifdef _OPENMP
197 // parallelize on threads
198 #pragma omp for
199 #endif
200  for (int k = i; k < COLUMNS; k++) {
201  for (int kk = 0; kk < ROWS; kk++) {
202  col_vector2[kk] = A[kk][k];
203  }
204  R[0][i][k] = (col_vector * col_vector2).sum();
205  }
206  }
207 }
208 } // namespace qr_algorithm
209 
210 #endif // NUMERICAL_METHODS_QR_DECOMPOSE_H_
save_exact_solution
void save_exact_solution(const double &X0, const double &X_MAX, const double &step_size, const std::valarray< double > &Y0)
Definition: ode_midpoint_euler.cpp:156
std::srand
T srand(T... args)
semi_implicit_euler
double semi_implicit_euler(double dx, double x0, double x_max, std::valarray< double > *y, bool save_to_file=false)
Compute approximation using the semi-implicit-Euler method in the given limits.
Definition: ode_semi_implicit_euler.cpp:103
semi_implicit_euler_step
void semi_implicit_euler_step(const double dx, const double &x, std::valarray< double > *y, std::valarray< double > *dy)
Compute next step approximation using the semi-implicit-Euler method.
Definition: ode_semi_implicit_euler.cpp:82
std::clock_t
eq
static double eq(double i)
Definition: newton_raphson_method.cpp:29
forward_euler_step
void forward_euler_step(const double dx, const double x, std::valarray< double > *y, std::valarray< double > *dy)
Compute next step approximation using the forward-Euler method.
Definition: ode_forward_euler.cpp:86
std::cos
T cos(T... args)
eq_der
static double eq_der(double i)
Definition: newton_raphson_method.cpp:39
exact_solution
void exact_solution(const double &x, std::valarray< double > *y)
Exact solution of the problem. Used for solution comparison.
Definition: ode_forward_euler.cpp:67
std::atof
T atof(T... args)
main
int main(int argc, char *argv[])
Definition: ode_forward_euler.cpp:189
std::setfill
T setfill(T... args)
midpoint_euler_step
void midpoint_euler_step(const double dx, const double &x, std::valarray< double > *y, std::valarray< double > *dy)
Compute next step approximation using the midpoint-Euler method.
Definition: ode_midpoint_euler.cpp:85
std::sqrt
T sqrt(T... args)
qr_algorithm
Functions to compute QR decomposition of any rectangular matrix.
std::clock
T clock(T... args)
save_exact_solution
void save_exact_solution(const double &X0, const double &X_MAX, const double &step_size, const std::valarray< double > &Y0)
Definition: ode_forward_euler.cpp:153
EPSILON
constexpr double EPSILON
system accuracy limit
Definition: newton_raphson_method.cpp:20
qr_algorithm::vector_proj
std::valarray< T > vector_proj(const std::valarray< T > &a, const std::valarray< T > &b)
Definition: qr_decompose.h:104
std::cout
std::ofstream
STL class.
qr_algorithm::vector_mag
double vector_mag(const std::valarray< T > &a)
Definition: qr_decompose.h:92
MAX_ITERATIONS
constexpr int16_t MAX_ITERATIONS
Maximum number of iterations.
Definition: newton_raphson_method.cpp:21
std::ofstream::close
T close(T... args)
std::perror
T perror(T... args)
std::valarray< double >
std::ofstream::open
T open(T... args)
qr_algorithm::operator<<
std::ostream & operator<<(std::ostream &out, std::valarray< std::valarray< T >> const &v)
Definition: qr_decompose.h:33
forward_euler
double forward_euler(double dx, double x0, double x_max, std::valarray< double > *y, bool save_to_file=false)
Compute approximation using the forward-Euler method in the given limits.
Definition: ode_forward_euler.cpp:102
std::rand
T rand(T... args)
problem
void problem(const double &x, std::valarray< double > *y, std::valarray< double > *dy)
Problem statement for a system with first-order differential equations. Updates the system differenti...
Definition: ode_semi_implicit_euler.cpp:53
std::sin
T sin(T... args)
problem
void problem(const double &x, std::valarray< double > *y, std::valarray< double > *dy)
Problem statement for a system with first-order differential equations. Updates the system differenti...
Definition: ode_forward_euler.cpp:54
std::endl
T endl(T... args)
qr_algorithm::vector_dot
double vector_dot(const std::valarray< T > &a, const std::valarray< T > &b)
Definition: qr_decompose.h:76
std::right
T right(T... args)
problem
void problem(const double &x, std::valarray< double > *y, std::valarray< double > *dy)
Problem statement for a system with first-order differential equations. Updates the system differenti...
Definition: ode_midpoint_euler.cpp:53
main
int main()
Definition: newton_raphson_method.cpp:44
midpoint_euler
double midpoint_euler(double dx, double x0, double x_max, std::valarray< double > *y, bool save_to_file=false)
Compute approximation using the midpoint-Euler method in the given limits.
Definition: ode_midpoint_euler.cpp:107
qr_algorithm::qr_decompose
void qr_decompose(const std::valarray< std::valarray< T >> &A, std::valarray< std::valarray< T >> *Q, std::valarray< std::valarray< T >> *R)
Definition: qr_decompose.h:146
main
int main(int argc, char *argv[])
Definition: ode_semi_implicit_euler.cpp:189
std::size_t
std::time
T time(T... args)
save_exact_solution
void save_exact_solution(const double &X0, const double &X_MAX, const double &step_size, const std::valarray< double > &Y0)
Definition: ode_semi_implicit_euler.cpp:153
exact_solution
void exact_solution(const double &x, std::valarray< double > *y)
Exact solution of the problem. Used for solution comparison.
Definition: ode_semi_implicit_euler.cpp:66
std::setw
T setw(T... args)
std::cin
std::numeric_limits
std::ofstream::is_open
T is_open(T... args)
exact_solution
void exact_solution(const double &x, std::valarray< double > *y)
Exact solution of the problem. Used for solution comparison.
Definition: ode_midpoint_euler.cpp:66
main
int main(int argc, char *argv[])
Definition: ode_midpoint_euler.cpp:192
std::ostream::precision
T precision(T... args)
machine_learning::sum
T sum(const std::vector< std::valarray< T >> &A)
Definition: vector_ops.hpp:232
std::pow
T pow(T... args)