Skip to content

Matrix exponentials in Rust

License

Apache-2.0, MIT licenses found

Licenses found

Apache-2.0
LICENSE-APACHE
MIT
LICENSE-MIT
Notifications You must be signed in to change notification settings

SuperFluffy/rust-expm

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

4 Commits
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Matrix exponentiation in Rust

This crate contains expm, an implementation of Algorithm 6.1 by Al-Mohy, Higham in the Rust programming language. It calculates the exponential of a matrix. See the linked paper for more information.

It uses the excellent rust-ndarray crate for matrix storage.

Example usage

The example below calculates the exponential of the unit matrix.

Important: You need to explicitly link to a BLAS + LAPACK provider such as openblas_src. See the explanations given at the blas-lapack-rs organization.

extern crate openblas_src;
use approx::assert_ulps_eq;
#[test]
fn exp_of_unit() {
    let n = 5;
    let a = ndarray::Array2::eye(n);
    let mut b = unsafe { ndarray::Array2::<f64>::uninitialized((n, n)) };

    crate::expm(&a, &mut b);

    for &elem in &b.diag() {
        assert_ulps_eq!(elem, 1f64.exp(), max_ulps=1);
    }
}

TODO

Care was taken to implement the algorithm with performance in mind. As such, no extra allocations after the initial setup of the Expm struct are done, and Expm can be reused to repeatedly calculate the exponential of matrices with the same dimension.

However, profiling the code might reveal ways to improve it.

  • Profile code and optimize;
  • Write more tests to verify the goodness of the algorithm for difficult matrices;
  • Tune the parameters t and itmax when calculating the 1-norms (right now, t=2, itmax=5, which is what Scipy and Matlab are doing);
  • Ensure that the compiler performs const propagation when calculating the Padé coefficients;
  • Test the Padé coefficient against their hard-coded results;
  • Evaluate, whether the unsafe blocks are really necessary in all instances or could be removed.

About

Matrix exponentials in Rust

Resources

License

Apache-2.0, MIT licenses found

Licenses found

Apache-2.0
LICENSE-APACHE
MIT
LICENSE-MIT

Stars

Watchers

Forks

Packages

No packages published

Languages