Struct rulinalg::matrix::decomposition::PartialPivLu

pub struct PartialPivLu<T> { /* fields omitted */ }

LU decomposition with partial pivoting.

For any square matrix A, there exist a permutation matrix P, a lower triangular matrix L and an upper triangular matrix U such that

PA = LU.

However, due to the way partial pivoting algorithms work, LU decomposition with partial pivoting is in general only numerically stable for well-conditioned invertible matrices.

That said, partial pivoting is sufficient in the vast majority of practical applications, and it is also the fastest of the pivoting schemes in existence.

Applications

Given a matrix x, computing the LU(P) decomposition is simple:

use rulinalg::matrix::decomposition::{PartialPivLu, LUP, Decomposition};
use rulinalg::matrix::Matrix;

let x = Matrix::<f64>::identity(4);

// The matrix is consumed and its memory
// re-purposed for the decomposition
let lu = PartialPivLu::decompose(x).expect("Matrix is invertible.");

// See below for applications
// ...

// The factors L, U and P can be obtained by unpacking the
// decomposition, for example by destructuring as seen here
let LUP { l, u, p } = lu.unpack();

Solving linear systems

Arguably the most common use case of LU decomposition is the computation of solutions to (multiple) linear systems that share the same coefficient matrix.

let b = vector![3.0, 4.0, 2.0, 1.0];
let y = lu.solve(b)
          .expect("Matrix is invertible.");
assert_vector_eq!(y, vector![3.0, 4.0, 2.0, 1.0], comp = float);

// We can efficiently solve multiple such systems
let c = vector![0.0, 0.0, 0.0, 0.0];
let z = lu.solve(c).unwrap();
assert_vector_eq!(z, vector![0.0, 0.0, 0.0, 0.0], comp = float);

Computing the inverse of a matrix

The LU decomposition provides a convenient way to obtain the inverse of the decomposed matrix. However, please keep in mind that explicitly computing the inverse of a matrix is usually a bad idea. In many cases, one might instead simply solve multiple systems using solve.

For example, a common misconception is that when one needs to solve multiple linear systems Ax = b for different b, one should pre-compute the inverse of the matrix for efficiency. In fact, this is practically never a good idea! A far more efficient and accurate method is to perform the LU decomposition once, and then solve each system as shown in the examples of the previous subsection.

That said, there are definitely cases where an explicit inverse is needed. In these cases, the inverse can easily be obtained through the inverse() method.

Computing the determinant of a matrix

Once the LU decomposition has been obtained, computing the determinant of the decomposed matrix is a very cheap operation.

assert_eq!(lu.det(), 1.0);

Methods

impl<T: 'static + Float> PartialPivLu<T>

fn decompose(matrix: Matrix<T>) -> Result<Self, Error>

Performs the decomposition.

Panics

The matrix must be square.

Errors

An error will be returned if the matrix is singular to working precision (badly conditioned).

impl<T> PartialPivLu<T> where T: Any + Float

fn solve(&self, b: Vector<T>) -> Result<Vector<T>, Error>

Solves the linear system Ax = b.

Here, A is the decomposed matrix satisfying PA = LU. Note that this method is particularly well suited to solving multiple such linear systems involving the same A but different b.

Errors

If the matrix is very ill-conditioned, the function might fail to obtain the solution to the system.

Panics

The right-hand side vector b must have compatible size.

Examples

let x = Matrix::identity(4);
let lu = PartialPivLu::decompose(x).unwrap();
let b = vector![3.0, 4.0, 2.0, 1.0];
let y = lu.solve(b)
          .expect("Matrix is invertible.");
assert_vector_eq!(y, vector![3.0, 4.0, 2.0, 1.0], comp = float);

fn inverse(&self) -> Result<Matrix<T>, Error>

Computes the inverse of the matrix which this LUP decomposition represents.

Errors

The inversion might fail if the matrix is very ill-conditioned.

fn det(&self) -> T

Computes the determinant of the decomposed matrix.

Trait Implementations

impl<T: Debug> Debug for PartialPivLu<T>

fn fmt(&self, __arg_0: &mut Formatter) -> Result

Formats the value using the given formatter.

impl<T: Clone> Clone for PartialPivLu<T>

fn clone(&self) -> PartialPivLu<T>

Returns a copy of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<T: Clone + One + Zero> Decomposition for PartialPivLu<T>

type Factors = LUP<T>

The type representing the ordered set of factors that when multiplied yields the decomposed matrix.

fn unpack(self) -> LUP<T>

Extract the individual factors from this decomposition.