Expand description
faer 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{thin}}$, then we get $$x = R_{\text{thin}}^{-1} Q_{\text{thin}}^H b$$
to translate this to code, we can proceed as follows:
use faer::dyn_stack::{MemBuffer, MemStack, StackReq};
use faer::reborrow::*;
use faer::linalg::qr::no_pivoting::factor;
use faer::linalg::{householder, triangular_solve};
use faer::{Conj, Mat, Par, mat};
// 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],
];
// approximate solution computed with numpy
let expected_solution = mat![
[0.33960324, -0.33812452, -0.8458301], //
[-0.25718351, 0.6281214, 1.07071764_f64],
];
let rank = Ord::min(a.nrows(), a.ncols());
// we choose the recommended block size for the householder factors of our problem.
let blocksize = factor::recommended_blocksize::<f64>(a.nrows(), a.ncols());
// we allocate the memory for the operations that we perform
let mut mem =
MemBuffer::new(StackReq::any_of(&[
factor::qr_in_place_scratch::<f64>(
a.nrows(),
a.ncols(),
blocksize,
Par::Seq,
Default::default(),
),
householder::apply_block_householder_sequence_transpose_on_the_left_in_place_scratch::<
f64,
>(a.nrows(), blocksize, b.ncols()),
]));
let mut stack = MemStack::new(&mut mem);
let mut qr = a;
let mut h_factor = Mat::zeros(blocksize, rank);
factor::qr_in_place(
qr.as_mut(),
h_factor.as_mut(),
Par::Seq,
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
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(),
Par::Seq,
stack.rb_mut(),
);
solution.truncate(rank, b.ncols());
// compute R_thin^{-1} Q_thin^H×B
triangular_solve::solve_upper_triangular_in_place(
qr.as_ref().split_at_row(rank).0,
solution.as_mut(),
Par::Seq,
);
for i in 0..rank {
for j in 0..b.ncols() {
assert!((solution[(i, j)] - expected_solution[(i, j)]).abs() <= 1e-6);
}
}Modules§
- col_
pivoting - the $QR$ decomposition with column pivoting 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.