Exponential

[sanderdemeyer]

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 Report

❌ Patch coverage is 0% with 109 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/implementations/exponential.jl	0.00%	87 Missing ⚠️
src/interface/exponential.jl	0.00%	9 Missing ⚠️
src/interface/matrixfunctions.jl	0.00%	8 Missing ⚠️
ext/MatrixAlgebraKitGenericSchurExt.jl	0.00%	3 Missing ⚠️
src/common/view.jl	0.00%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/MatrixAlgebraKit.jl	100.00% <ø> (ø)
src/common/view.jl	7.40% <0.00%> (-92.60%)	⬇️
ext/MatrixAlgebraKitGenericSchurExt.jl	0.00% <0.00%> (-100.00%)	⬇️
src/interface/matrixfunctions.jl	0.00% <0.00%> (ø)
src/interface/exponential.jl	0.00% <0.00%> (ø)
src/implementations/exponential.jl	0.00% <0.00%> (ø)

... and 36 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Jutho]

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]

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.

[lkdvos]

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

[sanderdemeyer]

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.

[lkdvos]

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.

[lkdvos]

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

[sanderdemeyer]

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]

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]

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]

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]

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]

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]

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.

[github-actions]

Your PR no longer requires formatting changes. Thank you for your contribution!

[sanderdemeyer]

The last few commits added the functions exponentiali and exponentiali!, for which it holds that exponentiali(t, A) = exponential(im*t*A). This can be useful for hermitian matrices A, where you still want to use the hermitian eigendecomposition for the matrix A, even though the matrix im*t*A is not per se hermitian.

[lkdvos]

Could we revisit the names for a second?
I think exponentiali looks a bit off, Julia tends to use im, and the combination of li! at the end all being characters that look very similar make it a bit hard to read too.
Just throwing some thoughts out there:

• Can we get away with a single function, and distinguish between the two cases using dispatch?
• Base uses cis (cos + i sin), is there a way to use something similar?
• It seems strange to support a scalar only in the imaginary case, and not for a real case, I would do either both or neither.
• There are still a number of suggestions and unaddressed comments, unless the github view is misleading me?

[sanderdemeyer]

• For me, exponentialim! would also be fine, although I really don't have any strong preference (you can assume this will be the case for all future naming discussions). If you have suggestions and/or a strong preference, feel free to either comment it here or just push it directly.
• The code works perfectly fine for both real and complex scalars. The reason we multiply with i explicitly, is that for real scalars, there would be no benefit in using the method with two arguments w.r.t. just using exponential(t*A). The only benefit stemming from the second method is the ability to use the hermitian eigendecomposition for the exponential of skew-hermitian matrices. The reason why two different names are used, is to make it fit with the current framework of @functiondef, where we need an @functiondef for exponential and an @functiondef n_args = 2 for exponentali, which do not work together. The current framework just doens't fit nicely with defining such constructions, and rather than throwing it all away and implementing it all ourselves for this function (which would not fit nicely with the rest of this package), defining the functions with 2 different names (1) allows us to use the current structure of this package, and (2) fits nicely with the argument of adding the imaginary unit. I would be okay with just not defining exponential(t, A), but @Jutho argued that this could be useful for real-time evolution (e.g. with cluster expansions), which is why I agree with him that this is the best approach.

There are a few comments that occur twice in the view, that I have reacted to once. The comments that I have taken into account and have reacted a 'thumbs up' to, I haven't resolved to make sure you saw they were taken into account. There are, I think, two comments left:

• the one where you said "Let me rephrase: I mean to first add that and then rebase this on top.". I don't really get what you mean with this comment. If you mean to first merge another PR (that hasn't been opened yet), I would really prefer to switch the order in view of time, since I would like to have the exponential function in the next version of TensorKit.
• the one where you said "This implementation is a bit strange, LinearAlgebra actually uses exp!(A) as its expert driver, so at least this is allowed to destroy the input matrix, but additionally you should also try and avoid allocating expA, and simply have expA === A". Github marks this comment as outdated, and I don't see which lines you commented this on (this is also already 3 weeks old). So could you either clarify what you mean or, like I said before, just push the desired changes directly if you have strong opinions on them.

[Jutho]

I think I finally managed to work away (some of) my backlog of PRs that I had to review, so this one is among the next on my list. I hope to go through it (and if necessary commit changes myself directly) before the end of the week.

[lkdvos]

I cleaned up part of the implementation, with two major changes:

• expA is only allocated if a scalartype change is required
• I refactored/added map_diagonal(!) in an attempt to have a specialization hook, thinking ahead a bit about extending this to different types and different matrix functions

I feel like with this change, the specific difference for exponentiali feels even smaller.
I would still argue that it is a bit strange to default to an imaginary scalar while not supporting a real one, and instead would argue that it could be more convenient to have a generic scalar for which you could simply fill in an imaginary value.
Thinking about that a bit more, I even would say that the reason we need to have the imaginary_evolution flag in MPSKit for time evolution is precisely because we defaulted to evolving with exp(im * t * H) instead of exp(t * H).
It is easy to detect real vs complex scalars in a type stable manner, but less so to detect purely imaginary vs complex.
(I hope I'm explaining this somewhat decently)

Trying to come up with some solutions I think could be better, I can think of some alternative approaches:

• We could introduce a ScaledOperator(A, tau) dedicated type, and then (re)use exponential(ScaledOperator(A, tau), ...). This might remove a little bit of overhead, and a default implementation could just instantiate the scaled operator. This avoids the need for @functiondef 2 ...
• We could simply always have a scale factor in the exponential function, which by default takes the value true or VectorInterface.One(), and just correctly handle the complex/real case by checking if this scalar type is real/complex.

---

Some more discussion points/things to consider, some I recall briefly bringing up but might be good to write down here:

These matrix functions should have the same output types, no matter whether they are ishermitian or not.
Instead, most matrix functions distinguish between having support for the entire \bbR -> \bbR, or require complex extensions to do so (this is of course not relevant for exp).
Therefore, it might actually be reasonable to instead of splitting MatrixFunctionViaEig and MatrixFunctionViaEigh into separate types, simply having a MatrixFunctionViaEigen with a runtime flag for hermitianness, rather than a compile time flag. (in essence, manually union-splitting with an if statement) This could potentially cut down on compile-times, and the number of required algorithms?

It might be reasonable to have a global LinearAlgebraAlgorithm backend, which functions similar to DiagonalAlgorithm in the sense that it simply dispatches to the respective LinearAlgebra routines.
With this, I think we could deprecate most of the GenericLinearAlgebra extension, as this would simply coincide.

[lkdvos]

TODO: do we want to keep this check here, knowing that it will again be checked in the iimplementation of eigh_full?

[sanderdemeyer]

Indeed, we should remove this test.

[lkdvos]

Same comment here

[lkdvos]

This stood out to me as a very specific choice, but it looks like it is hard to make this converge. Do we know what is going on here?

[sanderdemeyer]

The m = 2 choice was only made to make the tests run more quickly while running them locally, I must have forgotten to change it back. If I change it to m = 54 (like in the other tests), all tests still pass in my case (that is, in my latest version, without the commits of yesterday.)

[sanderdemeyer]

I don't get why you say the current code would not work with a real scalar. The argument τ is restricted to be a Number, not per se a real or imaginary one. Hence, it just corresponds to an irrelevant convention choice. If you really want to use it with a real argument (which would be useless, since you can just use the exponential function), you can still call exponentiali(-iτ, A). This would work perfectly fine. Indeed, this feels a bit unintuitive, so that's why we can just use the exponential function. I also think that the current implementation allows for both imaginary_evolution and real_evolution.

Apart from the fact that I don't see why the changes would make the difference between exponential! and exponentiali smaller (the fundamental problem remains), I feel like the proposed changes make everything more difficult than they need to be. You have always been an advocate for making an interface feel intuitive. I think that exp(A) = exponential(A) and exp(itA) = exponentiali(t, A)feel more intuitive than eitherexponential(ScaledOperator(it, A))orexp(A) = exponential(1, A) = ...`.

In terms of differentiating between eigh and eig during runtime: Then why do we not do this for eigh and eig themselves?

If you really feel strongly about any of this, feel free to tailor the implementation to how you see things. The only thing that I really care about is an implementation that works, preferably as generic as possible (again, because I would like to have this in the next TensorKit version). As long as the tests pass, this is all fine by me. I might not be a big fan of the choice of implementation, but naming conventions and interface choices are way less important to me than having either implementation.

[lkdvos]

I don't get why you say the current code would not work with a real scalar.

The difference is a bit subtle, but as I mentioned before this is about type stability.
(It is easy to detect real vs complex scalars in a type stable manner, but less so to detect purely imaginary vs complex.)
Maybe it is more obvious with an example:

function exponentialr(tau, A)
    expA = tau isa Complex ? complex(A) : A
    return exponentialr!(tau, A, expA) # this computes expA .= exp(A * tau) in whatever way
end
function exponentiali(t, A)
    expA = iszero(real(t)) ? A : complex(A)
    return exponentiali!(t, A, expA) # this computes expA .= exp(A * t * im) in whatever way
end

This first implementation is type stable for both purely real, purely imaginary and complex values of tau, while the second isn't.

---

In terms of differentiating between eigh and eig during runtime: Then why do we not do this for eigh and eig themselves?

The reason we opted for eigh and eig is that, unlike LinearAlgebra.eigen, this can be implemented in a type-stable manner for real input matrices. eigh will return real D and V, while eig will always return complex D and V.

For exponential we can revisit this, since both hermitian and non-hermitian real inputs get real outputs.

---

(again, because I would like to have this in the next TensorKit version)

I’m afraid that might be difficult: we’re trying to get a version released this week, so I don’t think integrating and polishing this in time is very realistic.

That said, if your main goal is to get something working for your use case in the next TensorKit release, you could already move ahead with something along the following lines:

function myexp(A::AbstractTensorMap)
    D, V = eig_full(A)
    D.data .= exp.(D.data)
    return V * D * inv(V)
end

From my side, I think it’s important that we take a bit of time here to come up with a good, consistent interface and to think through allocations and long-term maintenance, especially since this will likely influence all of the matrix functions. I’d really like to avoid having to introduce multiple breaking changes later.

But that shouldn’t block you from experimenting or using a simpler version for your immediate needs.