Consistency of return type for matrix functions

[mfasi]

I added specialized versions of expm, logm and sqrtm for symmetric and Hermitian matrices.

EDIT: The documentation is now in #12432

[KristofferC]

This fails the whitespace checker.

[mfasi]

I made the templates clearer, as suggested by @jiahao, and added a bunch of tests to cover a few missing lines in dense.jl.

[andreasnoack]

I think the convertions with full are right when the input is a StridedMatrix. That is also how we do it in inv.

Still if a matrix is either Hermitian or Symmetric, then f(A) will be Symmetric for sure

I think this is only correct for real elements. E.g.

julia> A = [1 im; -im 1.0]      
2x2 Array{Complex{Float64},2}:  
 1.0+0.0im  0.0+1.0im           
 0.0-1.0im  1.0+0.0im         

julia> ishermitian(A)        
true                         

julia> expm(A) |> t -> norm(t - t.')      
6.38905609893065                            

I'm not an expert on matrix functions, but we probably want for keep separate defition loops for entire and non-entire functions. I guess we can ensure type stability for entire functions with Symmetric and Hermitian arguments whereas this is not the the case for the logarithm and root functions.

I forgot to mention that it's a great update of the documentation.

[andreasnoack]

possible a nonprimary matrix function is returned

Should it be "non-principal"?

[mfasi]

@andreasnoack You're right, thank you for spotting this, I'll fix it and submit a separate PR just for the documentation (I'll have to adapt it to #11943).

[mfasi]

The documentation in now in #12432.

[mfasi]

@andreasnoack You're right, that condition is true only for real Hermitian, that is Symmetric. On the other hand, I think that you'd be right in real arithmetic, but numerically you cannot ensure that expm(A) is Hermitian by just knowing that A is Hermitian. An example that shows this (unless you're very lucky):

A1 = randn(4,4) + im*randn(4,4); A2 = A1' + A1; ishermitian(A2)

If you check diag(expm(A2)) you'll see that the imaginary part of the elements on the diagonal is small but not 0, so numerically the matrix is not Hermitian.

EDIT: it looks like this can happen even if A2 is positive definite, for example

A1 = randn(4,4) + im*randn(4,4); A2 = A1' + A1; ishermitian(A2)

which quite settles the question, I think. If A is Symmetric, so is f(A). If it is Hermitian, the result can be of any type, so returning a StridedMatrix seems appropriate.

[andreasnoack]

@mfasi Part of the point in having the Hermitian and Symmetric types is circumvent this numerical noice problem and exploit what we know mathematically, i.e. that expm(Hermitian) is Hermitian. We are doing an exact check of the imaginary part of the diagonal on construction which would cause an error in the example you present, but I thinks we should then simply set the imaginary parts to zero.

[mfasi]

@andreasnoack I submitted a version that enforces the return type when the perturbation of the Hermitian structure is due to rounding errors. As the exponential is the only entire function we have, I didn't implement a loop for it.

[andreasnoack]

Why is real necessary here?

[andreasnoack]

@mfasi Is this ready? I don't think I have more comment. Would you do a git rebase origin/master and forced push then we might be lucky enough to get AppVeyor running the tests.

[mfasi]

@andreasnoack Yes, I thought it was ready. I rebased on the master branch, but now the apple build on Travis fails as well. Not sure what happened, since I cannot access the whole log file, but I'm pretty sure the problem is not related to what I changed, it looks like it's trying to compile some fortran library. I'll try to rebase and force again as soon as a new commit is merged to the master branch.

[mfasi]

Both Travis and AppVeyor managed to build and check, feel free to merge.

[andreasnoack]

Great. Thanks.