[NDTensors] [BUG] Bug running DMRG on Metal with small cutoffs and large bond dimensions

using ITensors: MPO, OpSum, dmrg, random_mps, siteinds
using Metal: mtl
using Random: Random

function ham(n)
  os = OpSum()
  for j=1:n
    os += "Nup",j
    os += "Ndn",j
  end
  for j=1:n-1
    os += "Cdagup",j,"Cup",j+1
    os += "Cdagup",j+1,"Cup",j
    os += "Cdagdn",j,"Cdn",j+1
    os += "Cdagdn",j+1,"Cdn",j
  end
  return os
end

function main(; n, device, conserve_qns, nsweeps, maxdim)
  s = siteinds("Electron", n; conserve_nfparity=conserve_qns)
  Random.seed!(1234)
  ψ = device(random_mps(s, j -> isodd(j) ? "↑" : "↓"; linkdims=4))
  cutoff = [1e-3, 1e-13]
  outputlevel = 1
  H = device(MPO(ham(n), s))
  energy, ψ = dmrg(H, ψ; nsweeps, maxdim, cutoff, outputlevel)
end

Running either of these:

main(; device=mtl, conserve_qns=true, n=20, nsweeps=4, maxdim=200)
main(; device=mtl, conserve_qns=false, n=20, nsweeps=4, maxdim=200)

after some random small number of sweeps leads to the following error:

ERROR: operator does not appear to be hermitian: 0.0 vs NaN
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] expand!(iter::KrylovKit.LanczosIterator{ITensors.ProjMPO, ITensors.ITensor, KrylovKit.ModifiedGramSchmidt2}, state::KrylovKit.LanczosFactorization{ITensors.ITensor, Float32}; verbosity::Int64)
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/factorizations/lanczos.jl:250
  [3] eigsolve(A::ITensors.ProjMPO, x₀::ITensors.ITensor, howmany::Int64, which::Symbol, alg::KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2, Float64})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/lanczos.jl:76
...

It doesn't seem to occur at smaller bond dimensions, and also is fixed by user a larger cutoff, i.e.

  cutoff = [1e-3, 1e-7]

so seems to be related to the condition number of the MPS.