Skip to content

Exponential#94

Open
sanderdemeyer wants to merge 67 commits into
QuantumKitHub:mainfrom
sanderdemeyer:exponential
Open

Exponential#94
sanderdemeyer wants to merge 67 commits into
QuantumKitHub:mainfrom
sanderdemeyer:exponential

Conversation

@sanderdemeyer

@sanderdemeyer sanderdemeyer commented Nov 12, 2025

Copy link
Copy Markdown
Contributor

This implements the exponential of a matrix for both BLASFloats and BigFloats.

I have named these functions exponential and exponential!, instead of the usual exp and exp! from LinearAlgebra. Extending these methods while keeping the current structure using @algdef and @ functiondef results in some naming conflicts. The default for BLASFloats is to use LinearAlgebra.exp!. In TensorKit, we can still stick to the exp naming convention.

@codecov

codecov Bot commented Nov 12, 2025

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 86.17021% with 13 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
src/interface/exponential.jl 0.00% 9 Missing ⚠️
ext/MatrixAlgebraKitGenericSchurExt.jl 0.00% 3 Missing ⚠️
src/common/view.jl 50.00% 1 Missing ⚠️
Files with missing lines Coverage Δ
src/MatrixAlgebraKit.jl 100.00% <ø> (ø)
src/implementations/exponential.jl 100.00% <100.00%> (ø)
src/interface/matrixfunctions.jl 100.00% <100.00%> (ø)
src/common/view.jl 96.55% <50.00%> (-3.45%) ⬇️
ext/MatrixAlgebraKitGenericSchurExt.jl 81.25% <0.00%> (-18.75%) ⬇️
src/interface/exponential.jl 0.00% <0.00%> (ø)
🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/interface/exponential.jl Outdated
Comment thread test/exponential.jl Outdated
Comment thread test/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated

function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::ExponentialViaEigh)
D, V = eigh_full(A, alg.eigh_alg)
copyto!(expA, V * Diagonal(exp.(diagview(D))) * inv(V))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reduced allocation strategy:

Suggested change
copyto!(expA, V * Diagonal(exp.(diagview(D))) * inv(V))
iV = inv(V)
map!(exp, diagview(D))
mul!(expA, rmul!(V, D), iV)

@sanderdemeyer sanderdemeyer Nov 13, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It has to be map!(exp, diagview(D), diagview(D)) instead of map!(exp, diagview(D)), but good suggestion otherwise. I have also added it for the ExponentialViaEig.
EDIT: the suggested change works only for Julia 1.12 onwards. That's why I will keep the version with
3 arguments.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not just diagview(D) .= exp.(diagview(D))?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is that more efficient than the current code? If not, I'd prefer to keep it that way, since it feels a bit more natural to me.

Comment thread src/implementations/exponential.jl Outdated
Comment thread src/interface/decompositions.jl Outdated
change some input tests
remove redundant comment
include ComplexF16 in tests
fix unchanged test names and docs
improve allocations

@lkdvos lkdvos left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Overall I'm not fully convinced by the interface of exponential(!), especially in its current form and implementation this looks slightly strange.

LinearAlgebra uses an in-place version, i.e. it reuses the input array to return exp!, and looking at the different implementations you have here, it is not obvious that trying to fit this into a exponentiate!(A, expA, alg) signature is really helping us - on the contrary, all this is really doing is creating an additional copy at the end just to make sure that it is allocated in the provided output.
As we discussed for your previous PR, this really is not the purpose of being able to provide the output argument.

For the algorithms, thinking a bit ahead, it might be appropriate to just call these something along the lines of matrix functions via eig, since presumably these approaches are actually generic for all of these implementations.

Comment thread src/interface/matrixfunctions.jl
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated

function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::ExponentialViaEigh)
D, V = eigh_full(A, alg.eigh_alg)
copyto!(expA, V * Diagonal(exp.(diagview(D))) * inv(V))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not just diagview(D) .= exp.(diagview(D))?

Comment thread ext/MatrixAlgebraKitGenericSchurExt.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
@sanderdemeyer

sanderdemeyer commented Nov 13, 2025

Copy link
Copy Markdown
Contributor Author

Overall I'm not fully convinced by the interface of exponential(!), especially in its current form and implementation this looks slightly strange.

LinearAlgebra uses an in-place version, i.e. it reuses the input array to return exp!, and looking at the different implementations you have here, it is not obvious that trying to fit this into a exponentiate!(A, expA, alg) signature is really helping us - on the contrary, all this is really doing is creating an additional copy at the end just to make sure that it is allocated in the provided output. As we discussed for your previous PR, this really is not the purpose of being able to provide the output argument.

The idea of putting it in this framework is to allow for different sorts of algorithms, i.e. ones that would also work for BigFloats. I get that in the BLASFloat case, we should avoid extra allocations, but could you elaborate on your suggestion for ExponentialViaLA? Do you just want to get rid of the preallocated output?

For the algorithms, thinking a bit ahead, it might be appropriate to just call these something along the lines of matrix functions via eig, since presumably these approaches are actually generic for all of these implementations.

I may be missing your point here, but that was the idea behind the ExponentialViaEigh and ExponentialViaEig naming conventions? Or do you really want separate names for each eig_alg?

@lkdvos

lkdvos commented Nov 13, 2025

Copy link
Copy Markdown
Member

The idea of putting it in this framework is to allow for different sorts of algorithms, i.e. ones that would also work for BigFloats. I get that in the BLASFloat case, we should avoid extra allocations, but could you elaborate on your suggestion for ExponentialViaLA? Do you just want to get rid of the preallocated output?

I am indeed referring to the preallocated output, not the algorithms part. It really only makes sense to have the option of giving a preallocated output if we are actually able to use this, and for the current implementations you have this is not saving us any work, rather it is increasing it because you add an extra allocation at the beginning and an extra copy at the end.

While it is definitely possible to have initialize_output(exponentiate!, A) = A to just allocate the result in-place, I still wonder if it is then useful to have to implement the boiler plate for fitting it in a framework that we are mostly going to bypass anyways.

I may be missing your point here, but that was the idea behind the ExponentialViaEigh and ExponentialViaEig naming conventions? Or do you really want separate names for each eig_alg?

Sorry I should have explained that better, I meant that I would like to avoid having to also define SqrtViaEig, SinViaEig, ..., and simply have some algorithm that signifies "solve this by doing an eigenvalue decomposition".

@sanderdemeyer

sanderdemeyer commented Nov 13, 2025

Copy link
Copy Markdown
Contributor Author

I am indeed referring to the preallocated output, not the algorithms part. It really only makes sense to have the option of giving a preallocated output if we are actually able to use this, and for the current implementations you have this is not saving us any work, rather it is increasing it because you add an extra allocation at the beginning and an extra copy at the end.

While it is definitely possible to have initialize_output(exponentiate!, A) = A to just allocate the result in-place, I still wonder if it is then useful to have to implement the boiler plate for fitting it in a framework that we are mostly going to bypass anyways.

Is your suggestion then to just skip the whole @ functiondef and other general frameworks we have to define exponentiate and exponentiate!, or to keep the current framework somewhat, but just remove the preallocated output argument?

Sorry I should have explained that better, I meant that I would like to avoid having to also define SqrtViaEig, SinViaEig, ..., and simply have some algorithm that signifies "solve this by doing an eigenvalue decomposition".

Okay, I see. I agree and will change this.

@Jutho

Jutho commented Nov 13, 2025

Copy link
Copy Markdown
Member

Regarding the algorithm names, I agree with Lukas and also think we want to have a general
MatrixFunctionViaEig, though it is useful to know about which function exactly, because some functions map real matrices to real matrices (even if complex eigenvalues are involved), and others do not. I think these are the two major categories, so we could just have two different algorithms for that, but I don't know good names for that.

Regarding the role of the output arguments, I only partially agree. The whole point of why we started MatrixAlgebraKit.jl, is because in TensorKit we first want to define the output tensor, and then compute block per block the result, where we want to store the result in the corresponding block of the output tensor. Ideally, yes, the computation is such that we also use that output data as storage during the computation, in such a way that the end result "naturally" ends up there, but if that is difficult, a final copy! can still be useful. We can then always try to improve this later on behind the scenes, but at least TensorKit can be agnostic about this.

Note that the LinearAlgebra exp! is also cheating, and typically also ends up just copying the final result back into A. There is no way to actually compute exp(A) fully in-place. The typical approach via Pade approximations has a ton of allocations. I think I actually have a slightly allocation-friendlier implementation lying around for some things we did in CMPSKit.jl

@Jutho

Jutho commented Nov 13, 2025

Copy link
Copy Markdown
Member

Regardless of the comment about general matrix functions, it is a fact that the exponential is by far the most useful and common one that we need, so I am also not opposed to first thinking carefully about this one, and having some part of the implementation be specific for matrix exponentials.

In particular, one important consideration that we might want to include in this design, that is specific to our use case, is that we also might be interested in computing exp(-im* δt*H) for real time evolution, where it is probably useful to use eigh for the Hermitian matrix H, but still the end result will be complex. Julia has a cis function for cis(x) = exp(im * x), but unfortunately that comes without a minus sign in the argument. But maybe cis(-δt * H) is not too bad.

@lkdvos

lkdvos commented Nov 13, 2025

Copy link
Copy Markdown
Member

To comment on the TensorKit interaction, I definitely agree with the purpose, but this is not actually currently the design we ended up with.
Our interface very clearly states that we may use the provided output, but are in no way bound to do so.
This is reflected in the TensorKit implementation, see e.g. https://github.com/QuantumKitHub/TensorKit.jl/blob/feab2072b2909b9a127d871ef58c36e83623c2fc/src/factorizations/matrixalgebrakit.jl#L39
I'm of course happy to rediscuss that design decision, but I'm actually quite happy with how that is shaping out.

So basically there are two comments I have:

On the one hand, there is the question about whether or not there are implementations that benefit from providing an additional output array.
From looking at the implementations here, this is not obvious at all, actually it seems like instead we are simply allocating an additional array at the end to copy the results into, which we might as well have reused the A for, since we already "destroyed" that data anyways.
I am not saying that the implementations have to be fully allocation-free (the LAPACK routines still have workspaces too), I just think that if providing the output is a way to avoid having to do an allocation, it really shouldn't be increasing the allocations instead.
If we can't come up with a case where having this additional array actually helps us out, it might be reasonable to simplify define exponentiate!(A, alg) -> expA where expA === A instead, which would still work for TensorKit in the exact same way.
I definitely agree that this can be achieved by hacking it into the current interface as well, by simply having initialize_output(exponentiate!, A) = A, I just wanted to open up the discussion about whether or not that is something we like, or think this is just boilerplate we aren't using.
Giving it some thought, I'm definitely okay with the initialize_output(...) = A default and going from there, i.e. keeping @functiondef and "hacking" this into that.

On the other hand, given that interface, if there is no way of naturally making the output end up in the provided destination, I would really like to avoid ending up with a final copy!(provided_dest, computed_exp) at the end, even as a "for now" implementation that we can improve on later.
Given that we have to check for in-placeness in TensorKit anyways, I would in this case much rather have initialize_output() = nothing, and simply return the computed_exp instead.
This is the exact same comment and solution as for the generic linear algebra wrappers.

@Jutho

Jutho commented Nov 13, 2025

Copy link
Copy Markdown
Member

I guess I am a bit confused, because most of the implementations now do actually perform the final step in the calculation in such a way that the result is directly stored in the output array, no? It is only the algorithm that goes via Base/LinearAlgebra that requires the extra copy! step, and that implementation I will happily replace with my own Pade implementation.

But it is also true that, by the time the final step of the calculation is reached; the memory of A, which has already been destroyed and is no longer containing active information, can be reclaimed to store the result, so initialize_output(...) = A could indeed be the right design choice.

sanderdemeyer and others added 4 commits November 19, 2025 15:07
change name to `MatrixFunctionViaEig` etc
change `decompositions` to `matrixfunctions`
add default algorithm for Diagonal matrices
add input checks
add @testthrows to catch non-hermitian matrices being given to MatrixFunctionViaEigh
change default exponential algorithm to e.g. `MatrixFunctionViaEig` of the default `eig_alg`
Comment thread src/implementations/exponential.jl Outdated
Comment on lines +65 to +71
function exponential!(A, expA, alg::MatrixFunctionViaEigh)
check_input(exponential!, A, expA, alg)
D, V = eigh_full!(A, alg.eigh_alg)
expD = map_diagonal!(x -> exp(x / 2), D, D)
VexpD = rmul!(V, expD)
return mul!(expA, VexpD, V')
end

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function exponential!(A, expA, alg::MatrixFunctionViaEigh)
check_input(exponential!, A, expA, alg)
D, V = eigh_full!(A, alg.eigh_alg)
expD = map_diagonal!(x -> exp(x / 2), D, D)
VexpD = rmul!(V, expD)
return mul!(expA, VexpD, V')
end
exponential!(A, expA, alg::MatrixFunctionViaEigh) = exponential!((1, A), expA, alg)

Comment thread src/implementations/exponential.jl Outdated
Comment on lines +108 to +116
expD = map_diagonal!(x -> exp(x * τ), D, D)
iV = inv(V)
VexpD = rmul!(V, expD)
if eltype(A) <: Real && eltype(τ) <: Real
expA .= real.(VexpD * iV)
return expA
else
return mul!(expA, VexpD, iV)
end

@Jutho Jutho Jun 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
expD = map_diagonal!(x -> exp(x * τ), D, D)
iV = inv(V)
VexpD = rmul!(V, expD)
if eltype(A) <: Real && eltype(τ) <: Real
expA .= real.(VexpD * iV)
return expA
else
return mul!(expA, VexpD, iV)
end
if eltype(A) <: Real && eltype(τ) <: Real
VexpD = V .* exp.(transpose(diagview(D)) .* τ)
expAc = rdiv!(VexpD, LinearAlgebra.lu!(V))
return expA .= real.(expAc)
else
expA .= V .* exp.(transpose(D) .* τ)
return rdiv!(expA, LinearAlgebra.lu!(V))
end

Comment thread src/implementations/exponential.jl
Comment thread src/implementations/exponential.jl Outdated

function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::DiagonalAlgorithm)
check_input(exponential!, (τ, A), expA, alg)
return map_diagonal!(x -> exp(x * τ), expA, A)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return map_diagonal!(x -> exp(x * τ), expA, A)
diagview(expA) .= exp.(diagview(A) .* τ)
return expA

Comment thread src/interface/exponential.jl Outdated
Comment thread src/interface/matrixfunctions.jl Outdated
Comment thread src/interface/matrixfunctions.jl Outdated
Comment thread src/interface/matrixfunctions.jl Outdated

@Jutho Jutho left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there some more required changes that could improve allocations, code duplication and in one case correctness.

Comment thread src/implementations/exponential.jl Outdated
sanderdemeyer and others added 6 commits June 22, 2026 11:54
Co-authored-by: Jutho <Jutho@users.noreply.github.com>
Co-authored-by: Jutho <Jutho@users.noreply.github.com>
Co-authored-by: Jutho <Jutho@users.noreply.github.com>
Co-authored-by: Jutho <Jutho@users.noreply.github.com>
@sanderdemeyer

Copy link
Copy Markdown
Contributor Author

I think all comments are now taken care of (except for the discussion of whether to use 1 or one(eltype(A)), where I've chosen the latter.)

Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated
Comment on lines +73 to +74
VexpD = V .* exp.(transpose(diagview(D)) .* τ)
return mul!(expA, VexpD, V')

@lkdvos lkdvos Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was looking at this and slightly breaking my brain, if everything is complex than D is real (because hermitian), but we can still get away with rmul!(V, exponential!((τ / 2, D), D)) since V should be complex, and then we basically recover the path from above. Only in the case where A is real, tau is complex do we run into trouble because the in-place rmul! has to become a new allocation due to the element type changing.
However, even in that case it's probably better to use the tau/2 way because mul!(expA, VexpD, VexpD') can then use BLAS, since all matrices involved have complex entries. (I think in principle we can do mul!(complexC, complexA, realB), but not sure if that is automatically implemented?)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If tau is complex but A real and symmetric, we shouldn't be doing VexpD * VexpD' cause that gives something different, but VexpD * transpose(VexpD) might work.

It was to avoid these different eltype combinations that I used broadcasting, but I think the point is now to avoid broadcasting altogether, right?

@lkdvos lkdvos Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think by avoiding broadcasting it would work with TensorMap instances out of the box, but this might also really be taking it too far and it might not matter all that much for a first implementation, and we can always revisit to optimize

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the consensus then to keep the broadcasting here? There are a few other instances where broadcasting is used, and thus in order for exponential(!) to work with TensorMaps, there will anyway be some work needed in TensorKit.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But now it is completely inconsistent no? In one branch of an if else we do call exponential!(::Diagonal) and in the other we do broadcasting?

I think we either want the ::ExponentialViaEig(h) to either be fully written using exponential!(::Diagonal) and generic operations such as mul! and rmul!, so that it is directly useful for other types.

If we do not succeed, and anyway use broadcasting in part of it, I preferred my older version that just went all in with broadcasting.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed it such that it avoids broadcasting altogether. To avoid allocations as much as possible, I just consider the 4 different cases now.

Comment thread src/implementations/exponential.jl Outdated
Comment thread src/implementations/exponential.jl Outdated

@lkdvos lkdvos left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorry for being a bit slow on this, but I just realized that probably exponential!((tau, D)) is enough instead of the exponential!((tau, D), D, DiagonalAlgorithm()), which should make a lot of the code more generic and might allow us to just have it work for the TensorMap's instead.

Comment on lines +108 to +109
diagview(expA) .= exp.(diagview(A) .* τ)
return expA

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
diagview(expA) .= exp.(diagview(A) .* τ)
return expA
return map_diagonal!(x -> exp(x * τ), expA, A)

Comment on lines +9 to +12
exponential((τ,A); kwargs...) -> expτA
exponential((τ,A), alg::AbstractAlgorithm) -> expτA
exponential!((τ,A), [expA]; kwargs...) -> expτA
exponential!((τ,A), [expA], alg::AbstractAlgorithm) -> expτA

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
exponential((τ,A); kwargs...) -> expτA
exponential((τ,A), alg::AbstractAlgorithm) -> expτA
exponential!((τ,A), [expA]; kwargs...) -> expτA
exponential!((τ,A), [expA], alg::AbstractAlgorithm) -> expτA
exponential((τ, A); kwargs...) -> expτA
exponential((τ, A), alg::AbstractAlgorithm) -> expτA
exponential!((τ, A), [expA]; kwargs...) -> expτA
exponential!((τ, A), [expA], alg::AbstractAlgorithm) -> expτA

exponential!((τ,A), [expA]; kwargs...) -> expτA
exponential!((τ,A), [expA], alg::AbstractAlgorithm) -> expτA

Compute the exponential of the square matrix `A` or `τ*A`,

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Compute the exponential of the square matrix `A` or `τ*A`,
Compute the exponential of the square matrix `A` or `τ * A`,

return map_diagonal!(exp, expA, A)
end

function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::DiagonalAlgorithm)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA::AbstractMatrix, alg::DiagonalAlgorithm)
function exponential!((τ, A)::Tuple{Number, Any}, expA, alg::DiagonalAlgorithm)


# Diagonal logic
# --------------
function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::DiagonalAlgorithm)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::DiagonalAlgorithm)
function exponential!(A, expA, alg::DiagonalAlgorithm)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants