ndarray-householder-qr

A module for calculating the in-place QR decomposition of an ndarray using Householder triangularization

Introduction

The algorithm is the Householder QR Factorization algorithm as found on p. 73 of Trefethen and Bau's Numerical Linear Algebra. In pseudocode, the algorithm is:

for k = 1 to n
  x = A[k:m,k]
  v_k = sign(x_1) ||x||_2 e_1 + x
  v_k = v_k / ||v_k||_2
  A[k:m,k:n] = A[k:m,k:n] - 2 v_k (v_k^* A[k:m,k:n])

The specific implementation is based on the pseudocode from Walter Gander's Algorithms for the QR-Decomposition. This algorithm computes both R and the Householder reflectors in place, storing R in the upper-triangular portion of A, the diagonal of R in a separate vector and the Householder reflectors in the columns of A. To eliminate unnecessary operations, the Householder reflectors are normalized so that norm(v) = sqrt(2).

Example

A straightforward example of the usefulness of QR factorization is the solution of least squares problems. To fit the model y = a0 * x + a1 to the data points [x1,y1] = [0,1], [x2,y2] = [1,2], [x3,y3] = [2,3]:

var qr = require('ndarray-householder-qr'),
var vander = require('ndarray-vandermonde'),
    
var m = 3,
var n = 2,
 
var x = ndarray([0,1,2]), // independent variable
var y = ndarray([1,2,3]), // data points
 
var d = pool.zeros([n]),
var A = vander(x, n);
 
qr.factor(A, d);
qr.solve(A, d, y);
 
// result: y = ndarray([1, 1, 0]) --> y = 1 * x + 1

After this calculation, the factorization can be reused to solve for other inputs:

var y2 = ndarray([2, 3, 4]);
 
qr.solve(A, d, y2);
 
// result: y = ndarray([2, 1, 0]) --> y = 1 * x + 2

Usage

`qr.factor( A, d )`

Computes the in-place triangularization of A, returning the Householder reflectors in the lower-triangular portion of A (including the diagonal) and R in the upper-triangular portion of A (excluding diagonal) with the diagonal of R stored in d. d must be a one-dimensional vector with length at least n.

`qr.multiplyByQ(A, x)`

Compute the product Q * x in-place, replacing x with Q * x. A is the in-place factored matrix.

`qr.multiplyByQinv(A, x)`

A is the in-place factored matrix. Compute the product Q^-1 * x in-place, replacing x with Q^-1 * x. Since the product is shorter than x for m > n, the entries of x from n+1 to m will be zero.

`qr.constructQ(A, Q)`

Given the in-place factored matrix A (diagonal not necessary), construct the matrix Q by applying the reflectors to a sequence of unit vectors. The dimensions of Q must be between m x n and m x m. When the dimensions of Q are m x n, Q corresponds to the Reduced QR Factorization. When the dimensions are m x m, Q corresponds to the Full QR Factorization.

`qr.solve(A, d, x)`

Use the previously-calculated triangularization to find the vector x that minimizes the L-2 norm of (Ax - b). Note that the vector b is modified in the process.

A is the in-place factored matrix computed by factor
d is the diagonal of R computed by factor
x is the input vector of length m. The answer is computed in-place in the first n entries of x. The remaining entries are zero.

Benchmarks

$ npm run bench