`contract` for HermitianSeries

[antoine-levitt]

Hi @lxvm , many thanks for your package, it's very neat! I'm currently attempting to use it and need to contract a series along the first dimension, and hitting the error("not implemented") error. Can I ask if it would be easy to implement? Otherwise I tried converting the series to a non-hermitian one, contracting and converting back, but I can't seem to make it work; presumably I'd need to invert the operation c = f.c[CartesianIndices(ax_),_ax_[begin+div(length(_ax_),2):end],CartesianIndices(_ax)]? I'd also keep the original non-hermitian series around, but I'm getting it from AutoBZ.load_wannier90_data (re which also many thanks!) and I'm too scared to go into what's going on there (I spent many hours wrestling with wannier90 some years back and still have post-traumatic flashes)

[antoine-levitt]

(OK, I found out how to work around it by passing herm=false to AutoBZ)

[lxvm]

Hi @antoine-levitt, I think the workaround you already found is the only way to obtain a Fourier series with all of its coefficients using AutoBZ.jl. I made HermitianFourierSeries the default in AutoBZ.jl because it does speed up the Wannier interpolation ~x2 by halving the coefficient array along the innermost dimension as you noticed, although this speedup only applies if you contract the series along its last dimension so that contracting the series maintains memory contiguity. For constructing a FourierSeries from Wannier90 data outside of AutoBZ.jl, so you can experiment with it without worrying about existing optimizations, I think the simplest example is here

[antoine-levitt]

OK, thanks! It'd still be nice to just be able to contract along the first dimension, though: I'm interested in the case where I don't particularly care about the cost of contraction, but want to reuse the resulting fourierseries for more evals.

[lxvm]

Could you explain why you want to contract the series out of order with the numerical and type-system guarantees of a Hermitian output? If performance is not an issue you can already contract a non-Hermitian series out of order, or you can evaluate a Hermitian series point-wise without coding contractions yourself using the functor interface. Also taking a Hermitian view of the output of a non-Hermitian series should work well when the coefficient data has the Hermitian symmetry H_R = H_-R'. I mention at the beginning of the documentation that this package is intended for evaluation of series on hierarchical grids (with the hierarchy in starting from the outermost dimension and going inward because Julia is column-major) so to that end and for performance reasons the Hermitian series currently have some limitations.

The reason I can't easily change the code to contract the first dimension is because I implemented the Hermitian symmetry H_R = H_-R' by discarding coefficients for the negative R values along the first dimension. When I developed the package I also tested keeping coefficients for all R and discarding their upper triangle, eg using a regular series with SHermitianCompact coefficients, but this wasn't as fast. The latter approach is amenable to contracting any dimension but I chose the first for performance reasons and to not make assumptions about the types of the coefficients themselves.

I will document here that Hermitian series are limited to a standard order of contraction and also the link above showing how to get a series from Wannier90 data. In AutoBZ.jl I will document the herm keyword of load_wannier90_data to say that by default it returns a Hermitian series optimized for use in AutoBZ.jl but by setting herm=false a more user-friendly series is returned. However I will break that interface in AutoBZ.jl when I implement gauge-covariant derivatives, so if you are developing algorithms using these series I would suggest building series using the code snippet linked above where you can decide whether or not you need the series to be Hermitian.

[antoine-levitt]

That's fine, but you could still contract along the first dimension, build a new non-hermitian series then discard half of the (new) first dimension (old second one), right? But I understand that this may be a performance trap and it's an expert interface anyway, so working with non-hermitian series is reasonable. My application is that I want to evaluate on a hierarchical grid that may sometimes start from the first dimension (I can tell you about it in Lausanne if you're interested, it's pretty fun)

[lxvm]

That sounds like an interesting method and I'd like to hear about it in Laussane! I'd like to leave the issue open for now but I agree the change you described should work and I can implement it quickly in Laussane. Thank you for using the package!