Matrix exponential?

[ChristopherRabotin]

Is there a way to compute a matrix exponential using this library? If not, would that be of interest to anyone but myself? If there is some interest, I wouldn't mind giving a shot at implementing it (I used to have a reference book on implementing linear algebra computations, I just need to find it). However, I'm still a beginner in Rust, so there might need to be a thorough review process.

[sebcrozet]

Matrix exponential is not yet implemented but would be a great addition. If you which to try implementing it, it would be very appreciated. If you are a beginner, perhaps it would be easier to start implementing it for dynamically-sized real matrices first:

impl<N: Real> DMatrix<N> {
    fn exp(&self) -> DMatrix<N> {
        ...
    }
}

That way you wont have to struggle with the trait bounds that might be a bit intricate to handle the general case (the case where the input matrix can be a statically-sized matrix as well as a dynamically-sized one).

[ChristopherRabotin]

I'll be implementing it using the algorithm in "Matrix Computations 4th Edition" by Golub and Van Loan, cf. https://www.cs.cornell.edu/cv/GVL4/golubandvanloan.htm .

Edit: I will be using the Schur decomposition algorithm (I will detail the exact algorithm in the code). This was the one I initially wanted to implement but I couldn't find the reference after days of searching. Anyhow, I just found it: "Functions of Matrices: Theory and Computation" by Nicholas J. Higham, 2008 (for some odd reason I spent several weeks in December reading up on linear algebra algorithms). The attached screenshot is from "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later∗" page 29 (2003) by Moler and Van Loan.

[ChristopherRabotin]

I was going to put this in the core.ops.rs file: is that okay, or is there a better location for this function?

[sebcrozet]

ops.rs is ok as long as the implementation is not too long (say, less than 30 lines or so). Otherwise it is preferable to write it on its own file inside of the linalg folder.

[ChristopherRabotin]

Sorry, I'm having some issues understanding the current code, and that's surely due to my level of knowledge of Rust.

How does the following code work, from matrix_alga.rs? How does this not lead to an infinite loop? (For context, I was looking to see if the infinity L-norm for the matrix was already implemented.)

#[inline]
    fn norm(&self) -> N {
        self.norm()
    }

Thanks

[Ralith]

That function definition exists inside a trait impl, so it's defining a trait method. Inherent methods take precedence over trait methods with the same name, so it's calling the fn norm defined directly (i.e. not by a trait) on MatrixMN<N, R, C>.

[ChristopherRabotin]

I have yet another dumb question. I'm trying to run the tests (since I'm writing tests at the same time as I'm coding this up (the hardest so far is to get a good grasp of how to use this library from my small Rust knowledge)), and I get the following error:

> cargo test
   Compiling nalgebra v0.13.0 (file:///home/chris/Workspace/rust/nalgebra)
error[E0463]: can't find crate for `abomonation`
 --> tests/lib.rs:8:1
  |
8 | extern crate abomonation;
  | ^^^^^^^^^^^^^^^^^^^^^^^^^ can't find crate

error: aborting due to previous error

error: Could not compile `nalgebra`.
warning: build failed, waiting for other jobs to finish...
error: build failed

Moreover, is there a way to run a specific test from the test suite instead of the whole thing? I'd like to only tests my exponential code. Thanks.

[WaDelma]

You can run specific test by cargo test your_test_name

That extern crate abomonation; should be shielded by #[cfg(feature = "abomonation-serialize")] instead of #[cfg(feature = "abomonation")]. Cargo doesn't seem complain about it because the optional dependencies are also features...

[ChristopherRabotin]

I'm running version 1.20 of rust. Could that explain why the cfg shield doesn't seem to be working? Or, how do I enable or disable a cfg flag? Maybe did I enabled it by accident?

[ChristopherRabotin]

Re: cfg abomonation issue

I still have not been able to run the tests correctly. If I enable all the features with the --all-features flag, the tests actually fail (using rust 1.20). I also have not been able to disable that feature in the library. Hence, I'm currently testing the code from another project and writing all the implementation seperatly and will copy over that code to nalgebra when ready.

How do I use the Real trait?

I am trying to do some basic math with the elements of a matrix. However, rustc is complaining about me trying to divide by a floating point value instead of a parameter. How can I fix this?

fn exp_upper<N: Real>(mat: DMatrix<N>) -> DMatrix<N> {
    let n = mat.shape().0;
    let l1 = mat[(0, 0)];
    let l2 = mat[(1, 1)];
    let t12 = mat[(1, 0)];
    let delta = (l1 - l2) / 2.0;
    let e_off_diag = t12 * ((l1 + l2) / 2.0).exp() * delta.sinh() / delta;
    DMatrix::from_row_slice(n, n, &[l1.exp(), e_off_diag, l2.exp(), 0.0])
}

fn main() {

    let l1: f64 = rand::random();
    let l2: f64 = rand::random();
    //let mat = Matrix2::new(l1, 0.0, 0.0, l2);
    let mat = DMatrix::from_row_slice(2, 2, &[l1, 0.0, l2, 0.0]);

    let rslt;
    if mat.relative_eq(&mat.upper_triangle(), 1e-32, 1e-32) {
        // This is an upper triangular matrix, exponential computation is trivial.
        rslt = exp_upper(mat);
    } else if mat.relative_eq(&mat.lower_triangle(), 1e-32, 1e-32) {
        // Lower triangular, let's compute the exponential on the transpose.
        rslt = exp_upper(mat.transpose()).transpose();
    }
    print!("{:?}", mat);
    println!("{:?}", rslt);
}

[WaDelma]

I mean that there is currently bug in nalgebra:

#[cfg(feature = "abomonation")]
extern crate abomonation;

should be changed to

#[cfg(feature = "abomonation-serialize")]
extern crate abomonation;

Because the you are using generic real you cannot use concrete literals in your code. You can use to cast to it from literals. You can use N::one() and N::zero() to get one and zero of type N. You can also do let two: N = ::convert(2.0); to convert the literal.

[milibopp]

Uhm, that abomonation issue is my fault. Will fix the bug asap.

[sebcrozet]

@ChristopherRabotin Besides the issue addressed by #287 , the test suite won't work as it with a cargo test. Other feature flags are neccessary at various places of the test suite to fix this completely. Could you try the following:

cargo test --features "debug arbitrary serde-serialize abomonation-serialize"

And to run one specific test, e.g., tests with a name containing exp:

cargo test exp --features "debug arbitrary serde-serialize abomonation-serialize"

[ChristopherRabotin]

That works, thanks @sebcrozet .

[ChristopherRabotin]

Is there anyway to speed up the test building? I'm only getting back to coding this up now as this past week was incredibly busy, and even running just my two tests takes about three minutes.

[sebcrozet]

@ChristopherRabotin One option to speed up test building would be to split the existing tests into multiple crates (by putting severar .rs files at the root of the tests directory. But this will only speed up the compilation of tests when the tests themselves are modified. If nalgebra is modified, all the test will be recompiled anyway.

One thing I've been doing is either write tests on a separate project (on the same cargo workspace), or comment out any test (by commenting the mod statement on, e.g., tests/lib.rs) that is not necessary for my current developments. I don't know if there are better ways.

[harrysarson]

May I ask if there was any progress with this function?

[ChristopherRabotin]

I haven't done any work on it since the October 2017. I'm happy to provide some references if you'd like, but I'm currently too busy with other projects to work on this feature.

[harrysarson]

My matrix was symmetric so I used SymmetricEigen which worked nicely.

[fredrik-jansson-se]

I've implemented the exp function, heavily influenced by SciPy's implementation. It seems to work ok. It's implemented for static size matrices.

Would be very happy to contribute the code.

It could use some love from someone that know trait bounds better than I do and also some more test cases.

Please let me know.

[fredrik-jansson-se]

Hey, now checked in here: https://github.com/fredrik-jansson-se/nalgebra

Happy to create a PR if of interest.

[ChristopherRabotin]

I would be keen to see a PR!

Christopher Rabotin.

On Wed, Apr 1, 2020, 13:03 Fredrik Jansson notifications@github.com wrote:

Hey, now checked in here: https://github.com/fredrik-jansson-se/nalgebra

Happy to create a PR if of interest.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/rustsim/nalgebra/issues/283#issuecomment-607434935, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABEZV2AZQ2U45SGGO73YQX3RKOFYTANCNFSM4DZG4C2Q .

[fredrik-jansson-se]

Done, https://github.com/rustsim/nalgebra/pull/708

[fredrik-jansson-se]

PR is merged since a while, if you're good with it, maybe we can close this issue.

[ChristopherRabotin]

Awesome! Thanks.

I'm happy to close this issue.

Christopher Rabotin.

On Tue, Apr 28, 2020, 09:28 Fredrik Jansson notifications@github.com wrote:

PR is merged since a while, if you're good with it, maybe we can close this issue.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/rustsim/nalgebra/issues/283#issuecomment-620678052, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABEZV2HCGPMPHUR6HF4YWQ3RO3YZXANCNFSM4DZG4C2Q .

[dschrempf]

Hi! Just came across this thread and read a long while until I realized, hey, this should be closed!

Please close!