Support efficient products of `AbstractHamiltonian`s and non-uniform excitation generation

[joachimbrand]

The problem with the current AbstractHamiltonian interface

The interface for AbstractHamiltonian and AbstractObservable require amongst other things

• given an address be able to calculate the number of non-zero offdiagonals in the column with num_offdiagonals, and
• deterministically generate one of these offdiagonals (address and value) with get_offdiagonal, which essentially indexes each one with an integer.

This information is conveniently used for uniformly sampling the offdiagonals, or for generating the full column for deterministic matrix-vector products. A problem arises for operators where it is costly to calculate the number of non-zero off-diagonals ahead of time. For example this is the case if the operator is defined as the product of two operators A*B, where calculating num-offdiagonals accurately would require to first generate all offdiagonals from B and then calling num_offdiagonals on the resulting addresses with A. Similarly, for fermionic Hamiltonians it is often complicated and impractical to accurately determine the number of offdiagonals and algorithmically access them.

Avoiding the problem

Let us consider what is really required for Rimu's functionality:

1. For deterministic matrix-vector products is is necessary to generate an exhaustive sequence of offdiagonals (address and value) in addition to the diagonal matrix element to define the operator column defined by a given address.
2. For QMC or stochastic matrix-vector products is is necessary to generate offdiagonal address-value pairs with a known probability (in a WithReplacement fashion).
3. For the DynamicSemiStochastic spawning strategy it would be good to know the length of the column, i.e. num_offdiagonals - but a reasonable estimate may be sufficient.

These requirements are less than what is currently required. The first two requirements could be easily met for product or fermionic hamiltonians without calculating the number of offdiagonals ahead of time. Moreover, if num_offdiagonals is replaced by a cheaply available estimate then this could also be cheap to calculate for products of operators.

Suggested solution: An OperatorColumn type and constructor

OperatorColumn could both be an abstract type and constructor for a type instance (with potentially a custom type suitable for the type of Hamiltonian). Before computing any of the information in points 1-3 one would instantiate a struct

column = OperatorColumn(hamiltonian, address)

that could precompute and store the diagonal element (because that is needed in any scenario), the estimated column length, and other potentially relevant information (e.g. convert the address into an OccupiedModeMap for efficient generation of excitations). Then iterating over column would generate all address, value pairs (e.g. with the diagonal element as the first one) and offdiagonals can be sampled from the column:

addr, diagonal_matrix_element = first(column)   # first element from iteration
column_dvec = DVec(column)                               # fill a DVec by iterating `column`
new_address, probability, value = random_offdiagonal(column) # similar to the existing function but now acting on the `column` struct
estimated_length = num_offdiagonals(column) # now just an estimate that could be the same for every column

This suggests changing the interface for AbstractHamiltonian and AbstractOperator to require them to provide a method for OperatorColumn that instantiates a struct with the properties mentioned above. In this case, the Hamiltonian would "own" the excitation generator, i.e. the algorithm that provides the random selection of offdiagonals.

For existing Hamiltonians that already implement get_offdiagonal and the rest of the existing interface, default implementations for OperatorColumn can be provided. For new Hamiltonians, implementing get_offdiagonal would no longer be required.

In order to facilitate easy-to-compute estimates for num_offdiagonals, operator traits like OneBodyOperator and TwoBodyOperator could be introduced. Setting the trait would then result in a value of N*M, and N(N-1)*M^2 for BoseFS addresses, respectively, where N is the particle number and M the number of modes. Corresponding pre-computed values (possibly even computed at compile time) for other types of Fock addresses can be provided.

Required changes to existing code

Code related to spawning and for matrix-operator products and three-way dot products would have to be adapted to make use of OperatorColumn and be cleaned of the old interface functions. The SpawningStrategys WithoutReplacement and Bernoulli would either be completely deprecated, or the strategy could be passed as an additional argument into the random_offdiagonal function. In this case the sampler would have to take care of this.

Checklist

• New AbstractOperator interface: operator_column function and AbstractOperatorColumn type. diagonal_element, starting_address, random_offdiagonal, and offdiagonals operate on a column. Operator Column #317
• Implement OperatorProduct.
• Implement OperatorSum.
• Operator traits control the availability of random_offdiagonal and the offdiagonals iterator.
• Move the responsibility for StochasticStyle, SpawningStrategy, and CompressionStrategy out of the AbstractDVec and into the PMCAlgorithm.
• DynamicSemistochastic type memoizes column elements.

[mtsch]

I agree with most of the points here. Some thoughts:

• I've experimented with having a column struct that stores the diagonal element in the past as it's currently annoying that some information (usually the OccupiedModeMaps) needs to be computed twice for each address.
• Letting the column handle random generation is a good idea too as it would allow us to e.g. pick smaller elements less frequently (like with three-body terms in Transcorrelated1D) and is probably the only way to reasonably implement stuff like operator products.
• In some cases, we should be able to iterate over all offdiagonals more efficiently than calculating each of them one by one.
• Regarding WithoutReplacement and Bernoulli - I'm in favor of dropping them altogether as they don't provide any real benefit and I'm pretty sure nobody has touched them in a long time. If we decide to implement customisable samplers, I suggest we worry about them when we actually need them as to not overgeneralise the code before we actually know what generalisations we need.
• As we do this it might be worth packing some common functionality together. For example, we currently have HubbardReal1D, HubbardReal1DEP, ExtendedHubbardReal1D which all produce exactly the same offdiagonals.

[joachimbrand]

Sum of operators

If num_offdiagonals is replaced by max_offdiagonals(operator) to give an upper bound on the number of off-diagonal non-zero matrix elements then we could define a sum of operators A + B by

OperatorSum(A::AbstractOperator, B::AbstractOperator; weights = (max_offdiagonals(A), max_offdiagonals(B)))

where relative weights for random spawning can default to the respective number of max_offdiagonals. When a spawning attempt is made, the first step would be to decide to spawn from A with probability weights[1] / sum(weights), or to spawn from B otherwise.

The OperatorSum could be generalised to create arbitrary linear combinations of the form α * A + β * B, which would help streamline the definition of similar Hamiltonians.

For the number of off-diagonals we would simply have

max_offdiagonals(A + B) = max_offdiagonals(A) + max_offdiagonals(B)

Product of operators

For a product A * B = OperatorProduct(A, B) the drawing of a random offdiagonal could proceed simply with (in pseudocode)

function random_offdiagonal(OperatorColumn(A * B, address))
    b_address, b_element, b_prob = random_offdiagonal(OperatorColumn(B, address))
    a_address, a_element, a_prob = random_offdiagonal(OperatorColumn(A, b_address))
    return a_address, a_element * b_element, a_prob * b_prob
end

If the sum of probabilities of possible spawns b_prob and a_prob separately adds to 1 (as they should), then so should the products generated here.

Special handling would have to be applied for the diagonal element but this could be efficiently automated for the square of IsHermitian operators, where we know that every element of the column also defines the return path:

diagonal_matrix_element(B * B, address) = sum(OperatorColumn(B, address)) do b_address, b_element
    conj(b_element) * b_element
end

In this case the value of any offdiagonal that returns to the address should be set to zero to avoid double counting.

In the case of a general product of operators it may the easiest to not accurately calculate the diagonal element but instead define

diagonal_matrix_element(B * B, address) = first(OperatorColumn(A, address))[2] *  first(OperatorColumn(B, address))[2]

and allow for the fact that some of the off-diagonals will actually return to address. This will still result in accurate deterministic evaluations of the matrix-vector products and correct (but slightly noisier) stochastic ones.

[joachimbrand]

As suggested by @mtsch, the abstract type OperatorColumn could be replaced by duck-typing. In this case, the interface would instead demand that a method for operator_column(::AbstractOperator, ::AbstractFockAddress) is defined that returns a suitable object column = operator_column(hamiltonian, address), which has the properties previously demanded for OperatorColumn.

One advantage of that approach would be to be able to subtype column to AbstractVector and thus get the iteration functionality for free. A possible disadvantage is losing the ability to dispatch on the type ::OperatorColumn when building any functions that work with operator columns.

[joachimbrand]

Spawning strategies

With the proposed changes for the operator column, the responsibility for stochastically spawning off-diagonal elements of the column (excitation generation) is moved to the column struct, while it currently sits with the SpawningStratgy of the StochasticStyle inside the DVec or PDVec container. Along with this change, it would make sense to remove the StochasticStyle from DVec and PDVec and move the remaining functionality to the PMCAlgorithm, accessed by thealgorithm keyword argument of ProjectorMonteCarloProblem. Currently the only implemented algorithm is FCIQMC.

• The detail of spawning covered by WithReplacement, WithoutReplacement, Bernoulli should be handled by the operator column. It might be best to not worry about implementing all options and just implement what is most efficient.
• IsStochasticInteger could be removed and the corresponding algorithm simply triggered by integer number types for the DVec/PDVec container. (Alternatively, it could be removed completely and integer number types disallowed.)
• IsDynamicSemistochastic, Exact, IsStochasticWithThreshold should be handled by the PMCAlgorithm

Memoizing operator columns for DynamicSemistochastic

This is a suggestion for a possible code optimisation that could be added later.

In the DynamicSemistochastic spawning strategy, the full operator column is generated and multiplied with the coefficient if the magnitude of the coefficient is large enough (typically larger than the number of offdiagonals in the column) rather than spawning off-diagonals stochastically. As this will typically need to be done repeatedly for the same columns, it would make sense to compute the elements of the column only once and store the results to retrieve them later from storage (in the sense of memoizing the value of a function, see e.g. [ANN] Memoization.jl).

For the operator columns this could be achieved by a MemoizingOperator type that wraps the Hamiltonian and would compute and store those operator columns (e.g. in a dictionary) when all values are requested. The type instance could be generated during init(), i.e. while initialising a PMCSimulation, and passed on to the state structs instead of the Hamiltonian. This should work well with MPI and distribute the storage of the deterministic part of the Hamiltonian, as only the operator columns relevant to the addresses assigned to the MPI rank will be stored. Care has to be taken to assure thread safety when storing new columns (i.e. lock the dictionary).

[joachimbrand]

Reopening this issue as there are some open tasks left.