Crate faer_qr

This crate provides utilities for computing and manipulating the QR factorization with and without pivoting. The QR factorization decomposes a matrix into a product of a unitary matrix $Q$ (represented using block Householder sequences), and an upper trapezoidal matrix $R$, such that their product is equal to the original matrix (or a column permutation of it in the case where column pivoting is used).

Example

Assume we have an overdetermined system $AX = B$ with full rank, and that we wish to find the solution that minimizes the 2-norm.

This is equivalent to computing a matrix $X$ that minimizes the value $||AX - B||^2$, which is given by the solution $$X = (A^H A)^{-1} A^H B.$$

If we compute the QR decomposition of $A$, such that $A = QR = Q_{\text{thin}} R_{\text{rect}}$, then we get $$X = R_{\text{rect}}^{-1} Q_{\text{thin}}^H B.$$

To translate this to code, we can proceed as follows:

use assert_approx_eq::assert_approx_eq;
use dyn_stack::{PodStack, GlobalPodBuffer, StackReq};
use faer_core::{mat, solve, Conj, Mat, Parallelism};
use reborrow::*;

// we start by defining matrices A and B that define our least-squares problem.
let a = mat![
    [-1.14920683, -1.67950492],
    [-0.93009756, -0.03885086],
    [1.22579735, 0.88489976],
    [0.70698973, 0.38928314],
    [-1.66293762, 0.38123281],
    [0.27639595, -0.32559289],
    [-0.37506387, -0.13180778],
    [-1.20774962, -0.38635657],
    [0.44373549, 0.84397648],
    [-1.96779374, -1.42751757_f64],
];

let b = mat![
    [-0.14689786, -0.52845138, -2.26975669],
    [-1.00844774, -1.38550214, 0.50329459],
    [1.07941646, 0.71514245, -0.73987761],
    [0.1281168, -0.23999022, 1.58776697],
    [-0.49385283, 1.17875407, 2.01019076],
    [0.65117811, -0.60339895, 0.27217694],
    [0.85599951, -0.00699227, 0.93607199],
    [-0.12635444, 0.94945626, 0.86565968],
    [0.02383305, 0.41515805, -1.2816278],
    [0.34158312, -0.07552168, 0.56724015_f64],
];

// computed with numpy
let expected_solution = mat![
    [0.33960324, -0.33812452, -0.8458301],
    [-0.25718351, 0.6281214, 1.07071764_f64],
];

let rank = a.nrows().min(a.ncols());

// we choose the recommended block size for the householder factors of our problem.
let blocksize = faer_qr::no_pivoting::compute::recommended_blocksize::<f64>(a.nrows(), a.ncols());

// we allocate the memory for the operations that we perform
let mut mem =
    GlobalPodBuffer::new(StackReq::any_of(
        [
            faer_qr::no_pivoting::compute::qr_in_place_req::<f64>(
                a.nrows(),
                a.ncols(),
                blocksize,
                Parallelism::None,
                Default::default(),
            )
            .unwrap(),
            faer_core::householder::apply_block_householder_sequence_transpose_on_the_left_in_place_req::<
                f64,
            >(a.nrows(), blocksize, b.ncols())
            .unwrap(),
        ],
    ));
let mut stack = PodStack::new(&mut mem);

let mut qr = a.clone();
let mut h_factor = Mat::zeros(blocksize, rank);
faer_qr::no_pivoting::compute::qr_in_place(
    qr.as_mut(),
    h_factor.as_mut(),
    Parallelism::None,
    stack.rb_mut(),
    Default::default(),
);

// now the Householder bases are in the strictly lower trapezoidal part of `a`, and the
// matrix R is in the upper triangular part of `a`.

let mut solution = b.clone();

// compute Q^H×B
faer_core::householder::apply_block_householder_sequence_transpose_on_the_left_in_place_with_conj(
    qr.as_ref(),
    h_factor.as_ref(),
    Conj::Yes,
    solution.as_mut(),
    Parallelism::None,
    stack.rb_mut(),
);

solution.resize_with(rank, b.ncols(), |_, _| unreachable!());

// compute R_rect^{-1} Q_thin^H×B
solve::solve_upper_triangular_in_place(
    qr.as_ref().split_at_row(rank).0,
    solution.as_mut(),
    Parallelism::None,
);

for i in 0..rank {
    for j in 0..b.ncols() {
        assert_approx_eq!(solution.read(i, j), expected_solution.read(i, j));
    }
}

Modules

• col_pivoting
  The QR decomposition decomposes a matrix $A$ into the product $$AP^T = QR,$$ where $P$ is a permutation matrix, $Q$ is a unitary matrix (represented as a block Householder sequence), and $R$ is an upper trapezoidal matrix.
• no_pivoting
  The QR decomposition decomposes a matrix $A$ into the product $$A = QR,$$ where $Q$ is a unitary matrix (represented as a block Householder sequence), and $R$ is an upper trapezoidal matrix.