From acc50fcab73fd613566ea8fc89a54ff2e6879e08 Mon Sep 17 00:00:00 2001 From: Krishna Vedala <7001608+kvedala@users.noreply.github.com> Date: Mon, 8 Jun 2020 19:35:14 -0400 Subject: [PATCH] added QR decomposition algorithm Signed-off-by: Krishna Vedala <7001608+kvedala@users.noreply.github.com> --- numerical_methods/qr_decompose.h | 210 +++++++++++++++++++++++++ numerical_methods/qr_decomposition.cpp | 58 +++++++ 2 files changed, 268 insertions(+) create mode 100644 numerical_methods/qr_decompose.h create mode 100644 numerical_methods/qr_decomposition.cpp diff --git a/numerical_methods/qr_decompose.h b/numerical_methods/qr_decompose.h new file mode 100644 index 000000000..c9c369144 --- /dev/null +++ b/numerical_methods/qr_decompose.h @@ -0,0 +1,210 @@ +/** + * @file + * \brief Library functions to compute [QR + * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given + * matrix. + * \author [Krishna Vedala](https://github.com/kvedala) + */ + +#ifndef NUMERICAL_METHODS_QR_DECOMPOSE_H_ +#define NUMERICAL_METHODS_QR_DECOMPOSE_H_ + +#include +#include +#include +#include +#include +#include +#include +#ifdef _OPENMP +#include +#endif + +/** \namespace qr_algorithm + * \brief Functions to compute [QR + * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of any + * rectangular matrix + */ +namespace qr_algorithm { +/** + * operator to print a matrix + */ +template +std::ostream &operator<<(std::ostream &out, + std::valarray> const &v) { + const int width = 12; + const char separator = ' '; + + out.precision(4); + for (size_t row = 0; row < v.size(); row++) { + for (size_t col = 0; col < v[row].size(); col++) + out << std::right << std::setw(width) << std::setfill(separator) + << v[row][col]; + out << std::endl; + } + + return out; +} + +/** + * operator to print a vector + */ +template +std::ostream &operator<<(std::ostream &out, std::valarray const &v) { + const int width = 10; + const char separator = ' '; + + out.precision(4); + for (size_t row = 0; row < v.size(); row++) { + out << std::right << std::setw(width) << std::setfill(separator) + << v[row]; + } + + return out; +} + +/** + * Compute dot product of two vectors of equal lengths + * + * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ and + * \f$\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\f$ then + * \f$\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\f$ + * + * \returns \f$\vec{a}\cdot\vec{b}\f$ + */ +template +inline double vector_dot(const std::valarray &a, const std::valarray &b) { + return (a * b).sum(); + // could also use following + // return std::inner_product(std::begin(a), std::end(a), std::begin(b), + // 0.f); +} + +/** + * Compute magnitude of vector. + * + * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ then + * \f$\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\f$ + * + * \returns \f$\left|\vec{a}\right|\f$ + */ +template +inline double vector_mag(const std::valarray &a) { + double dot = vector_dot(a, a); + return std::sqrt(dot); +} + +/** + * Compute projection of vector \f$\vec{a}\f$ on \f$\vec{b}\f$ defined as + * \f[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\f] + * + * \returns NULL if error, otherwise pointer to output + */ +template +std::valarray vector_proj(const std::valarray &a, + const std::valarray &b) { + double num = vector_dot(a, b); + double deno = vector_dot(b, b); + + /*! check for division by zero using machine epsilon */ + if (deno <= std::numeric_limits::epsilon()) { + std::cerr << "[" << __func__ << "] Possible division by zero\n"; + return a; // return vector a back + } + + double scalar = num / deno; + + return b * scalar; +} + +/** + * Decompose matrix \f$A\f$ using [Gram-Schmidt + *process](https://en.wikipedia.org/wiki/QR_decomposition). + * + * \f{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}\\ + * \f} + */ +template +void qr_decompose( + const std::valarray> &A, /**< input matrix to decompose */ + std::valarray> *Q, /**< output decomposed matrix */ + std::valarray> *R /**< output decomposed matrix */ +) { + std::size_t ROWS = A.size(); // number of rows of A + std::size_t COLUMNS = A[0].size(); // number of columns of A + std::valarray col_vector(ROWS); + std::valarray col_vector2(ROWS); + std::valarray tmp_vector(ROWS); + + for (int i = 0; i < COLUMNS; i++) { + /* for each column => R is a square matrix of NxN */ + int j; + R[0][i] = 0.; /* make R upper triangular */ + + /* get corresponding Q vector */ +#ifdef _OPENMP +// parallelize on threads +#pragma omp for +#endif + for (j = 0; j < ROWS; j++) { + tmp_vector[j] = A[j][i]; /* accumulator for uk */ + col_vector[j] = A[j][i]; + } + for (j = 0; j < i; j++) { + for (int k = 0; k < ROWS; k++) { + col_vector2[k] = Q[0][k][j]; + } + col_vector2 = vector_proj(col_vector, col_vector2); + tmp_vector -= col_vector2; + } + + double mag = vector_mag(tmp_vector); + +#ifdef _OPENMP +// parallelize on threads +#pragma omp for +#endif + for (j = 0; j < ROWS; j++) Q[0][j][i] = tmp_vector[j] / mag; + + /* compute upper triangular values of R */ +#ifdef _OPENMP +// parallelize on threads +#pragma omp for +#endif + for (int kk = 0; kk < ROWS; kk++) { + col_vector[kk] = Q[0][kk][i]; + } + +#ifdef _OPENMP +// parallelize on threads +#pragma omp for +#endif + for (int k = i; k < COLUMNS; k++) { + for (int kk = 0; kk < ROWS; kk++) { + col_vector2[kk] = A[kk][k]; + } + R[0][i][k] = (col_vector * col_vector2).sum(); + } + } +} +} // namespace qr_algorithm + +#endif // NUMERICAL_METHODS_QR_DECOMPOSE_H_ diff --git a/numerical_methods/qr_decomposition.cpp b/numerical_methods/qr_decomposition.cpp new file mode 100644 index 000000000..237a5c946 --- /dev/null +++ b/numerical_methods/qr_decomposition.cpp @@ -0,0 +1,58 @@ +/** + * @file + * \brief Program to compute the [QR + * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given + * matrix. + * \author [Krishna Vedala](https://github.com/kvedala) + */ + +#include +#include +#include +#include +#include + +#include "./qr_decompose.h" + +using qr_algorithm::qr_decompose; +using qr_algorithm::operator<<; + +/** + * main function + */ +int main(void) { + unsigned int ROWS, COLUMNS; + + std::cout << "Enter the number of rows and columns: "; + std::cin >> ROWS >> COLUMNS; + + std::cout << "Enter matrix elements row-wise:\n"; + + std::valarray> A(ROWS); + std::valarray> Q(ROWS); + std::valarray> R(COLUMNS); + for (int i = 0; i < std::max(ROWS, COLUMNS); i++) { + if (i < ROWS) { + A[i] = std::valarray(COLUMNS); + Q[i] = std::valarray(COLUMNS); + } + if (i < COLUMNS) { + R[i] = std::valarray(COLUMNS); + } + } + + for (int i = 0; i < ROWS; i++) + for (int j = 0; j < COLUMNS; j++) std::cin >> A[i][j]; + + std::cout << A << "\n"; + + clock_t t1 = clock(); + qr_decompose(A, &Q, &R); + double dtime = static_cast(clock() - t1) / CLOCKS_PER_SEC; + + std::cout << Q << "\n"; + std::cout << R << "\n"; + std::cout << "Time taken to compute: " << dtime << " sec\n "; + + return 0; +}