Exponential#94
Conversation
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
|
|
||
| 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)) |
There was a problem hiding this comment.
Reduced allocation strategy:
| copyto!(expA, V * Diagonal(exp.(diagview(D))) * inv(V)) | |
| iV = inv(V) | |
| map!(exp, diagview(D)) | |
| mul!(expA, rmul!(V, D), iV) |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Why not just diagview(D) .= exp.(diagview(D))?
There was a problem hiding this comment.
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
left a comment
There was a problem hiding this comment.
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.
|
|
||
| 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)) |
There was a problem hiding this comment.
Why not just diagview(D) .= exp.(diagview(D))?
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
I may be missing your point here, but that was the idea behind the |
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
Sorry I should have explained that better, I meant that I would like to avoid having to also define |
Is your suggestion then to just skip the whole
Okay, I see. I agree and will change this. |
|
Regarding the algorithm names, I agree with Lukas and also think we want to have a general 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 Note that the LinearAlgebra |
|
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 |
|
To comment on the TensorKit interaction, I definitely agree with the purpose, but this is not actually currently the design we ended up with. 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. 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 |
|
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 But it is also true that, by the time the final step of the calculation is reached; the memory of |
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`
| 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 |
There was a problem hiding this comment.
| 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) |
| 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 |
There was a problem hiding this comment.
| 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 |
|
|
||
| function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::DiagonalAlgorithm) | ||
| check_input(exponential!, (τ, A), expA, alg) | ||
| return map_diagonal!(x -> exp(x * τ), expA, A) |
There was a problem hiding this comment.
| return map_diagonal!(x -> exp(x * τ), expA, A) | |
| diagview(expA) .= exp.(diagview(A) .* τ) | |
| return expA |
Jutho
left a comment
There was a problem hiding this comment.
I think there some more required changes that could improve allocations, code duplication and in one case correctness.
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>
|
I think all comments are now taken care of (except for the discussion of whether to use |
| VexpD = V .* exp.(transpose(diagview(D)) .* τ) | ||
| return mul!(expA, VexpD, V') |
There was a problem hiding this comment.
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?)
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
I changed it such that it avoids broadcasting altogether. To avoid allocations as much as possible, I just consider the 4 different cases now.
lkdvos
left a comment
There was a problem hiding this comment.
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.
Co-authored-by: Lukas Devos <ldevos98@gmail.com>
Co-authored-by: Lukas Devos <ldevos98@gmail.com>
Co-authored-by: Lukas Devos <ldevos98@gmail.com>
| diagview(expA) .= exp.(diagview(A) .* τ) | ||
| return expA |
There was a problem hiding this comment.
| diagview(expA) .= exp.(diagview(A) .* τ) | |
| return expA | |
| return map_diagonal!(x -> exp(x * τ), expA, 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 |
There was a problem hiding this comment.
| 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`, |
There was a problem hiding this comment.
| 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) |
There was a problem hiding this comment.
| 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) |
There was a problem hiding this comment.
| function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::DiagonalAlgorithm) | |
| function exponential!(A, expA, alg::DiagonalAlgorithm) |
This implements the exponential of a matrix for both
BLASFloatsandBigFloats.I have named these functions
exponentialandexponential!, instead of the usualexpandexp!fromLinearAlgebra. Extending these methods while keeping the current structure using @algdef and @ functiondef results in some naming conflicts. The default for BLASFloats is to useLinearAlgebra.exp!. InTensorKit, we can still stick to theexpnaming convention.