Paradox of normalizing the biggest eigenvalue module to one

General Usage Performance
, ,
1

I run into an interesting but nasty problem in normalizing a matrix so that the largest eigenvalue module is one. The λ in convertMPS_PBC2OBC should print one, but I got values slightly different from one.

The minimum code that reproduces my problem is

using SparseArrays
using LinearAlgebra

const ⊗ = kron
const L = 64

function getMPS(L::Int, N::Int)
    #------------------------------------------------------------------------------
    As = [spzeros(Float64, 2, 2) for i in 1:3]
    As[1][1, 2] = sqrt(2.0/3.0)
    As[2][1, 1], As[2][2, 2] = -1/sqrt(3.0), 1/sqrt(3.0)
    As[3][2, 1] = -sqrt(2.0/3.0)
    #------------------------------------------------------------------------------
    M = Matrix{SparseMatrixCSC{Float64,Int64}}(undef, 3, 3)
    for n in 1:3
        for m in 1:3
            M[n, m] = spzeros(Float64, N+1, N+1)
            for i in 1:N
                M[n, m][i, i]   = (n == m ? (-1)^i : 0)
                M[n, m][i, i+1] = ((n == 3 && m == 1) ? 2*(-1)^i : 0)
            end
            M[n,m][N+1, N+1]  = (n == m ? (-1)^(N+1) : 0)
        end
    end
    #------------------------------------------------------------------------------
    AList = [spzeros(Float64, 2*(N+1), 2*(N+1)) for i in 1:3]
    for n in 1:3
        for m in 1:3
            AList[n] .+= M[n, m] ⊗ As[m]
        end
    end
    return AList
end

#-----------------------------------------------------------------------------#
# convert PBC MPS to OBC MPS format
#-----------------------------------------------------------------------------#
function convertMPS_PBC2OBC(AList::T) where T <: Vector{SparseMatrixCSC{Float64,Int64}}
    q = length(AList); D = 2
    AList_OBC = [Matrix(I*1.0, D, D) ⊗ AList[n] for n in 1:q]
    #--------------------------------------------------------------------------
    # normalize A so that the norm of ψ does not blow up
    #--------------------------------------------------------------------------
    # E = Σₘ AList_OBC[m] ⊗ AList_OBC[m] then find the largest eigenvalue
    E = sum([A ⊗ A for A in AList_OBC])
    eigvals, = eigen(Matrix(E))
    λ = maximum(abs.(eigvals))
    AList_OBC = AList_OBC/sqrt(abs(λ))   # AList_OBC = [A/sqrt(abs(λ)) for A in AList_OBC]
    E = sum([A ⊗ A for A in AList_OBC])
    eigvals, = eigen(Matrix(E))
    λ = maximum(abs.(eigvals))
    println("---", λ, "---")   # should be one, but actually not   
    return AList_OBC
end

function test()
    for N in 1:10
        AList = getMPS(L, N)
        AList_OBC = convertMPS_PBC2OBC(AList)
    end
end

test()
2

This is just floating point math not being exact. I don’t think this is a bug.

3

This a typical result

---1.0000000000000009---
---1.0000015290562512---
---0.9999997349719375---
---0.9999880258403627---
---0.9995615547185717---
---0.9996769383879271---
---0.997801171956432---
---0.9980697350888111---
---0.9999470000617328---
---1.0015127247553473---

As you can see, the difference can be as large as 0.002 which is unacceptable in floating point calculations.

4

Can you give a more reduced example? (For example a literal matrix for which this happens).

5

It seems that the problem is with the matrix list AList returned from getMPS(): The matrices are very singular and have significant numerical instability (assuming the Julia library implementation has no bug). If I use other matrices as shown in the following code,

using LinearAlgebra
const ⊗ = kron

function normalizeA()
    AList = [rand(ComplexF64, 16, 16) for _ in 1:3]
    E = sum([conj(A) ⊗ A for A in AList])
    eigvals, = eigen(Matrix(E))
    λ = maximum(abs.(eigvals))
    AList_new = [A/sqrt(abs(λ)) for A in AList]
    F = sum([conj(A) ⊗ A for A in AList_new])
    eigvals, = eigen(Matrix(F))
    λ = maximum(abs.(eigvals))
    println("---", λ, "---")
end

for i in 1:10
    normalizeA()
end

The problem does not pop out anymore.

6

Yeah, that makes sense. Julia is just calling to Blas/Lapack for these types of operations, so in general, the problem is in the conditioning of the matrix.

1 Like