### 2 releases

0.1.1 | Jan 22, 2019 |
---|---|

0.1.0 | Jan 22, 2019 |

#**4** in #exponential

**MIT/Apache**

31KB

568 lines

# Matrix exponentiation in Rust

This crate contains

, 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.`expm`

It uses the excellent

crate for matrix storage.`rust-ndarray`

## 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

.
See the explanations given at the `openblas_src`

organization.`blas-lapack-rs`

`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`,` `1``f64``.``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

struct are done, and `Expm`

can be reused to repeatedly
calculate the exponential of matrices with the same dimension.`Expm`

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

and`t`

when calculating the 1-norms (right now,`itmax`

,`t``=``2`

, which is what Scipy and Matlab are doing);`itmax``=``5` - 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.

#### Dependencies

~5MB

~153K SLoC