Eigenpair accessors

[pablosanjose]

Currently, if we have a sp = Spectrum(ham), we can do energies(sp) to obtain the eigenenergies and states(sp) to obtain a matrix with the eigenstates. If we have a bs = bandstructure(ham, mesh), we can do bands(bs) to obtain a vector of Band objects, energies(bs) to obtain a vector of all energies of all bands and states(bs, i) to obtain all the eigenvectors of band i.

This is a bit of a mess. It also doesn't allow us to directly access one specific eigenvalue of a spectrum or bandstructure. If we want to plot the spatial profile of an eigenstate close to a given energy at a given momentum of bs, there is no easy way.

I propose redesigning the way we extract eigenvalues and eigenvectors from Spectrum, Band and Bandstructure by defining a convenient getindex indexing syntax. The following possibility comes to mind:

Spectrum sp

sp[1] -> ϵ₁, ψ₁  # first eigenpair in sp, with ϵ₁::Number and ψ₁::Vector
sp[1:3] -> ϵv, ψm # tuple of ϵv::Vector and  ψm::Matrix of first three eigenvalues and eigenvectors, respectively
sp[around = 0.2] -> ϵᵢ, ψᵢ  # single eigenpair closest to energy 0.2
sp[around = (0.2, 10)] -> ϵv, ψm  # 10 eigenpairs (encoded as ϵv::Vector and  ψm::Matrix) closest to energy 0.2

Band b

b[(ϕᵢ...)] = ϵ, ψ  # interpolated eigenpair at Bloch phases (ϕᵢ...)

Bandstructure bs

bs[1] -> band₁  # equivalent to bs.bands[1]
b2[2:3] -> bandvector # equivalent to bs.bands[2:3]
bs[1, (ϕᵢ...)] -> ϵ, ψ  # equivalent to bs.bands[1][(ϕᵢ...)]
bs[2:3, (ϕᵢ...)] -> ϵv, ψm  # ϵv::Vector and  ψm::Matrix generalization of the above for bands 2 and 3
bs[(ϕᵢ...), around = 0.2] -> ϵ, ψ  # interpolated eigenpair of the band with interpolated energy closest to 0.2
bs[(ϕᵢ...), around = (0.2,10)] -> ϵv, ψm  # interpolated eigenpairs of the 10 bands with interpolated energies closest to 0.2

The trickiest bit here is b[(ϕᵢ...)], which requires identifying the simplices which contain (ϕᵢ...) inside, and then using linear interpolation of eigenpairs on each. With our current definition of multivalued bands, there can be several such simplices, and as a result b[(ϕᵢ...)] can yield several eigenenergies and eigenstates. In order to keep this type-stable, we might need to make b[(ϕᵢ...)] a Tuple{Vector,Matrix}, even though the energy vector will most often have a single element. As a result, also bs[1, (ϕᵢ...)] will yield a Tuple{Vector,Matrix}. The only methods that may return a single Tuple{Number, Vector} are sp[1], sp[around=0.2] and bs[(ϕᵢ...), around = 0.2]. It might be worth considering making all these method return Tuple{Vector,Matrix} instead.

[pablosanjose]

As a side comment: using the b[(ϕᵢ...)::Tuple] syntax leaves open the possibility to define b[energy::Number] in the future to yield a cut of a band at a given energy (i.e. Fermi surfaces), in the form of a new Band of one dimension less. This will require generalizing Band to meshes with dimension independent of embedding space dimension.

[pablosanjose]

@fernandopenaranda brings up a good point in slack: in the presence of degeneracies it might be best to have sp[around = (0.2, 10)] return 10 degenerate subspaces, not 10 states. In other words, it should return 10 or more eigenpairs so that they have 10 different energies in total (and no eigenpair with any such energy is left out). The simpler "10 states" behavior is not desirable because it would result in unpredictable and often non-reproducible results if any degenerate state is excluded from the 10 closest to energy = 0.2. It is reasonable to assume you always need full subspaces when there are degeneracies.

I would argue that sp[1] should still return a single eigenpair, despite any possible degeneracies. Also sp[1:4] should return strictly 4 eigenpairs. Why? Because the syntax is too similar to the AbstractVector one to implement a different behavior, I think. It would be too surprising.

So, this would be the revised syntax, with notation ϵ::Number, ψ::Vector and ϵv::Vector, ψm::Matrix

sp[1] -> ϵ, ψ  # first eigenpair in sp, regardless of degeneracies.
sp[1:3] -> ϵv, ψm # first three eigenpairs, regardless of degeneracies
sp[around = 0.2] -> ϵv, ψm  # eigenpairs of the single (possibly degenerate) subspace of energy closest 0.2.
sp[around = (0.2, 10)] -> ϵv, ψm  # 10 energy subspaces closest to energy 0.2

b[(ϕᵢ...)] = ϵv, ψm  # all (interpolated) eigenpairs of band b at Bloch phases (ϕᵢ...)
                     # Note: they can be empty if interpolation is not possible, i.e. no simplex contains (ϕᵢ...)

bs[1] -> band₁  # equivalent to bs.bands[1]
b2[2:3] -> bandvector # equivalent to bs.bands[2:3]
bs[1, (ϕᵢ...)] -> ϵv, ψm  # equivalent to bs.bands[1][(ϕᵢ...)]
bs[2:3, (ϕᵢ...)] -> ϵv, ψm  #  analogous to the above but for bs.bands[2:3]
bs[(ϕᵢ...), around = 0.2] -> ϵm, ψm  # interpolated eigenpairs at (ϕᵢ...) of energy subspace closest to 0.2
bs[(ϕᵢ...), around = (0.2,10)] -> ϵv, ψm  # interpolated eigenpairs at (ϕᵢ...) of 10 energy subspaces closest to 0.2

For the moment we will define "degenerate" as we do in bandstructures: energies within sqrt(eps(T)), where T is the number type of the energies. In the future we might want to add a degtol keyword to bandstructure and to the above getindex accessors.

[pablosanjose]

I fear that indexing a bandstructure with interpolation bs[(ϕᵢ...)] is often going to be slower (and less precise) than directly computing the eigenstates at that point (ϕᵢ...), except perhaps for very large systems (for which you are unlikely to compute a full bandstructure in the first place). The reason is that the naive search for (ϕᵢ...) within all simplices takes O(M) time, where M is the number of simplices. A benchmark of a typical situation with M=2600 and a 32x32 Hamiltoninan yields bs[(ϕᵢ...)] 27 times slower than diagonalization of the Hamiltonian at a new point.

The obvious way out is to try to leverage tree search algorithms to locate the relevant simplices in O(log(M)) time. I'm not talking about NearestNeighbors.jl, however, but rather something like IntervalTrees.jl or similar that can deal with intervals. So that's one more dependence :-(. IntervalTrees is light, fortunately, but I'd have to see how to implement N-dimensional intervals with it (because our simplices are not 1D).

[pablosanjose]

What a rabbit hole this is turning out to be! It happens that there is no higher-dimensional version of IntervalTrees.jl in the Julia ecosystem yet (that I could find). I did find a number of intriguing packages that might do do searches of higher dimensional interval sets (e.g. RegionTrees.jl, SpatialIndexing.jl,...), but they are somewhat heavy, and are actually too general for what we need here.

The fact is that our search intervals (defined by the bounding boxes of simplices) are not completely arbitrary. Simplices stitch together to form a mesh. Actually, in the current version, simplicesform a structured mesh after projecting out the energy coordinate. (Structured means that vertices are the direct product of given sets of grid values along each axis). I think we can use this structured property to great effect for efficiency, as long as we give up the possibility of having bands with non-structured meshes. Note that a future Dirac-cone point splitting strategy does not need to violate the structured assumption. Locally adaptive band refinement, however, would not be possible.

It's important to note that even if we assume structured grids for the k-space, bands are not single values currently (and probably never will, because of the complicated topologies that are possible -- and common), so they are not simply matrices of energies on the structured grid. If bands could be written as matrices, indexing into them b[(ϕs...)] would be very efficient, because we would just need to convert (ϕs...) to the closest Cartesian indices. However, a Band needs to be stored as a 1D list of vertices due to its general multivalued nature.

So, to start, we need a way to efficiently go (ϕ1, ϕ2,...) to the Cartesian index nϕ = (n1,n2,....) of a certain parent vertex on the structured k-grid (using naive search or IntervalTree on each axis), and from that to the subset of vertices that share that (ϕ1, ϕ2,...). Taking a cue from the SparseArrayCSC implementation, we could: (1) sort vertices by their lexicographic order, so that the desired vertex indices are contiguous; (2) build a tensor T_nϕ = jn:jn´ that returns the index interval in vertices of the desired vertices. This gives the indices in O(1) time. What we actually need are the simplices that contain the original (ϕ1, ϕ2,...). Since the first vertex of any simplex is the one with smallest index in the simplex, we can lexicographically sort simplices so that all simplices with a first indices with the same (ϕ1, ϕ2,...) are adjacent. Doing that, and defining the "parent" vertex jn above to be that of smallest index, should allow us to (3) build a map S_jn -> sjn:sjn´ that returns the range of simplices with first vertex index equal to jn.

In summary, by adding two mapping arrays per band and keeping vertices and simplices lexicographically ordered, we could achieve O(1) retrieval of the simplices that enclose a given (ϕ1, ϕ2,...). The memory overhead should not be that important, in relative terms, given that we already store one eigenstate per vertex. The code complexity worries me slightly, but I guess it's a matter of trying it.

[pablosanjose]

A completely different alternative: Bandstructure{D} is rewritten to contain

• base::MarchingMesh{D} that is a structured grid of (ϕs...) vertices. Purely lazy, doesn't store vertices, but can generate vertices, edges and simplicies programmatically. Implements interpolative indexing base[(ϕs...)] -> simplex (produces just one simplex, or nothing -- this is the base, not the lifted band mesh).
• energies::Array{Vector{T},D} with all eigenvalues at each base vertex
• states::Array{Matrix{T´},D} with all eigenvectors at each base vertex
• bands::Vector{Array{Int,D}}, where bands[i] is a matrix of pointers indicating the energies and states of a given band on each base vertex. If bands[i][base_vertex]==0, then the i-th band does not contain that base_vertex.

The bands tensors are constructed by a minor variation of extranct_bands.

There is no need to copy the results of diagonalize into a contiguous array: we just copy the pointer of energies and states into the base tensor element.

Computing the interpolated energies of a Bandstructure at a given (ϕs...) proceeds as follows:

(1) obtain the unique simplex base[(ϕs...)] as a set of indices ss of base vertices;
(2) for each band in bands compute pointers ptrs=(band[s_i] for s in ss);
(3) If no ptr is zero, extract ϵs = (energies[s][p] for (s,p) in zip(ss,ptrs)) and ψs = (view(states[s],:,p) for (s,p) in zip(ss,ptrs));
(4) construct the interpolated energies from all the ϵs;
(5) filter out whichever energy is not needed if around is given;
(6) Construct the interpolated states from the ψs corresponding to surviving ϵs

This should avoid all copying until the very end (6).

It also requires reimplementing all the Bandstructure and Mesh internals :-(. But it does look like an improvement of the current approach.

[pablosanjose]

Oh, I just realized that bands::Vector{Array{Int,D}} assumes single-valued bands. We'd need to iterate that.

[pablosanjose]

I moved the discussion of the redesign of Bandstructure internals to the new issue #122. That would need to be addressed first, before tackling this issue