This post is about extracting target blocks of the Hamiltonian based on the flux of the initial state from which I want to make sparse array to speed up ED. My minimal example:
using SparseArray, ITensors
function heisenberg(n)
os = OpSum()
for j in 1:(n-1)
os += 1 / 2, "S+", j, "S-", j + 1
os += 1 / 2, "S-", j, "S+", j + 1
os += "Sz", j, "Sz", j + 1
end
return os
end
let n = 10
qn = true
s = siteinds("S=1/2", n; conserve_qns = qn)
H = MPO(heisenberg(n), s)
initstate(j) = isodd(j) ? "↑" : "↓"
ψ0 = random_mps(s, initstate; linkdims = 4)
cb = combiner(s...)
H_full = contract(H) * cb' * dag(cb)
ψ = contract(ψ0) * cb
effInd1 = findfirst(x -> x.first == flux(ψ), inds(H_full)[1].space)
effInd2 = nzblocks(ψ)[1].data[1] |> Int
@show effInd1 == effInd2
H_sp = sparse(ITensors.blockview(ITensors.tensor(H_full), nzblocks(H_full)[effInd1]))
ψ_eff = Array(ITensors.blockview(ITensors.tensor(ψ), nzblocks(ψ)[1]))
vals, vecs, info = @btime eigsolve(
$H_sp, $ψ_eff, 1, :SR; ishermitian = true, tol = 1e-6, krylovdim = 30, eager = true,
);
end
My Questions:
What is the right way to get the target symmetry sector? (My 2 approaches yield same results, named effInd1 and effInd2)
Is this the correct way to slice the target block and create a sparse array?
My approach is based from Slicing tensor with QN. Is this approach efficient and accurate for slicing the target block?
Naive Expectation: since the ψ contains only 1 non-zero block, explicitly extracting the effective block should not matter so much.
This approach is way faster than fusing the indices and eigsolve, even if I disable the sparsity:
Change H_sp = sparse(ITensors.blockview(ITensors.tensor(H_full), nzblocks(H_full)[effInd1]))
to H_arr = Array(ITensors.blockview(ITensors.tensor(H_full), nzblocks(H_full)[effInd1])))
I’m asking these questions because I couldn’t find functions like nzblocks documented from ITensors.jl or somewhere else. Maybe I just missed the part…
A lot of the functionality you are using is considered internal functionality that most users should not be using. A telltale sign that functionality is internal is that it is undocumented, not exported, or that you are accessing fields of types, such as x.space or x.data. You are of course free to use internal functionality, but be aware that internal functionality may change even in patch releases of ITensor without warning, and in fact I’m certain some of the functionality you are using will be changed in future releases because we are in the process of rewriting the entire NDTensors library ([NDTensors] Roadmap to removing `TensorStorage` types · Issue #1250 · ITensor/ITensors.jl · GitHub), though we don’t know the exact time frame for when that will be finished.
I might be able to guide you on implementing what you are trying to accomplish with documented, external functionality, however I’m not sure if that is worth the hassle for you or me since at the end of the day ITensor is meant more as a tensor network library rather than an exact diagonalization library. I think with the rewrite planned in [NDTensors] Roadmap to removing `TensorStorage` types · Issue #1250 · ITensor/ITensors.jl · GitHub as well as our release of non-abelian symmetries we will have more solid tools that we can make available for end users to use for performing exact diagonalization and analyzing symmetry sectors, though again the timeframe for that is unknown (I would definitely hope both of those would be completed within a year).
As an alternative I would suggest looking for devoted exact diagonalization libraries. Two that come up in searches or I was told about (but haven’t tried myself) are:
Something that could be helpful if you want to use ITensor’s OpSum interface and/or make it easier to compare ITensor’s DMRG results to exact diagonalization results would be to write a conversion tool from ITensor’s OpSum object to the symbolic operator format of those packages. That is something I have discussed with Alex Wietek, the author of the first package.
I’ll add that it is nice to see users exploring interesting applications of ITensor that push the envelope of the “standard” use cases like DMRG! Our goal is definitely to make it easier for users like you to reach down further below the “surface layer” of ITensor and do interesting things with block sparse tensors and other sparse and symmetric tensors, but unfortunately the library isn’t quite there yet. As I laid out in my post above, I think we have a good plan for a better design for NDTensors.jl that we can more comfortably document and expose to users, but that isn’t ready yet and mostly we just need the time and human resources to execute our plan for that part of the codebase.