update matrix exponential to use pade approximation algorithm for Issue #380
pkerichang
Branch: pkerichang/blaze:expm-fix
pkerichang/blaze:expm-fix
Branch: blaze-lib/blaze:master
blaze-lib/blaze:master
Open
#49 · Created · Last updated
Description
Hello,
I see that you already assigned someone to worked on Issue #380, but I have also started working on this problem on my own (mainly because I need this bug resolved soon). I basically just ported scipy’s expm() implementation here:
there are still some room for optimization, but at least the code works for now. Let me know if you’re interested in this (I understand if for some reason you need to keep the old algorithm. In that case I’ll probably just keep this alternative implementation on a fork).
Notes/summary of changes:
I updated the matrix exponential unit test to test for the case shown in Issue #380, and some other test cases from the paper “A New Scaling and Squaring Algorithm for the Matrix Exponential.” I also refactored the unit test code using templates, because the old unit tests have two copies of every test case where the only difference is the storage order, which is difficult to read/maintain.
I tried very hard to keep adapter types around in my code so matrix multiplication can benefit from speed optimizations blaze has built in. This is my first time dealing with adapters, so if I’m doing something bad please let me know.
Because of Issue #382, I have to store a (probably dense) identity matrix. once that is resolve we can have some memory improvement.
Compared to scipy’s implementation, the use of approximate induced-1-norms for determining approximation methods is not implemented yet, since they’re just speed optimization and won’t impact behavior.
One drawback of this algorithm versus the old one is that solve() is used to compute the final matrix. So hypothetically there could be some input matrices that work with the old algorithm and not this algorithm. Although in practice I never see that happen.
Scipy’s expm() function works with sparse matrices, so there’s no reason why this new ported algorithm can’t work on sparse matrices too. I just don’t have experience coding expressions for sparse matrices that I don’t know where to start.
Please let me know your comments. Obviously there’s quite a bit of code so I understand if you don’t have time immediately. Thanks.
Hello,
I see that you already assigned someone to worked on Issue #380, but I have also started working on this problem on my own (mainly because I need this bug resolved soon). I basically just ported scipy’s
expm()implementation here:https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.expm.html#scipy.sparse.linalg.expm
there are still some room for optimization, but at least the code works for now. Let me know if you’re interested in this (I understand if for some reason you need to keep the old algorithm. In that case I’ll probably just keep this alternative implementation on a fork).
Notes/summary of changes:
I updated the matrix exponential unit test to test for the case shown in Issue #380, and some other test cases from the paper “A New Scaling and Squaring Algorithm for the Matrix Exponential.” I also refactored the unit test code using templates, because the old unit tests have two copies of every test case where the only difference is the storage order, which is difficult to read/maintain.
I tried very hard to keep adapter types around in my code so matrix multiplication can benefit from speed optimizations blaze has built in. This is my first time dealing with adapters, so if I’m doing something bad please let me know.
Because of Issue #382, I have to store a (probably dense) identity matrix. once that is resolve we can have some memory improvement.
Compared to scipy’s implementation, the use of approximate induced-1-norms for determining approximation methods is not implemented yet, since they’re just speed optimization and won’t impact behavior.
One drawback of this algorithm versus the old one is that
solve()is used to compute the final matrix. So hypothetically there could be some input matrices that work with the old algorithm and not this algorithm. Although in practice I never see that happen.Scipy’s
expm()function works with sparse matrices, so there’s no reason why this new ported algorithm can’t work on sparse matrices too. I just don’t have experience coding expressions for sparse matrices that I don’t know where to start.Please let me know your comments. Obviously there’s quite a bit of code so I understand if you don’t have time immediately. Thanks.
EDIT: implemented code fragment 2.1