scijs/ndarray-givens-qr

Name: ndarray-givens-qr

Owner: scijs

Description: QR decomposition using Givens rotations.

Created: 2016-04-11 04:20:19.0

Updated: 2016-12-13 06:57:05.0

Pushed: 2016-05-03 23:16:06.0

Homepage: null

Size: 8

Language: JavaScript

GitHub Committers

User	Most Recent Commit	# Commits

Other Committers

User	Email	Most Recent Commit	# Commits

README

ndarray-givens-qr

QR decomposition using Givens rotations.

Usage

Solving Ax=b Problems

qrDecomp = require('ndarray-givens-qr');
ndarray = require('ndarray');

A = ndarray(new Float64Array([...]), [m, n]);
b = ndarray(new Float64Array([...]), [m]);
x = ndarray(new Float64Array([...]), [n]);

olve(A, b, x);

Alternatively, if no output vector is specified, then the vector is created and returned as shown:

qrDecomp = require('ndarray-givens-qr');
ndarray = require('ndarray');

A = ndarray(new Float64Array([...]), [m, n]);
b = ndarray(new Float64Array([...]), [m]);

x = qr.solve(A, b);

Matrix Factorization

qrDecomp = require('ndarray-givens-qr');
ndarray = require('ndarray');

A = ndarray(new Float64Array([...]), [m, n]);
Q = ndarray(new Float64Array([...]), [m, m]);
R = ndarray(new Float64Array([...]), [m, n]);

ecompose(A, Q, R);

Alternatively, if no output matrices are specified, then proper matrices are created and returned as shown:

qr = require('ndarray-givens-qr');
ndarray = require('ndarray');

A = ndarray(new Float64Array([...]), [m, n]);

results = qr.decompose(A);
Q = results.q;
R = results.r;

Background

For a matrix A with m rows and n columns, QR decompositions create an m x m matrix Q and an m x n matrix R, where Q is a unitary matrix and R is upper triangular.

The idea behind using Givens rotations is clearing out the zeros beneath the diagonal entries of A. A rotation matrix that rotates a vector on the X-axis can be rearranged to the following form:

   -+ +- -+   +- -+
  s | | a | = | r |
  c | | b |   | 0 |
   -+ +- -+   +- -+

This rearranging effectively “clears out” the last entry in the column vector. Givens rotations generalize this so that any row i and j can be “cleared out”.

Algorithm

The algorithm computes the Givens rotation using BLAS Level 1. The Givens rotation matrix is:

         +-                         -+
j,c,s) = | 1  ...  0  ...  0  ...  0 |
         |                           |
         | |   \   |       |       | |
         |                           |
         | 0  ...  c  ...  s  ...  0 |
         |                           |
         | |       |   \   |       | |
         |                           |
         | 0  ... -s  ...  c  ...  0 |
         |                           |
         | |       |       |   \   | |
         |                           |
         | 0  ...  0  ...  0  ...  1 |
         +-                         -+

where c and s are the values computed from the ROTG from BLAS Level 1. These values only appear at the intersection of the ith and jth rows and columns.

In pseudocode, the algorithm is:

A
I
j = 0 : 1 : n-1
r i = m-1 : -1 : j+1
a = A(i-1, j)
b = A(i, j)
[c, s, r] = BLAS1.rotg(a, b)
G = givens(i, j, c, s)
R = G * R
Q = Q * G^T
d

Future Work

This algorithm can be parallelized. Since each Givens rotation only affects the ith and jth rows of the R matrix, more than one column can be updated at a time. The stages at which a subdiagonal entry can be annihilated (“cleared out”) for an 8x8 matrix is given as:

                           -+
                            |
  *                         |
  8   *                     |
  7   9   *                 |
  6   8   10  *             |
  5   7   9   11  *         |
  4   6   8   10  12  *     |
  3   5   7   9   11  13  * |
                           -+

License

© 2015 Scijs. MIT License.

Author

Tim Bright