A module for calculating the in-place QR decomposition of an ndarray using Householder triangularization
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).
For an example, see ndarray-householder-qr. The only difference for the complex extension is that all arrays are doubled up with real and imaginary components in separate arrays. For some hints on working with complex numbers in ndarrays, see ndarray-blas-level1-complex.
factor( A_r, A_i, d_r, d_i )
Computes the in-place triangularization of
A_i (the real and imaginary components 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_i must be one-dimensional vectors with length at least
multByQ( A_r, A_i, x_r, x_i )
Compute the product Q * x in-place, replacing x with Q * x.
A_i are the real/imag components of the in-place factored matrix.
multByQinv( A_r, A_i, x_r, x_i )
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.
constructQ( A_r, A_i, Q_r, Q_i )
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.
factor( A_r, A_i, Q_r, Q_i )
Incomplete Compute the in-place QR factorization of A, storing R in A and outputting Q in Q.
solve( A_r, A_i, d_r, d_i, x_r, x_i )
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_iare the real/imag components of the in-place factored matrix computed by
d_iare the real/imag components of the diagonal of
x_iare the real/imag components of the input vector of length m. The answer is computed in-place in the first n entries of
x. The remaining entries are zero.
(c) 2015 Ricky Reusser. MIT License