Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added LU decomposition with complete pivoting #160

Merged
merged 13 commits into from Mar 24, 2017
Merged

Added LU decomposition with complete pivoting #160

Changes from 1 commit
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
59c268f
Added a "get" method to BaseMatrix.
Feb 25, 2017
4244708
Added a "get_mut" function to BaseMatrixMut.
Feb 25, 2017
9376d7d
Changed indentation and comments.
Feb 27, 2017
daaf922
Merge remote-tracking branch 'upstream/master'
Mar 7, 2017
f0be956
Added LU decomposition with complete pivoting, resolving #98
Mar 9, 2017
7d0c67b
Changed FullPivLu to only work with square matrices
Mar 12, 2017
2ffb026
Added FullPivLu to the decomposition table
Mar 12, 2017
572d45f
Fixed a warning I overlooked in a test
Mar 12, 2017
aeeb09c
Fix issues brought up in PR
Mar 18, 2017
ad9cdc4
Fix the tolerance for determining rank
Mar 21, 2017
e340619
Change is_singular to is_invertible, other minor changes
Mar 22, 2017
337405e
Fix bad 'is_invertible' behavior for empty matrix
Mar 23, 2017
5b54687
Fix a dumb typo in a comment
Mar 23, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 62 additions & 5 deletions src/matrix/decomposition/lu.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -477,8 +477,12 @@ impl<T> FullPivLu<T> where T: Any + Float {
let mut inv = Matrix::zeros(n, n);
let mut e = Vector::zeros(n);

// "solve" will return an error if the matrix is
// singular, so no need to check for singularity here
if self.is_singular() {
return Err(
Error::new(
ErrorKind::DivByZero,
"Singular matrix found while attempting inversion."));
}

for i in 0 .. n {
e[i] = T::one();
Expand Down Expand Up @@ -514,10 +518,63 @@ impl<T> FullPivLu<T> where T: Any + Float {
}

/// Computes the rank of the decomposed matrix.
/// # Examples
///
/// ```
/// # #[macro_use] extern crate rulinalg;
/// # use rulinalg::matrix::decomposition::FullPivLu;
/// # use rulinalg::matrix::Matrix;
/// # fn main() {
/// let x = matrix![1.0, 2.0, 3.0;
/// 4.0, 5.0, 6.0;
/// 5.0, 7.0, 9.0];
/// let lu = FullPivLu::decompose(x).unwrap();
/// assert_eq!(lu.rank(), 2);
/// # }
/// ```
pub fn rank(&self) -> usize {
self.lu.diag().fold(
0 as usize,
|x, &y| if y.abs() > T::epsilon() { x + 1 } else { x } )
let eps = self.epsilon();
let mut rank = 0;

for d in self.lu.diag() {
if d.abs() > eps {
rank = rank + 1;
} else {
break;
}
}

rank
}

/// Returns whether the matrix is singular.
///
/// # Examples
///
/// ```
/// # #[macro_use] extern crate rulinalg;
/// # use rulinalg::matrix::decomposition::FullPivLu;
/// # use rulinalg::matrix::Matrix;
/// # fn main() {
/// let x = Matrix::<f64>::identity(4);
/// let lu = FullPivLu::decompose(x).unwrap();
/// assert!(!lu.is_singular());
///
/// let y = matrix![1.0, 2.0, 3.0;
/// 4.0, 5.0, 6.0;
/// 5.0, 7.0, 9.0];
/// let lu = FullPivLu::decompose(y).unwrap();
/// assert!(lu.is_singular());
/// # }
/// ```
pub fn is_singular(&self) -> bool {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think perhaps the name is_singular is a little unfortunate if we wish to move forward with our plans of extending LU to rectangular matrices. In general, one only talks about singular matrices if they are also square. Hence, it would be a little strange to refer to a rectangular matrix as 'singular'. More precisely, is a rectangular matrix of full row/col rank singular, or non-singular? But if it's non-singular it's still not invertible!

The benefit of flipping it around and instead using is_invertible is that a matrix is invertible if and only if it is square and has no zeros on the diagonal of U. Thus it extends also nicely to rectangular matrices since they can never be invertible (and will always return false).

Did this make any sense?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, a relatively minor thing, but it is sufficient to check the last diagonal of U:

// Assuming is_invertible(...)
if diag_size > 0 {
    let diag_last = diag_size - 1;
    self.lu[[diag_last, diag_last]] > self.epsilon()
} else {
    // Are empty matrices invertible or not? Well, rank(A) = n, if A is nxn with n = 0, so I'd say it's invertible!
    // This seems to be the most consistent choice. It would be good to have a comment in the code about
    // this to make it clear that this is a deliberate choice.
    true
}

let diag_size = cmp::min(self.lu.rows(), self.lu.cols());

self.rank() != diag_size
}

fn epsilon(&self) -> T {
self.lu[[0, 0]].abs() * T::epsilon()
Copy link
Collaborator

@Andlon Andlon Mar 22, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We actually need to take into account the corner case of empty matrices here (in which case accessing [0, 0] will panic). The easiest is probably to use i.e. diag:

let max_element = self.lu.diag().next().unwrap_or(T::one());
max_element.abs() * T::epsilon()

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, the simplest is probably to use get! I.e. self.lu.get(0, 0).unwrap_or(T::one());.

}
}

Expand Down