ArgumentError: matrix contains Infs or NaNs in constructing MPO

[lin.zhang]

Hello,

I am trying to construct a density matrix with a weak domain wall-like imbalance in magnetization for the spin 1 system using MPO. The code is as follows:

using LinearAlgebra
using ITensors, ITensorMPS

L = 340
μ = 0.001

sites = siteinds(“S=1”, L)

rho = MPO(sites, n → (A = (n ≤ L ÷ 2 ? [1.0 + μ, 1.0, 1.0 - μ] : [1.0 - μ, 1.0, 1.0 + μ]);
A = A / sum(A);
diagm(0 => A)))

This code works for L <= 339 but breaks when L >= 340, and the following errors appear:

ArgumentError: matrix contains Infs or NaNs

Stacktrace:
[1] chkuplofinite(A::Matrix{Float64}, uplo::Char)
@ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.10.4+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:98
[2] syevr!(jobz::Char, range::Char, uplo::Char, A::Matrix{Float64}, vl::Float64, vu::Float64, il::Int64, iu::Int64, abstol::Float64)
@ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.10.4+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:5222
[3] eigen!
@ ~/.julia/juliaup/julia-1.10.4+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/symmetriceigen.jl:8 [inlined]
[4] eigen(A::Hermitian{Float64, Matrix{Float64}}; sortby::Nothing)
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/symmetriceigen.jl:13
[5] eigen
@ ~/.julia/juliaup/julia-1.10.4+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/symmetriceigen.jl:11 [inlined]
[6] eigen
@ ~/.julia/packages/NDTensors/gdQSG/src/lib/Expose/src/functions/linearalgebra.jl:24 [inlined]
[7] svd_recursive(M::Transpose{Float64, Matrix{Float64}}; thresh::Float64, north_pass::Int64)
@ NDTensors ~/.julia/packages/NDTensors/gdQSG/src/linearalgebra/svd.jl:40
[8] svd_recursive(M::Transpose{Float64, Matrix{Float64}})
@ NDTensors ~/.julia/packages/NDTensors/gdQSG/src/linearalgebra/svd.jl:29
[9] svd_recursive(M::Matrix{Float64}; thresh::Float64, north_pass::Int64)
@ NDTensors ~/.julia/packages/NDTensors/gdQSG/src/linearalgebra/svd.jl:32
[10] svd_recursive
@ ~/.julia/packages/NDTensors/gdQSG/src/linearalgebra/svd.jl:29 [inlined]
[11] svd(T::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}; mindim::Nothing, maxdim::Nothing, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, alg::Nothing, min_blockdim::Nothing)
@ NDTensors ~/.julia/packages/NDTensors/gdQSG/src/linearalgebra/linearalgebra.jl:108
[12] svd(A::ITensor, Linds::Vector{Index{Int64}}; leftdir::Nothing, rightdir::Nothing, lefttags::ITensors.TagSets.GenericTagSet{BitIntegers.UInt256, 4}, righttags::Nothing, mindim::Nothing, maxdim::Nothing, cutoff::Float64, alg::Nothing, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, min_blockdim::Nothing, utags::ITensors.TagSets.GenericTagSet{BitIntegers.UInt256, 4}, vtags::Nothing)
@ ITensors ~/.julia/packages/ITensors/oOwvi/src/tensor_operations/matrix_decomposition.jl:162
[13] truncate!(::NDTensors.BackendSelection.Algorithm{:frobenius, @NamedTuple{}}, M::MPO; site_range::UnitRange{Int64}, kwargs::@Kwargs{cutoff::Float64})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/abstractmps.jl:1692
[14] truncate!
@ ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/abstractmps.jl:1679 [inlined]
[15] truncate!(M::MPO; alg::String, kwargs::@Kwargs{cutoff::Float64})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/abstractmps.jl:1676
[16] MPO(::Type{Float64}, sites::Vector{Index{Int64}}, ops::Vector{Matrix{Float64}})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/mpo.jl:87
[17] MPO(::Type{Float64}, sites::Vector{Index{Int64}}, fops::var"#149#150")
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/mpo.jl:94
[18] MPO(sites::Vector{Index{Int64}}, ops::Function)
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/oOwvi/src/lib/ITensorMPS/src/mpo.jl:97

From the definition of density matrix (just being a tensor product of local density matrices), there should be no Infs and NaNs. But somehow they appear in constructing the MPO. Do you know how to fix this? Thanks very much in advance for your help.

[jzHan]

In numerical calculations, I believe it is advisable to avoid using such small parameters as much as possible, otherwise there is a high possibility of generating infinite numbers. I suggest you try increasing the \mu.

[miles]

Hi Lin,
I found the issue and it’s actually a bug or weakness in one of the ITensor MPO constructors. Please use the following code below as a workaround for your case. The code below calls a different constructor that avoids the issue.

using LinearAlgebra
using ITensors, ITensorMPS

let
  L = 340
  μ = 0.001

  sites = siteinds("S=1", L)

  function makeA(n)
    A = (n ≤ L ÷ 2 ? [1.0 + μ, 1.0, 1.0 - μ] : [1.0 - μ, 1.0, 1.0 + μ])
    A = A / sum(A)
    return diagm(0 => A)
  end

  ops = [makeA(n) for n=1:L]
  op_prod = Prod{Op}()
  for n in 1:L
    op_prod *= Op(ops[n], n)
  end
  rho = MPO(op_prod, sites)

  return
end

For more details about what the issue is, one of our ITensor MPO constructors does a numerical truncation of the MPO after constructing it. For some types of MPOs, even ones which contain all reasonably-sized numbers, doing a standard MPS-like truncation can lead to extremely large numerical values due to the difference between how states and operators are normalized. So really we should be truncating the MPO in a different way in that code. Please for now use the style of the code above.

Thanks, Miles

[lin.zhang]

Hi Miles,

Thanks very much for the code. It works!

Best wishes,
Lin