Reproducibility issue: identify source of randomness

Hi! I'm using the approximate function with the VUMPS algorithm and it seems to be giving a different result every time. I'm not sure whether there is a call somewhere to a randomized function (maybe one of the eigensolvers?).

MWE:

using MPSKit, TensorKit
ψ₀ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^10, 2)) # random 2-site infinite MPS with bond dimsenion 10 that we want to truncate
I = DenseMPO([MPSKit.add_util_leg(id(storagetype(MPSKit.site_type(ψ₀)), physicalspace(ψ₀, i))) for i in 1:length(ψ₀)]) # corresponding identity operator as an MPO
ψ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^5, 2)) # approximate by MPS of bond dimension 5
alg = VUMPS(; verbose=false) # variational approximation algorithm
ψ1, = approximate(ψ, (I, ψ₀), alg)
ψ2, = approximate(ψ, (I, ψ₀), alg)
norm(ψ1.AL[1] - ψ2.AL[1]) # should be zero if ψ1==ψ2

gives something greater than zero, different every time.

Is there actually some source of randomness hidden inside? If there is, where? Would it be possible to control it with a random number generator from the top call to approximate ?

Is it a phase? Infinite states are difficult to define up to a phase, and this ambiguity starts already in the canonicalization function.

On Tue, 26 Mar 2024, 17:59 Stefano Crotti, @.***> wrote:

Hi! I'm using the approximate function with the VUMPS algorithm and it seems to be giving a different result every time. I'm not sure whether there is a call somewhere to a randomized function (maybe one of the eigensolvers?).

MWE:

using MPSKit, TensorKit ψ₀ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^10, 2)) # random 2-site infinite MPS with bond dimsenion 10 that we want to truncate I = DenseMPO([MPSKit.add_util_leg(id(storagetype(MPSKit.site_type(ψ₀)), physicalspace(ψ₀, i))) for i in 1:length(ψ₀)]) # corresponding identity operator as an MPO ψ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^5, 2)) # approximate by MPS of bond dimension 5 alg = VUMPS(; verbose=false) # variational approximation algorithm ψ1, = approximate(ψ, (I, ψ₀), alg) ψ2, = approximate(ψ, (I, ψ₀), alg)norm(ψ1.AL[1] - ψ2.AL[1]) # should be zero if ψ1==ψ2

gives something greater than zero, different every time.

Is there actually some source of randomness hidden inside? If there is, where? Would it be possible to control it with a random number generator from the top call to approximate ?

— Reply to this email directly, view it on GitHub https://github.com/maartenvd/MPSKit.jl/issues/135, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJKVCR4D2547QYSAYOOY33Y2GLM5AVCNFSM6AAAAABFJKM52KVHI2DSMVQWIX3LMV43ASLTON2WKOZSGIYDQOBQGQZDMMI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

I think the calls to rand happens in the initializations of the environments, and the initializations of the canonization routines. I think you can control the former, but the latter would be a bit more difficult. This being said, I think a better defined way of comparing states is using dot(psi1, psi2). This gives the fidelity between the two, which in your example should give you exactly 1. Alternatively (but I don;t think we have this implemented), you could get rid of some of the gauge freedom by SVD decomposing the C matrix and absorbing the unitaries in AL and AR respectively. Is there a particular reason you would need the results to be elementwise equal?

Is it a phase? Infinite states are difficult to define up to a phase, and this ambiguity starts already in the canonicalization function. … On Tue, 26 Mar 2024, 17:59 Stefano Crotti, @.> wrote: Hi! I'm using the approximate function with the VUMPS algorithm and it seems to be giving a different result every time. I'm not sure whether there is a call somewhere to a randomized function (maybe one of the eigensolvers?). MWE: using MPSKit, TensorKit ψ₀ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^10, 2)) # random 2-site infinite MPS with bond dimsenion 10 that we want to truncate I = DenseMPO([MPSKit.add_util_leg(id(storagetype(MPSKit.site_type(ψ₀)), physicalspace(ψ₀, i))) for i in 1:length(ψ₀)]) # corresponding identity operator as an MPO ψ = InfiniteMPS(fill(ℂ^2, 2), fill(ℂ^5, 2)) # approximate by MPS of bond dimension 5 alg = VUMPS(; verbose=false) # variational approximation algorithm ψ1, = approximate(ψ, (I, ψ₀), alg) ψ2, = approximate(ψ, (I, ψ₀), alg)norm(ψ1.AL[1] - ψ2.AL[1]) # should be zero if ψ1==ψ2 gives something greater than zero, different every time. Is there actually some source of randomness hidden inside? If there is, where? Would it be possible to control it with a random number generator from the top call to approximate ? — Reply to this email directly, view it on GitHub <#135>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJKVCR4D2547QYSAYOOY33Y2GLM5AVCNFSM6AAAAABFJKM52KVHI2DSMVQWIX3LMV43ASLTON2WKOZSGIYDQOBQGQZDMMI . You are receiving this because you are subscribed to this thread.Message ID: @.>

Indeed it's a phase

julia> abs.(ψ1.AL[1].data ./ ψ2.AL[1].data)
10×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

I think the calls to rand happens in the initializations of the environments, and the initializations of the canonization routines. I think you can control the former, but the latter would be a bit more difficult. This being said, I think a better defined way of comparing states is using dot(psi1, psi2). This gives the fidelity between the two, which in your example should give you exactly 1. Alternatively (but I don;t think we have this implemented), you could get rid of some of the gauge freedom by SVD decomposing the C matrix and absorbing the unitaries in AL and AR respectively. Is there a particular reason you would need the results to be elementwise equal?

I'm using VUMPS to truncate a uMPS as a subroutine within a (fully deterministic) iterative scheme. The output of the outer program fluctuates a lot in between runs, that's why i'm looking for the sources of randomness. Comparing the A tensors element-wise was just to show that the output is different every time. I'm not sure i need the output of VUMPS to be exactly the same every time, maybe up to some gauge will be enough, but to find out i need to investigate, and being able to fix the RNG would be of great help.

The use case is a bit unusual: i'm using a uMPS to directly parametrize a "classical" probability distribution normalized in the 1-norm (instead of a wavefunction normalized in the 2-norm). I'm guessing this is the reason why non-uniqueness up to unitary transformation actually gives different results every time, while it wouldn't matter in a quantum context.

Why is it difficult to control the randomness in the canonicalization functions?

The quickest solution will be to fix the gauge so that the 1-norm of your truncated wavefunction is a positive real. Beware with your usecase that the approximation that vumps makes is not "optimal", in the sense that the 2-norm difference is not a natural distance for probability distributions (it might be possible to do some gradient scheme instead on the KL-divergence between the original and approximated probability distributions?)

The randomness is problematic to get rid of because the environment tensors are defined through an eigenvalue problem, and the solutions of eigenvalue problems are only defined up to a phase. I don't think the canonicalization introduces another source of randomness in the AL/AR tensors, as we determine a new AL tensor deterministically and then solve AL C = C AR. The C (and AC) tensors are however only determined up to a phase...

Op wo 27 mrt 2024 om 10:55 schreef Stefano Crotti @.***

:

I think the calls to rand happens in the initializations of the environments, and the initializations of the canonization routines. I think you can control the former, but the latter would be a bit more difficult. This being said, I think a better defined way of comparing states is using dot(psi1, psi2). This gives the fidelity between the two, which in your example should give you exactly 1. Alternatively (but I don;t think we have this implemented), you could get rid of some of the gauge freedom by SVD decomposing the C matrix and absorbing the unitaries in AL and AR respectively. Is there a particular reason you would need the results to be elementwise equal?

I'm using VUMPS to truncate a uMPS as a subroutine within a (fully deterministic) iterative scheme. The output of the outer program fluctuates a lot in between runs, that's why i'm looking for the sources of randomness. Comparing the A tensors element-wise was just to show that the output is different every time. I'm not sure i need the output of VUMPS to be exactly the same every time, maybe up to some gauge will be enough, but to find out i need to investigate, and being able to fix the RNG would be of great help.

The use case is a bit unusual: i'm using a uMPS to directly parametrize a "classical" probability distribution normalized in the 1-norm (instead of a wavefunction normalized in the 2-norm). I'm guessing this is the reason why non-uniqueness up to unitary transformation actually gives different results every time, while it wouldn't matter in a quantum context.

— Reply to this email directly, view it on GitHub https://github.com/maartenvd/MPSKit.jl/issues/135#issuecomment-2022348939, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJKVCRHZCTIV4GQNSQUVYTY2KCPDAVCNFSM6AAAAABFJKM52KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMRSGM2DQOJTHE . You are receiving this because you commented.Message ID: <maartenvd/MPSKit .@.***>

Beware with your usecase that the approximation that vumps makes is not "optimal", in the sense that the 2-norm difference is not a natural distance for probability distributions (it might be possible to do some gradient scheme instead on the KL-divergence between the original and approximated probability distributions?)

I agree that something KL-based would be ideal, it's just not that well suited for MPSs. It would require to compute averages of the log of a matrix product, which can't be simplified analytically and therefore, AFAIK, the best one can do is to approximate it by sampling (of course i'd be happy to know about a better solution!). Requiring that the truncated distribution be close to the original one in the 2-norm still makes sense and can be tackled analytically to a larger extent.

The randomness is problematic to get rid of because the environment tensors are defined through an eigenvalue problem

I see. It still looks like it's possible to propagate a random number generator as a keyword argument from the innermost call to rand(), up to outer calls to e.g. VUMPS(), isn't it? I just need help locating the call to rand in the code. I'd be happy to make a PR if people are interested, or otherwise just keep it in my local version of the package.

If all you want is to control the randomness, all calls to the random number generator should be generated by this line https://github.com/maartenvd/MPSKit.jl/blob/91af3a2ad6802eedfbcb9579dc0137c8fbd14126/src/utility/utility.jl#L137C1-L137C48

Every time mpskit needs a random tensor, we use calls like similar() to construct a tensor with the correct spaces/storagetype and only then randomize the entries.

Op wo 27 mrt 2024 om 12:23 schreef Stefano Crotti @.***

:

Beware with your usecase that the approximation that vumps makes is not "optimal", in the sense that the 2-norm difference is not a natural distance for probability distributions (it might be possible to do some gradient scheme instead on the KL-divergence between the original and approximated probability distributions?)

I agree that something KL-based would be ideal, it's just not that well suited for MPSs. It would require to compute averages of the log of a matrix product, which can't be simplified analytically and therefore, AFAIK, the best one can do is to approximate it by sampling (of course i'd be happy to know about a better solution!). Requiring that the truncated distribution be close to the original one in the 2-norm still makes sense and can be tackled analytically to a larger extent.

The randomness is problematic to get rid of because the environment tensors are defined through an eigenvalue problem

I see. It still looks like it's possible to propagate a random number generator as a keyword argument from the innermost call to rand(), up to outer calls to e.g. VUMPS(), isn't it? I just need help locating the call to rand in the code. I'd be happy to make a PR if people are interested, or otherwise just keep it in my local version of the package.

— Reply to this email directly, view it on GitHub https://github.com/maartenvd/MPSKit.jl/issues/135#issuecomment-2022529818, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJKVCSPPUZUPGDASVKHSZDY2KM3NAVCNFSM6AAAAABFJKM52KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMRSGUZDSOBRHA . You are receiving this because you commented.Message ID: <maartenvd/MPSKit .@.***>

QuantumKitHub / MPSKit.jl

Reproducibility issue: identify source of randomness #135