Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
qr_algorithm Namespace Reference

Functions to compute QR decomposition of any rectangular matrix. More...

Functions

template<typename T >
std::ostreamoperator<< (std::ostream &out, std::valarray< std::valarray< T >> const &v)
 
template<typename T >
std::ostreamoperator<< (std::ostream &out, std::valarray< T > const &v)
 
template<typename T >
double vector_dot (const std::valarray< T > &a, const std::valarray< T > &b)
 
template<typename T >
double vector_mag (const std::valarray< T > &a)
 
template<typename T >
std::valarray< T > vector_proj (const std::valarray< T > &a, const std::valarray< T > &b)
 
template<typename T >
void qr_decompose (const std::valarray< std::valarray< T >> &A, std::valarray< std::valarray< T >> *Q, std::valarray< std::valarray< T >> *R)
 
std::valarray< double > eigen_values (std::valarray< std::valarray< double >> *A, bool print_intermediates=false)
 

Detailed Description

Functions to compute QR decomposition of any rectangular matrix.

Function Documentation

◆ eigen_values()

std::valarray<double> qr_algorithm::eigen_values ( std::valarray< std::valarray< double >> *  A,
bool  print_intermediates = false 
)

Compute eigen values using iterative shifted QR decomposition algorithm as follows:

  1. Use last diagonal element of A as eigen value approximation \(c\)
  2. Shift diagonals of matrix \(A' = A - cI\)
  3. Decompose matrix \(A'=QR\)
  4. Compute next approximation \(A'_1 = RQ \)
  5. Shift diagonals back \(A_1 = A'_1 + cI\)
  6. Termination condition check: last element below diagonal is almost 0
    1. If not 0, go back to step 1 with the new approximation \(A_1\)
    2. If 0, continue to step 7
  7. Save last known \(c\) as the eigen value.
  8. Are all eigen values found?
    1. If not, remove last row and column of \(A_1\) and go back to step 1.
    2. If yes, stop.
Note
The matrix \(A\) gets modified
Parameters
[in,out]Amatrix to compute eigen values for
[in]print_intermediates(optional) whether to print intermediate A, Q and R matrices (default = false)
99  {
100  int rows = A->size();
101  int columns = rows;
102  int counter = 0, num_eigs = rows - 1;
103  double last_eig = 0;
104 
107 
108  /* number of eigen values = matrix size */
109  std::valarray<double> eigen_vals(rows);
110  for (int i = 0; i < rows; i++) {
111  Q[i] = std::valarray<double>(columns);
112  R[i] = std::valarray<double>(columns);
113  }
114 
115  /* continue till all eigen values are found */
116  while (num_eigs > 0) {
117  /* iterate with QR decomposition */
118  while (std::abs(A[0][num_eigs][num_eigs - 1]) >
120  // initial approximation = last diagonal element
121  last_eig = A[0][num_eigs][num_eigs];
122  for (int i = 0; i < rows; i++) {
123  A[0][i][i] -= last_eig; /* A - cI */
124  }
125 
126  qr_decompose(*A, &Q, &R);
127 
128  if (print_intermediates) {
129  std::cout << *A << "\n";
130  std::cout << Q << "\n";
131  std::cout << R << "\n";
132  printf("-------------------- %d ---------------------\n",
133  ++counter);
134  }
135 
136  // new approximation A' = R * Q
137  mat_mul(R, Q, A);
138 
139  for (int i = 0; i < rows; i++) {
140  A[0][i][i] += last_eig; /* A + cI */
141  }
142  }
143 
144  /* store the converged eigen value */
145  eigen_vals[num_eigs] = last_eig;
146  // A[0][num_eigs][num_eigs];
147  if (print_intermediates) {
148  std::cout << "========================\n";
149  std::cout << "Eigen value: " << last_eig << ",\n";
150  std::cout << "========================\n";
151  }
152 
153  num_eigs--;
154  rows--;
155  columns--;
156  }
157  eigen_vals[0] = A[0][0][0];
158 
159  if (print_intermediates) {
160  std::cout << Q << "\n";
161  std::cout << R << "\n";
162  }
163 
164  return eigen_vals;
165 }
Here is the call graph for this function:

◆ operator<<() [1/2]

template<typename T >
std::ostream& qr_algorithm::operator<< ( std::ostream out,
std::valarray< std::valarray< T >> const &  v 
)

operator to print a matrix

34  {
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 }
Here is the call graph for this function:

◆ operator<<() [2/2]

template<typename T >
std::ostream& qr_algorithm::operator<< ( std::ostream out,
std::valarray< T > const &  v 
)

operator to print a vector

53  {
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 }
Here is the call graph for this function:

◆ qr_decompose()

template<typename T >
void qr_algorithm::qr_decompose ( const std::valarray< std::valarray< T >> &  A,
std::valarray< std::valarray< T >> *  Q,
std::valarray< std::valarray< T >> *  R 
)

Decompose matrix \(A\) using Gram-Schmidt process.

\begin{eqnarray*} \text{given that}\quad A &=& *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\ \text{where}\quad\mathbf{a}_i &=& \left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column vectors)}\\ \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\ \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\ Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots & \mathbf{e}_{N-1}\end{bmatrix}\\ R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\ 0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ 0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}\\ \end{eqnarray*}

Parameters
Ainput matrix to decompose
Qoutput decomposed matrix
Routput 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 }
Here is the call graph for this function:

◆ vector_dot()

template<typename T >
double qr_algorithm::vector_dot ( const std::valarray< T > &  a,
const std::valarray< T > &  b 
)
inline

Compute dot product of two vectors of equal lengths

If \(\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\) and \(\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\) then \(\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\)

Returns
\(\vec{a}\cdot\vec{b}\)
76  {
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 }

◆ vector_mag()

template<typename T >
double qr_algorithm::vector_mag ( const std::valarray< T > &  a)
inline

Compute magnitude of vector.

If \(\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\) then \(\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\)

Returns
\(\left|\vec{a}\right|\)
92  {
93  double dot = vector_dot(a, a);
94  return std::sqrt(dot);
95 }
Here is the call graph for this function:

◆ vector_proj()

template<typename T >
std::valarray<T> qr_algorithm::vector_proj ( const std::valarray< T > &  a,
const std::valarray< T > &  b 
)

Compute projection of vector \(\vec{a}\) on \(\vec{b}\) defined as

\[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\]

Returns
NULL if error, otherwise pointer to output

check for division by zero using machine epsilon

105  {
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 }
Here is the call graph for this function:
std::setfill
T setfill(T... args)
mat_mul
void mat_mul(const std::valarray< std::valarray< double >> &A, const std::valarray< std::valarray< double >> &B, std::valarray< std::valarray< double >> *OUT)
Definition: qr_eigen_values.cpp:54
std::sqrt
T sqrt(T... args)
std::printf
T printf(T... args)
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
qr_algorithm::vector_mag
double vector_mag(const std::valarray< T > &a)
Definition: qr_decompose.h:92
std::valarray
STL class.
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)
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
std::size_t
std::setw
T setw(T... args)
std::numeric_limits
std::ostream::precision
T precision(T... args)
machine_learning::sum
T sum(const std::vector< std::valarray< T >> &A)
Definition: vector_ops.hpp:232