Exponential of a matrix

[loiseaujc]

Following a comment by Meromorphic on the discourse here, this PR implements the matrix exponential function expm. It does so in a new stdlib_linalg_matrix_functions submodule which might eventually accomodate later on implementations of other matrix functions (e.g. logm, sqrtm, cosm, etc).

Proposed interface

• E = expm(A [, order, err]) where A is the input n x n matrix, order (optional) is the order of the Pade approximation, and err of type linalg_state_type for 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

• Base implementation.
• Tests
• in-code Documentation
• Specifications
• Example

Possible improvements

The implementation is fully working and validated. Computational performances and/or robustness could potentially be improved eventually by considering:

• Automatic estimation of the best order for the Pade approximation based on the data. This is used for instance in scipy.linalg.expm.
• Balancing is not used currently, but it might improve the robustness.
• At the end of the algorithm, the squaring step is implemented using a simple loop involving matmul. Additional performances can be gained by implementing a fast matrix_power function (see np.matrix_power for instance). (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)

In practice, it has however never failed me so far.

cc @perazz, @jvdp1, @jalvesz

[loiseaujc]

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.

[jvdp1]

Thank you @loiseaujc ! LGTM. I have some minor suggestions and some additional ones that might lead to some discussions.

[jvdp1]

Should this procedure not provided by stdlib_lapack (or similar)? Otherwise, developers and users might replicate this type of procedures quite often (cc @perazz )

[loiseaujc]

I think it should. I've sent a PM on the discourse to discuss this issue actually. I was suggesting to create a new stdlib_lapack_handling module where all of these functions could be defined and accessed by other modules.