Skip to content

add det(Symmetric) and method= to expm(Matrix)#13093

Closed
alyst wants to merge 3 commits intoJuliaLang:masterfrom
alyst:symm_fixes
Closed

add det(Symmetric) and method= to expm(Matrix)#13093
alyst wants to merge 3 commits intoJuliaLang:masterfrom
alyst:symm_fixes

Conversation

@alyst
Copy link
Copy Markdown
Contributor

@alyst alyst commented Sep 12, 2015

  • There was no det(Symmetric) (see det(Symmetric) fails LinearAlgebra.jl#252). Now it just falls back to the dense case, don't know if there's better solution.
  • Added keyword argument method=:auto/:higham to expm(AbstractMatrix) to avoid switching to Symmetric branch if the dense matrix is symmetric. Otherwise, tests for exponential of symmetric matrix are essentially comparing expm(Symmetric) to itself. Also, given many numerical stability problems of expm(), it would be convenient to have several numerical methods to choose from.

det(expm(mtx)) should be 1.0, this should be a good test for numerical
stability
@kshyatt kshyatt added the linear algebra Linear algebra label Sep 12, 2015
@stevengj
Copy link
Copy Markdown
Member

If the matrix is symmetric, isn't the Symmetric branch the best algorithm? There shouldn't be any stability problems with expm for symmetric matrices, where the eigenvectors are orthogonal. I don't think expm(AbstractMatrix) should default to a suboptimal algorithm just to make testing easier.

@stevengj
Copy link
Copy Markdown
Member

I guess for Symmetric we could transform to a SymTridiagonal matrix using the standard similarity transform and then use ldltfact to get the determinant, rather than LU factorizing the dense matrix. But maybe that is slower.

@alyst
Copy link
Copy Markdown
Contributor Author

alyst commented Sep 12, 2015

For example, if

A = [-0.8599805912689542 44.28830587812244
     44.28830587812245    0.8599805912689469]
julia> eigvals(A)
JuliaLang/julia#2-element Array{Float64,1}:
# -44.2967
#  44.2967

Then trace(A) is exactly zero and so det(expm(A)) == 1, but because of a big difference between A eigenvalues different numerical methods can give exponentials with either 0 or very big determinants.

@alyst
Copy link
Copy Markdown
Contributor Author

alyst commented Sep 12, 2015

The alternative to expm!(method=) would be, of course, extracting expm_nonsymm!() method from expm!().

It would also be possible to eliminate two matrix allocations when expm(M) is called for symmetric dense matrix:

  • one if check for symmetry is done before M is copied
  • another one if full(expm(Symmetric(M))) is replaced by full!(expm(Symmetric(M))). Destructive full!() doesn't exist, but it could be defined as copytri!(A.data, A.uplo)

@Keno
Copy link
Copy Markdown
Member

Keno commented Feb 28, 2020

I think it's fair to say that we don't like these kinds cases where keyword arguments are used to disable a specific piece of functionality. If this is still of interest, I would suggest extracting exp_higham! from the exp! method. I don't think it needs to be propagated to the non-destructive versions, because the use case for this is quite advanced, so such users can use the destructive versions directly. (I also don't think exp_higham! should be exported).

@Keno Keno closed this Feb 28, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

linear algebra Linear algebra

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants