# Paradox of normalizing the biggest eigenvalue module to one

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()
``````

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

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.

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

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.

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