NaN Values for j > 1700 in Quantum Fisher Information Calculation using QuantumOptics

I’m calculating the Quantum Fisher Information (QFI) for the kicked top system using the QuantumOptics package in Julia. The code works fine for initial states with j up to 1700, but when I increase j to 1800 or higher, I start getting NaN values in the results. Even the norm of the initial state comes as NaN. Below is the code I am using:

#The QFI of the kicked tip system using the QuantumOptics pkg and the Krylov subspace

using QuantumOptics
using LinearAlgebra
using KrylovKit  
using Expokit

# Define the system parameters
j = 1800
J = SpinBasis(j)
Jx = 0.5 * sigmax(J)
Jy = 0.5 * sigmay(J)
Jz = 0.5 * sigmaz(J)

# Initial state
initial_state = coherentspinstate(J, Ο€/2, Ο€/2)

# Ensure the initial state is properly normalized
println("Initial state norm: ", norm(initial_state))

Ξ± = Ο€/2   # Precession angle
k = 7  # Kicking strength for chaotic regime
Ξ΅ = 0.001  # Small perturbation

# Construct the Hamiltonians
H_precession = Ξ± * Jz
H_kick = (k / (2*j)) * Jy^2


# Safe normalization function using BigFloat
function safe_norm(state)
    nrm = norm(state)
    if nrm < BigFloat(1e-10)
        return BigFloat(1e-10)
    end
    return nrm
end

#check if your Hamiltonians contain any NaNs or Infs
function check_matrix_validity(H)
    if any(isnan, H) || any(isinf, H)
        println("Hamiltonian Matrix contains NaNs or Infs")
        return false
    else
        println("Hamiltonian Matrix is valid")
        return true
    end
end
# Check the validity of the Hamiltonians
is_valid = check_matrix_validity(sparse((H_kick + H_precession).data)) &&
           check_matrix_validity(sparse((H_kick + (Ξ± + Ξ΅) * Jz).data))

if !is_valid
    error("Invalid hamiltonian matrix detected, stopping execution.")
end


# Define a function to apply the exponential of a matrix to a state vector
function apply_exp_H(state, H::Operator, t)
    H_matrix = sparse(H.data)  # Extract the sparse matrix from the Operator
    expmv(-1im * t, H_matrix, state.data)  # Pass the matrix and vector (state.data) to expmv
end

# Evolve the state using the Krylov subspace method
function krylov_fidel(U_fwd, U_bwd, state, t_max)
    F_Ξ΅_values = []

    for t in 1:t_max
        # Evolve state using U_fwd
        state_vec = apply_exp_H(state, U_fwd, 1)
        state = Ket(state.basis, state_vec)  # Recreate the state object with the new vector
        state /= safe_norm(state)  # Normalize the state

        # Calculate Loschmidt echo
        echo_state_vec = apply_exp_H(state, U_bwd, 1)
        echo_state = Ket(state.basis, echo_state_vec)
        F_Ξ΅ = abs(dagger(initial_state) * echo_state)^2
        push!(F_Ξ΅_values, F_Ξ΅)
        println("F_Ξ΅ at t=$t: ", F_Ξ΅)
    end

    return F_Ξ΅_values
end

# Assuming U_Ξ± and U_Ξ±_Ξ΅ are your Hamiltonians or precomputed unitary matrices
t_max = 10
F_Ξ΅_values = krylov_fidel(H_kick + H_precession, H_kick + (Ξ± + Ξ΅) * Jz, initial_state, t_max)


# Compute Loschmidt echo and QFI
QFI_values = [4 * (1 - F_Ξ΅) / (Ξ΅^2) for F_Ξ΅ in F_Ξ΅_values]
println("Loschmidt echo: ", F_Ξ΅_values)
println("QFI: ", QFI_values)

#Plot the QFI values
using PyPlot

plot(1:t_max, QFI_values, linewidth=2)
xlabel("Time (t)")
ylabel("QFI I_Ξ±(t)")
title("QFI using Loschmidt Echo for j = 1800 and k = 7 ")
savefig("QFI_Loschmidt_Fig.png")
println("Plot saved as QFI_Loschmidt_Fig.png.")

This issue can be reproduced by running the provided code with j = 1800 or higher. The NaN values appear during the display of the norm of the initial state and computation of the Loschmidt echo and QFI.

  • Julia version: 1.7.2
  • Packages and versions:
    • QuantumOptics: QuantumOptics v1.1.1
    • KrylovKit: KrylovKit v0.6.1
    • Expokit: Expokit v0.2.0
  • The code works fine for j <= 1700
1 Like

Hey there :slight_smile:
Please have look at

And enclose the code you posted in triple backticks ```
Also your Julia version is rather old. I’d recommend updating to the current version 1.10.5 (I don’t think this causes your problems but it’s better to use a supported version anyways).

Can you please also provide some sample output where the vector breaks down? Is it just the vector that (over time) gets NaNs or do the Hamiltonians contains NaNs from the beginning?

Try updating your KrylovKit.jl since your are 2 versions behind. Perhaps that helps already. Otherwise try increasing the Krylov dimension.

The renormalization step in your loop should in principle be uncessary. If the norm changes (more than a couple eps) during the Krylov propagation then went wrong already because it was not unitary.

Hi,

Thank you for your inputs. I have included the program within the backticks now.
The initial state generation itself has a problem it seems, the norm of it comes out as NaN.
This is the output I am getting:

Initial state norm: NaN
Hamiltonian Matrix is valid
Hamiltonian Matrix is valid
ERROR: LoadError: ArgumentError: matrix contains Infs or NaNs.
Stacktrace:
  [1] chkfinite(A::Matrix{ComplexF64})
    @ LinearAlgebra.LAPACK /opt/julias/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:97
  [2] gebal!(job::Char, A::Matrix{ComplexF64})
    @ LinearAlgebra.LAPACK /opt/julias/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:227
  [3] exp!(A::Matrix{ComplexF64})
    @ LinearAlgebra /opt/julias/julia-1.7.2/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:618
  [4] expmv!(w::Vector{ComplexF64}, t::Complex{Int64}, A::SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, vec::Vector{ComplexF64}; tol::Float64, m::Int64, norm::typeof(norm), anorm::Float64)
    @ Expokit ~/.julia/packages/Expokit/fQbKs/src/expmv.jl:129
  [5] #expmv!#2
    @ ~/.julia/packages/Expokit/fQbKs/src/expmv.jl:42 [inlined]
  [6] #expmv#1
    @ ~/.julia/packages/Expokit/fQbKs/src/expmv.jl:34 [inlined]
  [7] expmv
    @ ~/.julia/packages/Expokit/fQbKs/src/expmv.jl:33 [inlined]
  [8] apply_exp_H(state::Ket{SpinBasis{1800//1, Int64}, Vector{ComplexF64}}, H::Operator{SpinBasis{1800//1, Int64}, SpinBasis{1800//1, Int64}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, t::Int64)
    @ Main ~/julia/QFI2.jl:61
  [9] krylov_fidel(U_fwd::Operator{SpinBasis{1800//1, Int64}, SpinBasis{1800//1, Int64}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, U_bwd::Operator{SpinBasis{1800//1, Int64}, SpinBasis{1800//1, Int64}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, state::Ket{SpinBasis{1800//1, Int64}, Vector{ComplexF64}}, t_max::Int64)
    @ Main ~/julia/QFI2.jl:70
 [10] top-level scope
    @ ~/julia/QFI2.jl:87

The Load Error comes at this line-

F_Ξ΅_values = krylov_fidel(H_kick + H_precession, H_kick + (Ξ± + Ξ΅) * Jz, initial_state, t_max)

I have updated the KrylovKit pkg and tried again, but the problem persists.

So what breaks is clearly the state vector generated by coherentspinstate. You can observe that easily by:

julia> using LinearAlgebra, QuantumOptics;
julia> coherentspinstate(SpinBasis(1600), Ο€/2, Ο€/2) |> norm
1.0000000000000115

julia> coherentspinstate(SpinBasis(1700), Ο€/2, Ο€/2) |> norm
1.000000000000014

julia> coherentspinstate(SpinBasis(1750), Ο€/2, Ο€/2) |> norm
1.0000000000000078

julia> coherentspinstate(SpinBasis(1800), Ο€/2, Ο€/2) |> norm
NaN

This is ultimately a bug in QuantumOptics.jl and you should open an issue on their github.
I think it must be some kind of under-/overflow problem. However, looking at their code I don’t immediately see why.

I tried using big floats for the computation but that didn’t help interestingly enough:

julia> coherentspinstate(SpinBasis(1800), big(Ο€)/2, big(Ο€)/2) |> norm
NaN

I also don’t quite see why we get NaN… Or rather I see the NaN likely stems from a product of Inf*0 but I don’t see why there is an Inf :thinking:

EDIT: Okay I think I understood the problem. Their code is:

function coherentspinstate(b::SpinBasis, theta::Real, phi::Real)

    result = Ket(b)
    data = result.data

    N = length(b)-1
    Ξ± = sin(theta / 2) * exp(1im * phi / 2)
    Ξ² = cos(theta / 2) * exp(-1im * phi / 2)

    # forward pass: `c_n = sqrt(binomial(N, n)) * Ξ±^n` with `n  β‰₯ 0`,
    # using recursive `binomial(N, n) = ((N+1-n)/n) * binomial(N, n-1)`
    coefficient = 1.0
    @inbounds for n = 1:N+1
        data[n] = coefficient
        coefficient *= Ξ± * sqrt((N + 1 - n) / n)
    end

    # backward pass:  c_n *= Ξ²^(N-n)
    factor = 1.0
    @inbounds for n = N+1:-1:1
        data[n] *= factor
        factor *= Ξ²
    end

    return result

end

Now if N is too large, then it overflows to Inf at some point and so after the forward pass the data array contains only Inf after some point because the magnitude just gets too large. The NaNs are created just because that’s what happens in Complex arithmetic, when there are zeros and Infs due to the cross multiplication.

Now I didn’t follow your physical problem to much, but can you perhaps start in a different state? Because this works:

julia> coherentspinstate(SpinBasis(1800), 0, 0) |> norm
1.0

Ok found a simple workaround:

julia> function mycoherentspinstate(b::SpinBasis, theta::Real, phi::Real)
       
           result = Ket(b)
           data = result.data
       
           N = length(b)-1
           Ξ± = sin(theta / 2) * exp(1im * phi / 2) |> log
           Ξ² = cos(theta / 2) * exp(-1im * phi / 2) |> log
       
           # forward pass: `c_n = sqrt(binomial(N, n)) * Ξ±^n` with `n  β‰₯ 0`,
           # using recursive `binomial(N, n) = ((N+1-n)/n) * binomial(N, n-1)`
           coefficient = zero(Ξ±)
           @inbounds for n = 1:N+1
               data[n] = coefficient
               coefficient += Ξ± + 0.5*log((N + 1 - n) / n)
           end
       
           # backward pass:  c_n *= Ξ²^(N-n)
           factor = zero(Ξ²)
           @inbounds for n = N+1:-1:1
               data[n] = exp(factor+data[n])
               factor += Ξ²
           end
       
           return result
       
       end

julia> mycoherentspinstate(SpinBasis(1800), Ο€/2, Ο€/2) |> norm
1.0000000000291729

It’s a bit slower though:

julia> @benchmark mycoherentspinstate($(SpinBasis(1400)), $(Ο€/2), $(pi/2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  72.425 ΞΌs …  1.743 ms  β”Š GC (min … max): 0.00% … 91.73%
 Time  (median):     74.869 ΞΌs              β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   77.613 ΞΌs Β± 31.335 ΞΌs  β”Š GC (mean Β± Οƒ):  0.83% Β±  2.00%

   β–‚β–„β–‡β–ˆβ–ˆβ–†β–…β–„β–„β–„β–„β–ƒβ–ƒβ–β–                   ▁ ▁▁▁   ▁                β–‚
  β–…β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–‡β–‡β–‡β–‡β–…β–†β–…β–†β–†β–‡β–†β–†β–†β–†β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–ˆβ–‡β–‡β–‡β–‡β–†β–‡β–†β–†β–†β–†β–† β–ˆ
  72.4 ΞΌs      Histogram: log(frequency) by time      99.2 ΞΌs <

 Memory estimate: 43.91 KiB, allocs estimate: 3.

julia> @benchmark coherentspinstate($(SpinBasis(1400)), $(Ο€/2), $(pi/2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  21.720 ΞΌs …  1.381 ms  β”Š GC (min … max): 0.00% … 92.78%
 Time  (median):     23.537 ΞΌs              β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   24.997 ΞΌs Β± 28.233 ΞΌs  β”Š GC (mean Β± Οƒ):  2.36% Β±  2.07%

     β–‚ β–ˆβ–‚                                                      
  β–‚β–‚β–ƒβ–ˆβ–ˆβ–ˆβ–ˆβ–…β–…β–„β–ƒβ–„β–ƒβ–ƒβ–ƒβ–ƒβ–‚β–ƒβ–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–β–‚β–‚β–β–β–‚β–β–β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚ β–ƒ
  21.7 ΞΌs         Histogram: frequency by time        40.3 ΞΌs <

 Memory estimate: 43.91 KiB, allocs estimate: 3.

But seems similarily accurate (for cases it works):

julia> coherentspinstate(SpinBasis(1500), Ο€/2, pi/2) β‰ˆ mycoherentspinstate(SpinBasis(1500), Ο€/2, pi/2)
true

EDIT: Just saw that this issue arose previously (though with a bit more moderate spin length): Binomial overflow when constructing spin coherent state with large S Β· Issue #43 Β· qojulia/QuantumOpticsBase.jl Β· GitHub

1 Like

Thanks a lot!
The program works just fine now.
It’s alright even if it’s a little slow :grin:
There’s one thing I don’t understand still. I took the inner product of the state generated from the coherentspinstate() and your function for j= 1500 and they come as near orthogonal :

Inner product: -3.3306690738754696e-16 + 2.4524326690180108e-20im
Magnitude of inner product: 3.3306690829043253e-16

and small j = 5 too:

Inner product: 1.1102230246251565e-16 + 7.601180602262292e-18im
Magnitude of inner product: 1.1128220698128932e-16

I may be wrong to compare them though :sweat_smile:

That cannot be since I checked that the results are identical:

What did you do precisely?

Yeah, that’s why I was confused as to why the inner product came as near 0.
Here’s the code:

state1 = coherentspinstate(J, Ο€/2, Ο€/2).data
state2 = mycoherentspinstate(J, Ο€/2, Ο€/2).data

# Calculate the inner product
inner_product = dot(conj(state1), state2)

# Output the inner product
println("Inner product: ", inner_product)
println("Magnitude of inner product: ", abs(inner_product))

Dot applies conj to its first argument in this form. Try:

dot(state1, state2)

Oh yes,
Now its 1 :tada:

Magnitude of inner product: 1.000000000018987

Thanks!