Exponential of a matrix #968
Conversation
|
I think it is pretty much ready for review. Note that the optimizations I've mentioned (i.e. variable order for the Pade approximation, balancing, etc) wouldn't change the signature of the function. Either we take it as is and gradually improve it over time, or we improve it right away (but I won't have much time to do it right now). As you like. |
There was a problem hiding this comment.
Thank you @loiseaujc ! LGTM. I have some minor suggestions and some additional ones that might lead to some discussions.
Co-authored-by: Jeremie Vandenplas <[email protected]>
Update working branch with latest master features.
This reverts commit 65ad5f2.
|
Seems like I have issues with |
|
|
||
| contains | ||
|
|
||
| !> schur decomposition tests |
There was a problem hiding this comment.
I think this is a copy-paste typo: "Matrix exponential tests" ?
There was a problem hiding this comment.
Yes indeed. Will update shortly.
| contains | ||
|
|
||
| #:for rk,rt,ri in RC_KINDS_TYPES | ||
| module function stdlib_expm_${ri}$(A, order, err) result(E) |
There was a problem hiding this comment.
For the base implementation I would have done it using a subroutine that pre-supposes the matrix is already allocated in memory, e.g.: real(<>), intent(inout) :: E(:,:). The function interface could be simply built on top of it to have automatic (costly) allocations.
There was a problem hiding this comment.
I mimicked what had been for the Cholesky factorization. The lowest-level driver is a subroutine doing in-place computation. I've also added the matrix_exp interface if anyone wants to do a subroutine call rather than a function one.
| allocate(tests(0)) | ||
|
|
||
| #:for rk,rt,ri in RC_KINDS_TYPES | ||
| tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)] |
There was a problem hiding this comment.
Maybe the current fails with msys CI is related to the changes in #1008
There was a problem hiding this comment.
I merged the latest additions from master into my branch but it does not seem to solve the problem. Really not sure what's going on since all other tests pass with flying colors. Will investigate further.
There was a problem hiding this comment.
Looking a bit online, it seems that the exit code 0xC0000374 returned when the test fails indicates a heap corruption problem.
Following a comment by Meromorphic on the discourse here, this PR implements the matrix exponential function
expm. It does so in a newstdlib_linalg_matrix_functionssubmodule which might eventually accomodate later on implementations of other matrix functions (e.g.logm,sqrtm,cosm, etc).Proposed interface
E = expm(A [, order, err])whereAis the inputn x nmatrix,order(optional) is the order of the Pade approximation, anderrof typelinalg_state_typefor error handling.Key facts
The implementation makes use of a standard "scaling and squaring" approach to compute the Pade approximation with a fixed order. It is adapted from the implementation by John Burkardt.
Progress
Possible improvements
The implementation is fully working and validated. Computational performances and/or robustness could potentially be improved eventually by considering:
orderfor the Pade approximation based on the data. This is used for instance inscipy.linalg.expm.At the end of the algorithm, the squaring step is implemented using a simple loop involving(irrelevant for this particular task since the scaling factor is chosen as a power of 2, might still be a useful function in the grand scheme of things)matmul. Additional performances can be gained by implementing a fastmatrix_powerfunction (seenp.matrix_powerfor instance).In practice, it has however never failed me so far.
cc @perazz, @jvdp1, @jalvesz