Julia: Implement more matrix functions

[jiahao]

Higham and Deadman have published a list of software implementing algorithms for matrix functions, showing that Julia's library for these functions could be improved greatly relative to the state of the art.

The purpose of this issue is to discuss and track the implementation of matrix functions.

From Higham and Deadman's list:

General function evaluations

• [ ] General matrix function with derivatives of the underlying scalar function available: Schur–Parlett algorithm (Davies and Higham, 2003)
• [ ] Function of a symmetric or Hermitian matrix by diagonalization
• [ ] Condition number estimate for general functions: cond(f,A)

• [ ] matrix function - vector product evaluation: f(A) * b

  Specific functions

• [x] expm: scaling and squaring algorithm (Al-Mohy and Higham, 2009)

• [x] logm: inverse scaling and squaring algorithm (Al-Mohy, Higham, and Relton, 2012,
  2013)
• [ ] sqrtm:
  
  
  
  • [x] Schur algorithm (Björck and Hammarling, 1983)
  
  
  • [ ] real version for real matrices (Higham, 1987)
  
  
  • [ ] blocking algorithm for performance improvements (Deadman, Higham, and Ralha, 2013)
  
  
• [x] A^t for real matrix powers: Schur– Padé algorithm (Higham and Lin, 2013) (#21184)

• [ ] Matrix unwinding function

  Fréchet derivatives

• [ ] Exponential: scaling and squaring algorithm (Al-Mohy and Higham, 2009)

• [ ] Logarithm: inverse scaling and squaring algorithm (Al-Mohy, Higham, and Relton, 2013)
• [ ] General matrix function: complex step algorithm (Al-Mohy and Higham, 2010) or using block 2 × 2 matrix formula.

• [ ] A^t for real matrix powers: Schur– Padé algorithm (Higham and Lin, 2013)

  Miscellaneous functions

• [x] inv and det for Tridiagonal and SymTridiagonal matrices (Usmani, 1994) (#5358)

[StefanKarpinski]

It would be amazing to include all of these algorithms. Thanks for putting this issue together, @jiahao.

[jiahao]

It's not clear to me what interface should be implemented for the very generic things like "function of a symmetric or Hermitian matrix".

We could define methods for elementary functions until the cows come home, but that is both tedious and namespace polluting (we'd need sinm, cosm, atanhm, etc. to disambiguate from the ordinary forms which broadcast over the elements) , and furthermore wouldn't scale well for linear combinations of functions.

Maybe we want something like applym, so that a call like applym(A->exp(A)+sin(2A)-cosh(4A), A) would evaluate the matrix expm(A)+sinm(2A)-coshm(4A)?

[StefanKarpinski]

This is one of the reasons why having all these functions be vectorized is kind of shitty. I suspect that @JeffBezanson will vehemently agree. There could be a MatrixFunctions module that you have to import from to get the matrix variants with all the same names.

[jiahao]

Sure, but that won't help with the fact that a complicated matrix function like A->expm(A)+sinm(2A)-coshm(4A) can be computed more quickly all at once than computing (expm(A))+(sinm(2A))-(coshm(4A)) by working out the individual elementary functions separately, since the underlying Schur-Parlett algorithm is essentially the same.

[StefanKarpinski]

Yes, that's a very good point.

[stevengj]

Trefethen's contour-integration method is also neat for mostly analytic functions (with known singularities), though I don't know how it compares to the state of the art (in performance and generality).

[jiahao]

@stevengj it looks like that paper is Ref. 32 of the Higham and Deadman report, and isn't currently available in any major package (although it has Matlab code).

[simonbyrne]

This seems like a good candidate for a GSOC project.

As far as an interface goes, I would suggest a macro

B = @matfun exp(A)+sin(2A)-cosh(4A)

A naive approach could simply replace the functions (e.g. exp -> expm), a more advanced approach could try to identify the common decompositions.

[GunnarFarneback]

The item "[sqrtm] real version for real matrices (Higham, 1987)" is checked. I can't see that algorithm in base. I made an attempt on it about a year ago but concluded that the complexity just wasn't worth it. I've put up the code as a gist if someone wants to have a look.

[jiahao]

Thanks for catching that. I must have ticked it by accident when writing up this issue.

[dlfivefifty]

Just tried to call logm, but wasn't implemented. Had to go to Mathematica

[fabianlischka]

Whoever tackles this issue might find this paper helpful for testing:
Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014

[StefanKarpinski]

Great reference – thanks, @fabianlischka.

[mfasi]

The logm function has been added by pull request #12217.

[ViralBShah]

Also cc: @higham @alanedelman

[ararslan]

Would sinm and cosm (matrix sine and cosine, respectively, see here) be useful? They're easily derived from expm. I can put together a PR if that sounds desirable.

[higham]

Dear Alex,

On 06/05/2016 00:11, Alex Arslan wrote:

Would |sinm| and |cosm| be useful? They're easily derived from |expm|. I can
put together a PR if that sounds desirable.

They are worth having.
It's possible to compute cosm(A) as (X+inv(X))/2, where X = expm(i*A), but
this requires complex arithmetic even when A is real, so is expensive.
It would be better to use the algorithms in

Awad Al-Mohy, Nicholas Higham & Samuel Relton, New Algorithms for
Computing the Matrix Sine and Cosine Separately or Simultaneously,
SIAM J. Sci. Comput. 37, A456-A487, 2015

https://github.com/sdrelton/cosm_sinm

---Nick

[ararslan]

Hi Nick, thanks for the heads up, I'll look at that paper.

While I have you here, would you be willing to license the MATLAB code for the unwinding matrix under MIT so it can be ported to Julia? I wasn't able to find a license for the code that accompanies the paper.

[higham]

Dear Alex,

On 06/05/2016 21:55, Alex Arslan wrote:

While I have you here, would you be willing to license the MATLAB code for the
unwinding matrix under MIT so it can be ported to Julia? I wasn't able to find
a license for the code that accompanies the paper.

Ah, yes - there's a reason why we haven't made that available in the usual
licensed way on GitHub: the code uses subfunction

function [mm,ind] = swapping(m)

copied directly from the MATLAB funm. Since it wasn't an option to write
another version of that subfunction we just made the code available from

http://eprints.ma.man.ac.uk/2094/

I notice that we made the original 2003 funm codes available at

http://www.maths.manchester.ac.uk/~higham/NAMF/

with a GPL license and there is an earlier version of swapping.m there.
I don't remember how that version compared with the one that eventually
appeared in funm. Perhaps you can work from that older version?
If you want to do so and the GPL license is a problem I can change the
license on the above page.

---Nick

[ararslan]

Please correct me if I'm wrong, but my understanding is that MIT and GPL are incompatible; GPL is a much more stringent license that requires that derived works, which includes direct ports of the code, be licensed under GPL as well. Further, if one portion of a project is licensed under GPL, the rest of the project must also be GPL. Since Julia's core is entirely MIT, that wouldn't work; the port would have to be its own GPL-licensed package. So if you're willing to license any of your code under MIT to facilitate porting to Julia, that would be fantastic!

[Keno]

The MIT and the GPL are not incompatible. The result however, would be GPL licensed as a whole distribution. Since the restrictions of the GPL are too stringent for the standard library of a programming language however, we generally do not accept new contributions of GPL licensed code into base julia and are actively trying to move existing GPL licensed code out of base (there's already a makefile flag that disables such code).

[higham]

On 07/05/2016 19:56, Alex Arslan wrote:

So if you're willing to license any of your code under MIT to facilitate
porting to Julia, that would be fantastic!

Done - the page

http://www.maths.manchester.ac.uk/~higham/NAMF/

has been updated with the MIT license.

---Nick

[ViralBShah]

@higham Thank you.

[ararslan]

@higham That's fantastic, thank you very much!

[stevengj]

Note that we are planning on deprecating all of the elementwise "vectorized" operations like exp(A) in favor of exp.(A) in 0.6 (see #17302). So, at some point in the future we may be able to reclaim those names to use for matrix functions.

[dlfivefifty]

👍🏽

[ararslan]

@stevengj I thought the convention for matrix functions was to append m, e.g. expm?

[stevengj]

@ararslan, that's a Matlabism that is imposed on the language because exp(A) is elementwise; the underlying mathematical matrix function has no "m" suffix in the math literature. See #8450 for a long discussion of why/how we moved away from implicitly vectorized functions like exp(A) in Julia.

[ararslan]

I've read through there, I just didn't realize there was any desire to do away with the m carryover from Matlab.

[stevengj]

I don't think we should rush into ditching the m suffix; it's just something we can consider in the future. We will want to have exp(A) give a deprecation message to use exp.(A) for quite a while, to teach users about the dot syntax.

[dlfivefifty]

I think dropping the m is more consistent with the Taylor series definition of exp. For example,

exp(dual(1.,2.)) 

returns the dual corresponding to the Taylor series definition (which in this case is equivalent to expm([1. 2.; 0. 1.]).

That's ignoring the Matlab and Mathematica precedents.

[dlfivefifty]

Whether or not exp(A) replaces expm(A), I think down the line exp.(A) should replace exp(A), as the current exp(A) is confusing for non-MATLAB users

[tkelman]

@dlfivefifty sooner than you think, there's already an open PR for it: https://github.com/JuliaLang/julia/pull/17302/files#diff-8278b779f2ea681192ba5b020a2c3e2bL137

[felixrehren]

Version 2 of the Deadman-Higham list has been out for a year and Julia is the same as it was first time. Let's get better before version 3 comes out. I'm going to look at Schur-Pade -- that together with the top-4 tasks "General function evaluations" in the OP would cover most of this. Anyone else want to join?

P.S. we already beat R in the DH list, with built-in Matlab and SciPy looking vulnerable as well!

[simonbyrne]

Is the "Deadman-Higham ranking" our Weissman Score?

[dlfivefifty]

Could there be a function matrixbroadcast(f,A) that calculates the matrix function f(A) for general f? This would be easy to implement for diagonalisable A, not sure about otherwise

[stevengj]

Even for diagonalizable A, the possibility of nearly defective A would be an issue. Of course, it could be restricted to Hermitian A, but that case is so easy that we would hardly add much value by implementing a function for it in Base.

[alanedelman]

There is a pretty good chance that matrix derivatives would work with ForwardDiff.jl, at least
if there is a corresponding generic linear algebra function. I've been pretty impressed that
one can get derivatives of powers and svd's and inverses, etc.

[ChrisRackauckas]

Should this all be moved out to a MatrixFunctions.jl?

[felixrehren]

Seems like a good idea -- I think everything on the list other than expm logm sqrtm is a clear candidate for moving out (e.g. I have a rootm function that takes arbitrary matrix pth roots -- probably doesn't fit in base, but could be nice somewhere). Whether expm logm sqrtm are worth keeping in base or can move out as well, I don't know.

@ararslan You did the great work moving SpecialFunctions.jl. What do you think about MatrixFunctions.jl? Other people that were involved there IIRC but are not in this thread: @musm, @andreasnoack

[ViralBShah]

Should we close this, and perhaps there ought to be a better place for this in a new repo?

[ChrisRackauckas]

If it's extending sqrt and log like exp in #23233, then it should be in Base since otherwise it's type piracy.

[dlfivefifty]

we could do something like this:

for f in (:sqrt, ...)
   @eval $f(A::AbstractArray) = matrixfunction($f, A)
end
matrixfunction(f, A) = error(“Install MatrixFunctions.jl”)
end

Then MatrixFunctions.jl can safely override Base.matrixfunction(::typeof(sqrt), A)

[robzan8]

I have implemented a couple of algorithms for computing matrix functions.
For dense matrix functions: Schur-Parlett with automatic differentiation done with TaylorSeries.jl.
For sparse matrix functions: the rational Krylov/Arnoldi method, with poles found automatically using AAA (the AAA algorithm for rational approximation of functions).

You can check them out here https://github.com/robzan8/MatFun.jl