WIP: Improve accuracy of ^ for ill-conditioned matrices

[mfasi]

Add the algorithm for A^t from Higham and Lin, 2013. See #5840.

I split the old version of logm into several functions to reuse some of the code.

[StefanKarpinski]

"Improve" is a little vague. Does this change behavior or just performance? Some explanation and examples of the difference would be nice.

[mfasi]

@StefanKarpinski You're right, an explanation might help. The implementation that we have now in base is not reliable for general matrices because it fails when the eigenvectors of the matrix are ill-conditioned. A classical example is the comparison of the (relative) residuals for

julia> e= 1e-16; A = [1+e 1; 0 1];
julia> cond(A)
2.618033988749896
julia> F = eigfact(A); cond(F.vectors)
9.007199254740992e15

On master I get

julia> norm((A^(1/2))^2 - A) / norm(A)
0.6180339887498948

and on powm

julia> norm((A^(1/2))^2 - A) / norm(A)
6.861555643110582e-17

To answer your question, I'd say that the focus is on the behaviour here.

[stevengj]

Are we worried about type-stability of this function? Or is it slow enough that we don't care?

[StefanKarpinski]

If we can be type stable without sacrificing too much convenience, it's better to, even if things aren't performance critical.

[andreasnoack]

So far we haven't been type stable for matrix functions like this, e.g. eig and sqrtm return Float64 for Float64 input if possible. If we want type stability, we could

1. Mimic scalar function and return DomainError
2. Always return Complex

I think option 1 is a no go because the consequence would be to solve real eigenvalue problems in complex arithmetic which is much more expensive. The main problem with option two is that some people might be annoyed by the complex return matrix when all entries are real valued, but that might be a reasonable price for type stability.

[mfasi]

@stevengj Thank you for checking the patch thoroughly. I addressed most of the points you raised, still to do:

• add tests
• add documentation
• push cond(A) after the Schur decomposition

I'll make this PR a WIP and will probably add one or two other commits (tests and docs). Not sure I found the least expensive way to do the in place operations, so feel free to check the loops again (I'll mark them with a line comment for convenience).

[mfasi]

I managed to push the use of cond(A) after the Schur factorization. I added inline documentation and tests as well.

[stevengj]

I think it would probably be clearer, and consistent with standard practice in linear algebra, to name your upper-triangular matrices U instead of A, here and below.

[mfasi]

Agreed, but I'd rather prefer not to do it with this PR, because it would break the consistency with the notation of sqrtm. I'll put it in my TODO list.

[stevengj]

Uh oh, looks like something went terribly wrong with a rebase...

[mfasi]

@stevengj Yes, you're right. No Idea how to fix this (I'm trying but it looks like I went a git push -f too far), I think I will have to start over with a new branch

[alyst]

@mfasi To avoid making a new PR you can start a new local branch from the current master, cherry-pick the relevant commits from this branch, make sure that it works for you locally, then replace this github branch with your new local one.

[andreasnoack]

@mfasi Did you ever manage to recover this work on a new branch? As mentioned, we'd be very interested in getting this into Julia.

[stevengj]

Merged in #21184.