ndarray-givens-qr
QR decomposition using Givens rotations.
Usage
Solving Ax=b Problems
var qrDecomp = require('ndarray-givens-qr');
var ndarray = require('ndarray');
var A = ndarray(new Float64Array([...]), [m, n]);
var b = ndarray(new Float64Array([...]), [m]);
var x = ndarray(new Float64Array([...]), [n]);
qr.solve(A, b, x);
Alternatively, if no output vector is specified, then the vector is created and returned as shown:
var qrDecomp = require('ndarray-givens-qr');
var ndarray = require('ndarray');
var A = ndarray(new Float64Array([...]), [m, n]);
var b = ndarray(new Float64Array([...]), [m]);
var x = qr.solve(A, b);
Matrix Factorization
var qrDecomp = require('ndarray-givens-qr');
var ndarray = require('ndarray');
var A = ndarray(new Float64Array([...]), [m, n]);
var Q = ndarray(new Float64Array([...]), [m, m]);
var R = ndarray(new Float64Array([...]), [m, n]);
qr.decompose(A, Q, R);
Alternatively, if no output matrices are specified, then proper matrices are created and returned as shown:
var qr = require('ndarray-givens-qr');
var ndarray = require('ndarray');
var A = ndarray(new Float64Array([...]), [m, n]);
var results = qr.decompose(A);
var Q = results.q;
var 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:
+- -+ +- -+ +- -+
| c -s | | a | = | r |
| s 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:
+- -+
G(i,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:
R = A
Q = I
for j = 0 : 1 : n-1
for 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^T * R
Q = Q * G
end
end
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:
+- -+
| * |
| 7 * |
| 6 8 * |
| 5 7 9 * |
| 4 6 8 10 * |
| 3 5 7 9 11 * |
| 2 4 6 8 10 12 * |
| 1 3 5 7 9 11 13 * |
+- -+
License
© 2015 Scijs. MIT License.
Author
Tim Bright