DMRG runs slow

[mikpizza]

Hi all,
I am trying to perform the DMRG on a Transverse Field Ising Model of 8 spins. I'm coming from ITensors, and was pretty shocked to see the DMRG running that slow on Quimb. Also, sometimes the DMRG algorithm breaks, stopping its execution because of the energy having an imaginary part.

Question 1: what did I do wrong such that running DMRG is ineficcient and, sometimes, incorrect?

Here I post my function to build the Hamiltonian:

def tfim_mpo(N:int, J, g, pbc=True):
    """
    Build the transverse-field Ising model (TFIM) Hamiltonian as a SpinHam1D object in quimb.

    Parameters
    ----------
    N : int
        Number of sites in the 1D chain.
    J : float
        Coupling strength for the nearest-neighbor interaction (X_i X_{i+1}).
    g : float
        Transverse field strength for the on-site Z term.
    cyclic : bool, optional
        Whether to apply periodic boundary conditions (default True).

    Returns
    -------
    ham : quimb.tensor.tensor_1d_tebd.SpinHam1D
        Spin Hamiltonian object representing the TFIM Hamiltonian.
    """
    builder = SpinHam1D(S=1 / 2,  cyclic=pbc)
    builder += -4 * J, 'Z', 'Z' # the factor of 4 is due to 'Z' being 'S_z' and not just a Pauli z operator.
    builder += 2 * g, 'X'       # the factor of 2 is due to 'X' being 'S_x' and not just a Pauli x operator.
    return builder.build_mpo(N)

and the code I execute:

from quenches_utilities import *
import time
from qutip import *
import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":
    # parameters of the Hamiltonian
    N = 8
    J = 0.0
    h_x = -0.3
    g = -h_x

    # Build the Hamiltonian using Numba for efficiency
    H_mpo = tfim_mpo(N, J, g)
    dmrg = qtn.DMRG(H_mpo, bond_dims=[10, 20, 100, 100, 200], cutoffs=1e-8, bsz=2)

    dmrg.solve(tol=1e-6, max_sweeps = 20, verbosity=1)

The output I get is:

1, R, max_bond=(10/10), cutoff:1e-08
100%|#############################################| 8/8 [00:00<00:00,  9.73it/s]
Energy: -1.965925856914151 ... not converged.
2, R, max_bond=(10/20), cutoff:1e-08
100%|#############################################| 8/8 [00:01<00:00,  5.28it/s]
Energy: -3.624313843572727 ... not converged.
3, R, max_bond=(20/100), cutoff:1e-08
100%|#############################################| 8/8 [00:03<00:00,  2.21it/s]
Energy: -0.27179179958155036 ... not converged.
4, R, max_bond=(40/100), cutoff:1e-08
100%|#############################################| 8/8 [00:03<00:00,  2.31it/s]
Energy: -0.014501582237419066 ... not converged.
5, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:01<00:00,  6.76it/s]
Energy: -3.10069795645483e-05 ... not converged.
6, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:08<00:00,  1.08s/it]
Energy: -3.1655960431686883e-06 ... not converged.
7, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:04<00:00,  1.86it/s]
Energy: -1.5305246805731275e-07 ... not converged.
8, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:03<00:00,  2.64it/s]
Energy: -7.51583421841056e-09 ... converged!

Also, I am pretty worried about the bond dimension. In this model, when J=0 (as I specified in the code) the ground state is a product state.
Question 2: Why am I getting a state with such a big bond dimension instead?

The same code, in a different run, also produces the following error:

1, R, max_bond=(10/10), cutoff:1e-08
100%|#############################################| 8/8 [00:00<00:00,  9.99it/s]
Energy: -1.3550175495183845 ... not converged.
2, R, max_bond=(10/20), cutoff:1e-08
100%|#############################################| 8/8 [00:01<00:00,  5.15it/s]
Energy: -4.957536842639521 ... not converged.
3, R, max_bond=(20/100), cutoff:1e-08
100%|#############################################| 8/8 [00:02<00:00,  3.03it/s]
Energy: -4.9442928062058655 ... not converged.
4, R, max_bond=(40/100), cutoff:1e-08
100%|#############################################| 8/8 [00:05<00:00,  1.53it/s]
Energy: -2.468799611054961 ... not converged.
5, R, max_bond=(53/200), cutoff:1e-08
100%|#############################################| 8/8 [00:02<00:00,  3.78it/s]
Energy: -0.9298142204425284 ... not converged.
6, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:04<00:00,  1.67it/s]
Energy: -0.0604223659292594 ... not converged.
7, R, max_bond=(64/200), cutoff:1e-08
100%|#############################################| 8/8 [00:02<00:00,  2.83it/s]
Energy: -0.02059533778184725 ... not converged.
8, R, max_bond=(62/200), cutoff:1e-08
100%|#############################################| 8/8 [00:02<00:00,  3.01it/s]
Energy: (0.12260636648328359-1.501494942688413e-17j) ... not converged.
9, R, max_bond=(59/200), cutoff:1e-08
  0%|                                                     | 0/8 [00:00<?, ?it/s]Traceback (most recent call last):
  File "/home/mik/codespace/projects/Quenches/test.py", line 23, in <module>
    dmrg.solve(tol=1e-6, max_sweeps = 20, verbosity=1)
    ~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 1127, in solve
    energy = self.sweep(direction=direction, **sweep_opts)
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 957, in sweep
    self.ME_eff_norm = MovingEnvironment(self.TN_norm, **nm_opts)
                       ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 276, in __init__
    self.init_segment(begin, start, stop)
    ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 289, in init_segment
    self.init_non_segment(start, stop + self.bsz // 2)
    ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 354, in init_non_segment
    callback(start, stop, self.begin)
    ~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_dmrg.py", line 452, in cyclic_canonizer
    k.canonize_cyclic(slice(start, stop), bra=b, inv_tol=inv_tol)
    ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_1d.py", line 972, in canonize_cyclic
    kbc = kb.replace_section_with_svd(
        start,
    ...<6 lines>...
        rtags="_RIGHT",
    )
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_core.py", line 5885, in replace_section_with_svd
    return self.replace_with_svd(
           ~~~~~~~~~~~~~~~~~~~~~^
        where=slice(start, stop),
        ^^^^^^^^^^^^^^^^^^^^^^^^^
    ...<4 lines>...
        **replace_with_svd_opts,
        ^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_core.py", line 5842, in replace_with_svd
    TL, TR = tensor_split(
             ~~~~~~~~~~~~^
        A,
        ^^
    ...<9 lines>...
        rtags=rtags,
        ^^^^^^^^^^^^
    )
    ^
  File "/home/mik/miniconda3/lib/python3.13/functools.py", line 934, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^

  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/tensor_core.py", line 629, in tensor_split
    left, s, right = _SPLIT_FNS[method](array, **opts)
                     ~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/quimb/tensor/decomp.py", line 609, in isvd
    U, s, V = sli.svd(x, k)
              ~~~~~~~^^^^^^
  File "/home/mik/codespace/projects/Quenches/venv/lib/python3.13/site-packages/scipy/linalg/interpolative.py", line 927, in svd
    U, S, V = _backend.idzr_rsvd(A, k, rng=rng)
              ~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^
  File "_decomp_interpolative.pyx", line 1925, in scipy.linalg._decomp_interpolative.idzr_rsvd
ValueError: could not broadcast input array from shape (64,) into shape (256,)
  0%|                                                     | 0/8 [00:00<?, ?it/s]

[mikpizza]

It turns out the problem was just occuring in Periodic Boundary Conditions

[mikpizza]

Closing the discussion to not waste the time of others :)

[jcmgray]

Just to add a few pointers:

• As you found PBC DMRG is not so efficient, especially not with small L, it actually gets better as the length gets longer.
• Generally for 1D PBC, either using Simple Update (/ imaginary TEBD), or directly optimizing using autodiff will be more stable.

And as a side note, the _decomp_interpolative error was a bug in scipy scipy/scipy#22451 that should be fixed now scipy/scipy#22834.