Reduced density matrices

[lkdvos]

Replace the expectation value calculation for PEPS with computing tr(H * rho).

I checked the contraction costs and the leading cost is the same for both local operators, local norms and local density matrices. (In fact, the optimal contraction order for a single site local operator is exactly first computing the reduced densitymatrix and then contracting with the operator).
The main difference here is that instead of computing both the operator and the norm network, we can simply compute the density matrix once, divide it by its trace, and then evaluate the operator.

I'll try and include some benchmarks as well later, just to prove that everything is behaving as desired, but overall it seems to be between 25-40% faster to evaluate expectation values like this, and the memory usage scales accordingly as well.
Note that this also affects the gradient calculations, and even the compile times since the contraction functions are quite demanding for the compiler.

We should also try and hardcode the "easy" cases to optimize the intermediate permutations, but that might have to wait for a future PR.

[codecov]

Codecov Report

Attention: Patch coverage is 90.19608% with 10 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/contractions/localoperator.jl	89.01%	10 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/algorithms/toolbox.jl	98.29% <100.00%> (+0.01%)	⬆️
src/utility/util.jl	58.47% <100.00%> (+3.42%)	⬆️
src/algorithms/contractions/localoperator.jl	73.93% <89.01%> (-21.02%)	⬇️

🚀 New features to boost your workflow:

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

[Yue-Zhengyuan]

reduced_densitymatrix can be useful when users want to measure different physical quantities on the same sites (e.g. hopping energy, spin correlation, total energy on the same bond) by avoiding repeated contractions. I suggest to export it for end users.

[lkdvos]

Good for me, are we happy with the interface in that case?

[Yue-Zhengyuan]

Current interface is already good for me 👍 though one minor improvement is to add a default case in which bra = ket:

reduced_densitymatrix(inds, ket, env) = reduced_densitymatrix(inds, ket, ket, env)

[Yue-Zhengyuan]

In the same spirit, I think it's time to resolve the TODO in start_correlator?

PEPSKit.jl/src/algorithms/correlators.jl

Lines 83 to 102 in 47e6675

	# TODO: part of these contractions is duplicated between the two output tensors,
	# so could be optimized further
	@autoopt @tensor Vn[χSE Detop Debot; χNE] :=
	E_south[χSE Dstop Dsbot; χSW2] *
	C_southwest[χSW2; χSW] *
	E_west[χSW Dwtop Dwbot; χNW] *
	C_northwest[χNW; χN] *
	conj(bra(sandwich)[d; Dnbot Debot Dsbot Dwbot]) *
	ket(sandwich)[d; Dntop Detop Dstop Dwtop] *
	E_north[χN Dntop Dnbot; χNE]
	@autoopt @tensor Vo[χSE Detop dstring Debot; χNE] :=
	E_south[χSE Dstop Dsbot; χSW2] *
	C_southwest[χSW2; χSW] *
	E_west[χSW Dwtop Dwbot; χNW] *
	C_northwest[χNW; χN] *
	conj(bra(sandwich)[d1; Dnbot Debot Dsbot Dwbot]) *
	removeunit(O, 1)[d1; d2 dstring] *
	ket(sandwich)[d2; Dntop Detop Dstop Dwtop] *
	E_north[χN Dntop Dnbot; χNE]

The improvement is quite similar to what we are doing here with RDMs, although I haven't check the cost and memory usage there:

• We first contract Vo without inserting the physical operator O.
• Then Ve is obtained by taking the trace.
• Finally Vo is obtained by contracting with O.

BTW, can you inform me automated tools used to such checks? Anyway we can do this in the PR that resolves #231.

[lkdvos]

I'd have to think a bit more about that one though, it really does alter the contractions quite a bit unfortunately.

I'm not sure what you mean by automatic tools for checks, what kind of checks do you have in mind?

[pbrehmer]

Assuming that this PR is finished: Looks ready to go for me! Perhaps in the future (and in a different PR) we could have some example that uses the reduced_densitymatrix function to show how it can be used in action.

[lkdvos]

Lets merge #229 first, I'll update this one and see if I find the energy to actually test the performances a bit.

[ogauthe]

As a comment, if you plan to support PEPO CTMRG and tensors with ancilla in the future, you may consider implementing it here.
You can dispatch over leg number, add a dummy leg for PEPS tensor (=pure wavefunction) and then use a single contraction function, written with an ancilla. The dummy leg would not add any cost.

[lkdvos]

While I fully agree with this, unless someone feels like tackling that right now I might leave this for a future PR 😉