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

**URL:** <https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923>\
**Category:** Quantum\
**Tags:** question, package\
**Created:** [September 2, 2024, 3:51am UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923 "2024-09-02T03:51:24Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Dawn](https://avatars.discourse-cdn.com/v4/letter/d/71e660/32.png) [@Dawn](https://discourse.julialang.org/u/Dawn)\
**Post date:** [September 2, 2024, 3:51am UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/1 "2024-09-02T03:51:24Z")

</div>

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:

```julia
#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

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 2, 2024, 4:53am UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/2 "2024-09-02T04:53:48Z")

</div>

Hey there 🙂  
Please have look at

> [@Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

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.

---

<div class="post-metadata">

**Author:** ![Dawn](https://avatars.discourse-cdn.com/v4/letter/d/71e660/32.png) [@Dawn](https://discourse.julialang.org/u/Dawn)\
**Post date:** [September 2, 2024, 10:37am UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/3 "2024-09-02T10:37:17Z")

</div>

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:

```julia
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-

```julia
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.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 2, 2024, 11:28am UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/4 "2024-09-02T11:28:29Z")

</div>

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

```julia-repl
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](https://github.com/qojulia/QuantumOptics.jl).  
I think it must be some kind of under-/overflow problem. However, looking at [their code](https://github.com/qojulia/QuantumOptics.jl/blob/v1.1.1/src/phasespace.jl#L290-L312) I don’t immediately see why.

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

```julia
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` 🤔

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

```julia
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 `NaN`s are created just because that’s what happens in Complex arithmetic, when there are zeros and `Inf`s 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-repl
julia> coherentspinstate(SpinBasis(1800), 0, 0) |> norm
1.0

```

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 2, 2024, 12:00pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/5 "2024-09-02T12:00:03Z")

</div>

Ok found a simple workaround:

```julia-repl
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-repl
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
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](https://github.com/qojulia/QuantumOpticsBase.jl/issues/43)

---

<div class="post-metadata">

**Author:** ![Dawn](https://avatars.discourse-cdn.com/v4/letter/d/71e660/32.png) [@Dawn](https://discourse.julialang.org/u/Dawn)\
**Post date:** [September 3, 2024, 1:10pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/6 "2024-09-03T13:10:03Z")

</div>

Thanks a lot!  
The program works just fine now.  
It’s alright even if it’s a little slow 😁  
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 :

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

```

and small j = 5 too:

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

```

I may be wrong to compare them though 😅

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 3, 2024, 1:34pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/7 "2024-09-03T13:34:32Z")

</div>

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

> [@abraemer](#):
>
> ```julia
> julia> coherentspinstate(SpinBasis(1500), π/2, pi/2) ≈ mycoherentspinstate(SpinBasis(1500), π/2, pi/2)
> true
> 
> ```

What did you do precisely?

---

<div class="post-metadata">

**Author:** ![Dawn](https://avatars.discourse-cdn.com/v4/letter/d/71e660/32.png) [@Dawn](https://discourse.julialang.org/u/Dawn)\
**Post date:** [September 3, 2024, 2:34pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/8 "2024-09-03T14:34:28Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 3, 2024, 2:35pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/9 "2024-09-03T14:35:54Z")

</div>

> [@Dawn](#):
>
> `inner_product = dot(conj(state1), state2)`

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

```julia
dot(state1, state2)

```

---

<div class="post-metadata">

**Author:** ![Dawn](https://avatars.discourse-cdn.com/v4/letter/d/71e660/32.png) [@Dawn](https://discourse.julialang.org/u/Dawn)\
**Post date:** [September 3, 2024, 2:38pm UTC](https://discourse.julialang.org/t/nan-values-for-j-1700-in-quantum-fisher-information-calculation-using-quantumoptics/118923/10 "2024-09-03T14:38:00Z")

</div>

Oh yes,  
Now its 1 🎉

```julia
Magnitude of inner product: 1.000000000018987

```

Thanks!