Global Subspace Expansion

Hi everyone,

I am simulating the real-time dynamics of a 1D Transverse Field Ising Model (TFIM) using 1-site TDVP. The Hamiltonian is given by:

H = 4J \sum_{j=1}^{N-1} \hat{S}^z_j \hat{S}^z_{j+1} + 2\Gamma \sum_{j=1}^N \hat{S}^x_j

I am computing two time-dependent observables: 
  1. The Spin Autocorrelation Function:

    C(t) = \frac{1}{N} \sum_{j=1}^N \zeta_j(0) \cdot 2 \langle \psi(t) | \hat{S}^z_j | \psi(t) \rangle

    where \\zeta_j(0) = \\pm 1 represents the initial classical spin configuration at site j.

  2. The Mid-Chain Von Neumann Entanglement Entropy:

S_{\text{mid}}(t) = -\text{Tr}(\rho_A \ln \rho_A) = -\sum_k p_k \ln p_k

Because 1-site TDVP cannot dynamically increase the bond dimension (\\chi), I am using the Global Subspace Expansion (GSE) algorithm via expand(..., alg="global_krylov") intermittently during the time loop to adaptively grow the Krylov subspace.

My questions are:

  1. Is this inline deployment of expand mathematically sound when evaluating history-dependent observables like C(t) and bipartite entanglement entropy S(t)?

  2. Does the state perturbation introduced during the expand step cause an unphysical jump or artificial dissipation in the Schmidt spectrum that distorts the true Von Neumann entropy tracking?

Below is a working example showcasing my exact measurement and expansion architecture:

using ITensors
using ITensorMPS
using LinearAlgebra

# ========================================================
# 1. MID-CHAIN ENTANGLEMENT ENTROPY
# ========================================================
function entropy_mid(psi, mid)
    psi_copy = copy(psi)
    orthogonalize!(psi_copy, mid)
    _, S = svd(psi_copy[mid], (linkind(psi_copy, mid - 1), siteind(psi_copy, mid)))
    svals = diag(S)
    p = abs2.(svals)
    p = p[p .> 1e-14]
    isempty(p) && return 0.0
    p ./= sum(p)
    return -sum(p .* log.(p))
end

# ========================================================
# 2. GLOBAL SUBSPACE EXPANSION (GSE)
# ========================================================
function gse_expand(psi, H; krylovdim=2, gse_cutoff=1e-8, gse_maxdim=128)
    psi = expand(
        psi,
        H;
        alg="global_krylov",
        krylovdim=krylovdim,
        cutoff=gse_cutoff,
        apply_kwargs=(; maxdim=gse_maxdim)
    )
    orthogonalize!(psi, 1)
    normalize!(psi)
    return psi
end

# ========================================================
# 3. MAIN EXECUTION BLOCK (MWE)
# ========================================================
function main()
    # System Parameters
    N = 12
    T_max = 20.0       # Shortened for quick forum reproduction
    dt = 0.1
    J = 1.0
    Gamma = 1.0
    
    # MPS/TDVP Truncation Settings
    cutoff = 1e-8
    maxdim = 128
    
    # GSE Settings
    use_gse = true
    gse_every_steps = 20
    gse_maxdim = 128
    gse_cutoff = 1e-8
    krylovdim = 2

    # Construct Site Indices and Hamiltonian (Pure TFIM, No hz potential)
    sites = siteinds("S=1/2", N)
    os = OpSum()
    for j in 1:(N - 1)
        os += 4.0 * J, "Sz", j, "Sz", j + 1
    end
    for j in 1:N
        os += 2.0 * Gamma, "Sx", j
    end
    H = MPO(os, sites)

    # Initial State: Fully aligned Ferromagnetic State |Up Up ... Up>
    psi = MPS(sites, ["Up" for j in 1:N])
    z0 = ones(N) # Reference vector for autocorrelation calculation
    
    t_steps = collect(0.0:dt:T_max)
    n_steps = length(t_steps)
    mid = div(N, 2)

    println("Starting TDVP time evolution...")
    
    for i in 1:n_steps
        t = t_steps[i]
        
        # Measurement Step (Every 10 steps)
        if i == 1 || mod(i - 1, 10) == 0 || i == n_steps
            expZ = real.(expect(psi, "Sz"))
            Cval = sum(j -> z0[j] * 2.0 * expZ[j], 1:N) / N
            Sval = entropy_mid(psi, mid)
            chi = maxlinkdim(psi)
            
            println("t=$(round(t, digits=2)) | C(t)=$(round(Cval, digits=4)) | Entropy=$(round(Sval, digits=4)) | Max Bond Dim=$chi")
        end

        if i == n_steps
            break
        end

        # Dynamic Bond Expansion step
        if use_gse && mod(i, gse_every_steps) == 0
            psi = gse_expand(psi, H; krylovdim=krylovdim, gse_cutoff=gse_cutoff, gse_maxdim=gse_maxdim)
        end

        # Single-site TDVP time step
        psi = tdvp(H, -im * dt, psi; nsteps=1, nsite=1, cutoff=cutoff, maxdim=maxdim, normalize=true)
    end
    println("Evolution completed successfully.")
end

main()

Any insights regarding best practices for combining GSE with 1-site TDVP measurements would be highly appreciated!

Thanks in advance.

Expand only acts on the gauge degrees of freedom of the MPS and should leave all the physical properties unchanged. So if \ket{\psi} is the state before calling expand and \ket{\psi'} is the state afterwards then \braket{\psi | \psi'} = 1.

Consider the straightforward two site MPS for the state \ket{\psi} = \ket{00} in the “matrix of vectors” format

\ket{\psi} = \begin{pmatrix} \ket{0} \\ \end{pmatrix} \begin{pmatrix} \ket{0} \\ \end{pmatrix}

If you are restricted to modifying a single site tensor at a time then the state will always be a product state. If however you represent the state as

\ket{\psi'} = \begin{pmatrix} \ket{0} & \ket{1} \\ \end{pmatrix} \begin{pmatrix} 1 \ket{0} \\ 0 \ket{1} \\ \end{pmatrix}

then by modifying the second site tensor you can exactly represent any two-qubit state. Yet \ket{\psi} and \ket{\psi'} are physically indistinguishable.

Thank you for your reply @corbett5 ,

I have another question about this:

Does using expand function helps us keep a smaller bond dimension and get accurate results for real time evolution? I am interested in long time evolution of MPS. In order to do this the entanglement entropy grows very fast. Is there a way to control this bond dimension using GSE?

GSE is a way to get more accurate simulation at a cheaper cost, namely using one site TDVP instead of two site. Since one site TDVP does not modify the bond dimensions at all, if you are doing repeated GSE and TDVP you will likely want to truncate the state before performing GSE to slow the growth of the bond dimension. But at the end of the day real time evolution rapidly increases the entanglement, and no MPS technique is going to allow you to efficiently capture these highly entangled states.