feat: information and measurement

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
@@ -23,6 +24,7 @@ CheckConcreteStructs = "0.1"
 DataInterpolations = "8"
 DiffEqNoiseProcess = "5.27.0"
 ExplicitImports = "1.11"
+ForwardDiff = "1"
 FunctionWrappers = "1"
 JET = "0.9.18, 0.10, 0.11"
 LinearAlgebra = "1.10"

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -19,6 +20,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 [compat]
 DataInterpolations = "8"
 Documenter = "1"
+ForwardDiff = "1"
 LaTeXStrings = "1"
 MacroTools = "0.5"
 ModelingToolkitBase = "1"

docs/make.jl

@@ -30,6 +30,7 @@ pages = [
         "examples/07-2_hong-ou-mandel__quantum-pulse.md",
         "examples/08-1_pulse-delay__simple.md",
         "examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.md",
+        "examples/10-1_quantum-fisher-information__automatic-differentiation.md",
     ],
 ]
 

docs/src/api.md

@@ -94,6 +94,28 @@ solve_mode_evolution_symmetric
 correlation_matrix
 ```
 
+## [Information](@id API: Information)
+
+```@docs
+quantum_fisher_information
+```
+
+```@docs
+classical_fisher_information
+```
+
+```@docs
+povm_probabilities
+```
+
+```@docs
+projective_measurement
+```
+
+```@docs
+parameter_derivative
+```
+
 ## [Utilities](@id API: Utilities)
 
 ```@docs

examples/10-1_quantum-fisher-information__automatic-differentiation.jl

@@ -0,0 +1,205 @@
+# # Quantum Fisher Information with Automatic Differentiation
+#
+# Quantum Fisher information (QFI) answers an optimistic question: if a quantum
+# state ``\rho(\theta)`` depends on a parameter ``\theta``, how much information
+# about ``\theta`` is available in the state before we commit to a particular
+# measurement? For one copy of the state it gives the quantum Cramér-Rao bound
+#
+# ```math
+# \mathrm{Var}(\hat\theta) \geq \frac{1}{F_Q}.
+# ```
+#
+# The bound is useful only if we can also ask a more practical question: how much
+# information is obtained by the measurements we can actually perform? For a POVM
+# or projective measurement with outcomes ``M_i``,
+#
+# ```math
+# p_i(\theta) = \mathrm{Tr}[M_i \rho(\theta)],
+# \qquad
+# F_C = \sum_i \frac{(\partial_\theta p_i)^2}{p_i}.
+# ```
+#
+# The classical Fisher information (CFI) satisfies ``F_C \leq F_Q``. Comparing
+# the two tells us whether a measurement is close to optimal or whether useful
+# information remains hidden in coherences that the measurement does not access.
+#
+# In this example we estimate the detuning of a driven Kerr parametric
+# oscillator. The derivative ``\partial_\Delta \rho(t)`` is obtained by
+# differentiating the SLH time evolution with `ForwardDiff.jl`. We then compare
+# the QFI with the CFI of three different projective readouts.
+
+using QuantumInputOutput
+using SecondQuantizedAlgebra
+using QuantumOptics
+using LinearAlgebra
+using Plots
+
+#
+
+# ## SLH Model
+#
+# The Kerr parametric oscillator has Hamiltonian
+#
+# ```math
+# H = -\Delta\,a^\dagger a
+#     + K\,a^\dagger a^\dagger a a
+#     - G\,(a^\dagger a^\dagger + a a),
+# ```
+#
+# and one decay channel ``L = \sqrt{\gamma}\,a``. We keep the parameters
+# symbolic in the SLH model and provide numerical values only when we translate
+# and evolve it.
+
+h = FockSpace(:kpo)
+a = Destroy(h, :a)
+
+@variables Δ::Real K::Real G::Real γ::Real
+
+H = -Δ * a' * a + K * a' * a' * a * a - G * (a' * a' + a * a)
+Gkpo = SLH(1, √(γ) * a, H)
+nothing # hide
+
+#
+
+const K_ = 0.001
+const G_ = 0.002
+const γ_ = 0.01
+const N = 20
+
+b = FockBasis(N - 1)
+ψ0 = fockstate(b, 0)
+ρ0 = dm(ψ0)
+
+params(Δ_) = Dict(Δ => Δ_, K => K_, G => G_, γ => γ_)
+nothing # hide
+
+# ## Candidate Measurements
+#
+# `classical_fisher_information` works with any list of POVM elements. For
+# projective readout, [`projective_measurement`](@ref) builds the projectors from
+# a Hermitian observable. This keeps the package API independent of the physical
+# system: the readout below happens to use a truncated oscillator basis, but the
+# Fisher-information code only sees projectors.
+#
+# We compare:
+#
+# - occupation readout, resolving each basis state;
+# - parity readout, which intentionally coarse-grains occupation into even and
+#   odd outcomes;
+# - quadrature readout, here represented by the eigenprojectors of ``X=a+a^\dagger``
+#   in the same truncated basis.
+
+aqo = destroy(b)
+
+occupation = Operator(b, b, Diagonal(collect(0:(N - 1))))
+parity = Operator(b, b, Diagonal((-1.0) .^ collect(0:(N - 1))))
+quadrature = aqo + aqo'
+
+measurements = (
+    occupation = projective_measurement(occupation),
+    parity = projective_measurement(parity),
+    quadrature = projective_measurement(quadrature),
+)
+
+function cfi_for(measurement, ρ, dρ)
+    return classical_fisher_information(
+        ρ,
+        dρ,
+        measurement;
+        probability_floor = 1e-12,
+        derivative_floor = 1e-6,
+    )
+end
+nothing # hide
+
+# ## Differentiating the SLH Trajectory
+#
+# [`parameter_derivative`](@ref) evolves the density matrix to a chosen time and
+# returns both ``\rho(t)`` and ``\partial_\Delta\rho(t)``. The `estimate` keyword
+# names the symbolic parameter to differentiate, while `parameter` fixes the
+# numerical working point.
+
+ρ, dρ = parameter_derivative(Gkpo, b, ρ0, 2000; estimate = Δ, parameter = params(0.0))
+
+final_information = (
+    QFI = quantum_fisher_information(ρ, dρ),
+    occupation = cfi_for(measurements.occupation, ρ, dρ),
+    parity = cfi_for(measurements.parity, ρ, dρ),
+    quadrature = cfi_for(measurements.quadrature, ρ, dρ),
+)
+@info "Fisher information at final time" final_information
+nothing # hide
+
+# ## Readout Comparison
+#
+# Repeating the differentiated solve over a grid of readout times gives a direct
+# comparison between the measurement-independent QFI bound and concrete
+# measurement choices. The full occupation measurement keeps more information
+# than parity because it resolves the individual outcomes that parity combines.
+# Quadrature probes a different basis and can reveal parameter sensitivity that
+# is not visible in occupation probabilities alone.
+
+ts = range(0, 2000; length = 100)
+
+F_t = map(ts) do t
+    ρ, dρ = parameter_derivative(Gkpo, b, ρ0, t; estimate = Δ, parameter = params(0.0))
+    (
+        quantum = quantum_fisher_information(ρ, dρ),
+        occupation = cfi_for(measurements.occupation, ρ, dρ),
+        parity = cfi_for(measurements.parity, ρ, dρ),
+        quadrature = cfi_for(measurements.quadrature, ρ, dρ),
+    )
+end
+
+QFI_t = getproperty.(F_t, :quantum)
+CFI_occupation_t = getproperty.(F_t, :occupation)
+CFI_parity_t = getproperty.(F_t, :parity)
+CFI_quadrature_t = getproperty.(F_t, :quadrature)
+nothing # hide
+
+# ## Visualization
+
+p_info = plot(
+    ts,
+    QFI_t;
+    lw = 2.5,
+    label = "QFI bound",
+    xlabel = "time",
+    ylabel = "Fisher information",
+    grid = true,
+    size = (700, 400),
+)
+plot!(p_info, ts, CFI_occupation_t; lw = 2, ls = :dash, label = "occupation CFI")
+plot!(p_info, ts, CFI_parity_t; lw = 2, ls = :dot, label = "parity CFI")
+plot!(p_info, ts, CFI_quadrature_t; lw = 2, ls = :dashdot, label = "quadrature CFI")
+p_info
+
+#
+
+# At the final time we can also inspect the measurement probabilities. A readout
+# with high CFI is not merely one with a broad distribution: it must have outcome
+# probabilities that change strongly when the estimated parameter changes.
+
+occupation_probabilities = povm_probabilities(ρ, measurements.occupation)
+
+p_probs = bar(
+    0:(N - 1),
+    occupation_probabilities;
+    xlabel = "occupation outcome",
+    ylabel = "probability",
+    label = "",
+    grid = false,
+    size = (700, 300),
+)
+p_probs
+
+# ## Package versions
+
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(
+    ["QuantumInputOutput", "SecondQuantizedAlgebra", "QuantumOptics", "ForwardDiff", "Plots"],
+    mode = PKGMODE_MANIFEST,
+)

examples/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -18,6 +19,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Documenter = "1"
+ForwardDiff = "1"
 LaTeXStrings = "1"
 MacroTools = "0.5"
 ModelingToolkitBase = "1"
@@ -26,6 +28,6 @@ Plots = "1"
 QuantumCumulants = "0.5"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.5"
-SecondQuantizedAlgebra = "0.6"
+SecondQuantizedAlgebra = "0.8.2"
 SymbolicUtils = "4"
 Symbolics = "7"

src/QuantumInputOutput.jl

@@ -12,6 +12,7 @@ using LinearAlgebra: LinearAlgebra, I, mul!
 using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, Tsit5, solve
 using StaticArrays: StaticArrays, SMatrix, SVector
 using FunctionWrappers: FunctionWrappers, FunctionWrapper
+using ForwardDiff: ForwardDiff
 
 const SQA = SecondQuantizedAlgebra
 
@@ -42,6 +43,12 @@ export SLH,
     solve_mode_evolution_symmetric,
     # Correlations
     correlation_matrix,
+    # Information
+    quantum_fisher_information,
+    classical_fisher_information,
+    povm_probabilities,
+    projective_measurement,
+    parameter_derivative,
     # Operators
     substitute_operators
 
@@ -51,5 +58,6 @@ include("utils.jl")
 include("pulses.jl")
 include("correlations.jl")
 include("interaction_picture.jl")
+include("information.jl")
 
 end