Complex matrices

[TrueDoctor]

Do you plan to implement complex matrices?

[SuperFluffy]

I don't recall if the paper that this implementation is based on only defines expm for reals, or if one can just drop in complex. Looking at scipy.linalg.sparse.expm, they don't seem to make any distinction.

I am happy to accept (and review) a implementation for complex numbers, but I don't have the capacity to provide one myself (or do research in how far this algorithm is applicable).

[TrueDoctor]

Thanks a lot, I already started testing complex numbers, but it seem like rust-condest has to be extended as well. I might add support for f32 and c32 as well…
I'll keep you updated if I make any progress

[SuperFluffy]

Excellent. Good luck! :)

Just a word of caution: before you use this library for anything serious, cross check your results with matlab. Matlab's implementation was provided by Al-Mohy and Higham, so that's the one I have most trust in.

When I was writing rust-expm I discovered that scipy's expm was not true to the paper and there was an error in the calculation of the backward error function, see my PR here:
scipy/scipy#9705

Scipy's expm didn't seem to outright produce any wrong results; but there is probably some edge case where it might have done so.

What I am trying to say is: I don't have a sufficiently well developed mathematical understanding to really grok their paper and come up with a test suite that verifies my implementation, in line with the error bounds they provide in the paper.

So while I think I didn't mess up rust-expm, it's not really tested...

[TrueDoctor]

I only managed to refactor about half of the code with generics because I just didn't have the time. For now, I'll just call scipy in my rust code, because it gets the job done and it is not the performance critical part. I'll let you know, should I pick up work again, but until then, it might be appropriate to close the issue

[TrueDoctor]

Uh I just saw https://github.com/hmunozb/rust-expm/tree/generic-scalar That looks promising! Should have looked at the forks earlier 🙄
It would be great if could integrate the changes into rust-expm and rust-condest to help others avoid the same pitfall as I

[SuperFluffy]

Thanks for the hint, I haven't even seen.

If the maintainer of the fork wants to submit a PR, I will gladly incorporate it. Looks like a lot of changes were done to support their use on MacOS though (changing to lapacke, for example).

[TrueDoctor]

As far as I know, all of the changes on the generic-scalar branch work without changing the backend. I think lapack is only used on the with-lapack branch.
But the crate does not seem to be stable yet, I ran into an issue while trying it out today.

(1, 0)                  
(0, 1)

works, but

(1, 0)
(0, 2)

does not. I don't know whether that is an issue on my side, but I've written a mail to the maintainer since the fork has issues disabled.
I'll let you know when I hear back from him

[hmunozb]

Hi everyone. As it appears, my forked repository in theory enables complex numbers, but I put off a PR here for now because it relies on a uniform interface for blas/lapack routines generic over f32/f64/c32/c64 scalars, which I whipped up in a frenzy, so it's unstable and not yet listed on crates.io. I also went to more lengths to support running on macOS. (The lapacke dependency stayed the same in the final generic-scalar branch, but I found it more feasible to link to system openblas as a dev dependency during testing)

As for @TrueDoctor's issue, which concerns the panic message

'tests::complex_exp' panicked at 'assertion failed: `(left == right)`
 left: `4`,
right: `7`', src/lib.rs:249:9

I've only ran into it when the matrix exponential was attempted on matrices with fairly large norms, I assumed this was just an upstream issue with the Pade sums, but I've never seen it occur with such a simple matrix. Still, I can offer some time to look into it if this is an issue with c64 that doesn't occur with f64 for that example.

I've enabled issues on my expm and condest forks (I wasn't aware they were disabled). My intention was to stabilize my blas-traits crate and publish it once I was more or less done with the research project on my end, and then do PRs here. For the most part, the only poor folks with the misfortune to need complex exponentials are those of us in the quantum mechanics/quantum computing fields, so I'm all for helping make sure things are stable in rust for any increasing interest in that community.

[TrueDoctor]

This seems to be a upstream issue: #4
It seems to work correctly by just removing the assertion…

[emmatyping]

I'd love to see generic (well, specifically complex) support added. Is there anything I can do to help move this along?

[TrueDoctor]

I'd love to see generic (well, specifically complex) support added. Is there anything I can do to help move this along?

@hmunozb has already put an amazing effort into porting the function for complex numbers.
SuperFluffy/rust-condest#7 seems to be a prerequisite
Reviewing the pr would probably help @SuperFluffy with merging the code.

I think #4 is still open but I haven't checked that in a while. But yeah the best way to help would probably be to help @hmunozb and maybe also do some testing to confirm the impl produces the expected results

[hmunozb]

Yes, an extra pair of eyes on SuperFluffy/rust-condest#7 would be helpful. Once that gets merge then I'll make a PR for the generic-scalar branch right away.

[emmatyping]

Okay I'll try to review that PR later!

[emmatyping]

Alright, I gave it a look over. Happy to review the generic-scalar PR when that goes up! Thank you for working on these :)