BayesianInference(povm::AbstractArray{Matrix{T}},
    basis=gell_mann_matrices(size(first(povm), 1), complex(T))) where {T}

Create a Bayesian inference object from a POVM.

This is passed to the prediction method in order to perform the Bayesian inference.

GinibreEnsamble(dim::Int)

A type representing a Ginibre ensamble of complex matrices of dimension dim.

HaarUnitary(dim::Int)

A type representing a Haar-random unitary matrix of dimension dim.

HaarVector(dim::Int)

A type representing a Haar-random unitary vector of dimension dim.

History{T<:Integer}

A type that represents a history of outcomes.

Fields

• history::Vector{T}: A vector of outcomes. history[i] is the outcome of the i-th measurement.

LinearInversion(povm, basis=gell_mann_matrices(size(first(povm), 1)))

Construct a linear inversion method for quantum state tomography.

ProductMeasure(dim::Int)

A type representing a measure on the density states. It is a product Haar measure on the unitary group and a uniform (Lebesgue) measure on the simplex.

Simplex(dim::Int)

A type representing a random point on the simplex embeded in a space of dimension dim.

X_matrix(j, k, d, ::Type{T}=ComplexF32) where {T<:Union{Real,Complex}}

Compute the real off diagonal matrix of the generalized Gell-Mann matrices in dimension d.

The type of the matrix elements is T, which defaults to ComplexF32. The only non-zero elements are X[j, k] = 1 and X[k, j] = 1. The matrices are normalized to have unit Hilbert-Schmidt norm.

Examples

julia> X_matrix(1,2,2)
2×2 Matrix{ComplexF32}:
      0.0+0.0im  0.707107+0.0im
 0.707107+0.0im       0.0+0.0im

Y_matrix(j, k, d, ::Type{T}=ComplexF32) where {T<:Complex}

Compute the imaginary off diagonal matrix of the generalized Gell-Mann matrices in dimension d.

The type of the matrix elements is T, which defaults to ComplexF32. The only non-zero elements are Y[j, k] = im and Y[k, j] = -im. The matrices are normalized to have unit Hilbert-Schmidt norm.

Examples

julia> Y_matrix(1,2,2)
2×2 Matrix{ComplexF32}:
 0.0+0.0im       0.0+0.707107im
 0.0-0.707107im  0.0+0.0im

Z_matrix(j, d, ::Type{T}=ComplexF32) where {T<:Union{Real,Complex}}

Compute the j'th diagonal matrix of the generalized Gell-Mann matrices in dimension d.

The type of the matrix elements is T, which defaults to ComplexF32. The matrices are normalized to have unit Hilbert-Schmidt norm. The identity matrix is returned when j == 0.

Examples

julia> Z_matrix(0, 3)
3×3 Matrix{ComplexF32}:
 0.57735+0.0im      0.0+0.0im      0.0+0.0im
     0.0+0.0im  0.57735+0.0im      0.0+0.0im
     0.0+0.0im      0.0+0.0im  0.57735+0.0im

julia> Z_matrix(1, 3)
3×3 Matrix{ComplexF32}:
 0.707107+0.0im        0.0+0.0im  0.0+0.0im
      0.0+0.0im  -0.707107+0.0im  0.0+0.0im
      0.0+0.0im        0.0+0.0im  0.0+0.0im

julia> Z_matrix(2, 3)
3×3 Matrix{ComplexF32}:
 0.408248+0.0im       0.0+0.0im        0.0+0.0im
      0.0+0.0im  0.408248+0.0im        0.0+0.0im
      0.0+0.0im       0.0+0.0im  -0.816497+0.0im

acceptance!(x₀, x, ℓπ₀, ∇ℓπ₀, ∇ℓπ, f, σ)

Accept or reject the proposed state x given the current state x₀. If accepted, the state x₀ is updated to x and the gradient ∇ℓπ₀ is updated to ∇ℓπ. Returns a tuple with the updated log-likelihood ℓπ and a boolean indicating if the state was accepted.

augment_povm(povm::AbstractArray{Matrix{T}}, unitaries...; 
    weights=fill(one(T) / (length(unitaries) + 1), length(unitaries) + 1) where {T}

Augment a POVM (Positive Operator-Valued Measure) by applying a set of unitary transformations to it.

Arguments

• povm: The POVM to be augmented. It is represented as an array of matrices.
• unitaries: Variable number of unitary matrices representing the transformations to be applied.
• weights: An optional array of weights associated with each unitary transformation. If not provided, it defaults to a uniform distribution.

Returns

• A new POVM that is the result of applying the unitary transformations to the input POVM.

Example

bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
quater_wave_plate = [1 im; im 1] / √2
povm = augment_povm(bs_povm, half_wave_plate, quater_wave_plate, weights=[1 / 2, 1 / 4, 1 / 4])

basis_decomposition(Ω, basis=gell_mann_matrices(d))

Decompose the array Ω in the provided orthonormal basis.

If no basis is provided, the Gell-Mann matrices of appropriate dimension are used.

If Ω has dimension d, then basis should be an array with dimesnion d+1 with the last dimension indexing the basis elements.

complete_representation(history::History{T}, size) where {T}

Create a complete representation of the given history.

Arguments

• history::History: A History object which contains a history of outcomes.
• size: The size of the resulting representation.

Returns

• result: An array of size size where the i-th element is the number of times the i-th event occurred in the history.

Example

julia> h = History([1,1,1,2,1])
History{Int64}([1, 1, 1, 2, 1])

julia> complete_representation(h,(2,2))
2×2 Matrix{Int64}:
 4  0
 1  0

complete_representation(outcomes::Matrix{T}, size) where {T}

Create a complete representation of the given outcomes.

outcomes is a matrix where the first row contains the indices of non-zero elements from of the complete representation and the second row contains the corresponding non-zero values.

Returns a vector of size size where the i-th element is the value of the pair whose index is i in outcomes. If there is no such pair, the value is 0.

This function has an inverse reduced_representation.

Example

julia> outcomes = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> complete_representation(outcomes, (2,2))
2×2 Matrix{Int64}:
 3  0
 4  0

compose_povm(povms::AbstractArray{Matrix{T}}...; weights=fill(one(T) / length(povms), length(povms))) where {T}

Compose a POVM (Positive Operator-Valued Measure) from a set of given POVMs.

Arguments

• povms: Variable number of POVMs. Each POVM is represented as an array of matrices.
• weights: An optional array of weights associated with each POVM. If not provided, it defaults to a uniform distribution.

Returns

• A new POVM that is a composition of the input POVMs, weighted by their respective weights.

Example

povm1 = [rand(2,2) for _ in 1:3]
povm2 = [rand(2,2) for _ in 1:3]
composed_povm = compose_povm(povm1, povm2)

fidelity(ρ::AbstractMatrix, σ::AbstractMatrix)
fidelity(ψ::AbstractVector, φ::AbstractVector)

Calculate the fidelity between two quantum states.

The states can be pure or mixed, and they are represented by vectors ψ and φ or density matrices ρ and σ, respectively.

gell_mann_matrices(d, ::Type{T}=ComplexF32; include_identity=true) where {T<:Complex}

Generate a set of Gell-Mann matrices of dimension d.

The Gell-Mann matrices are a set of d^2 - 1 linearly independent, traceless, Hermitian matrices that, when augmented with the identity, form a basis for the space of d × d hermitian matrices.

The matrix order is real off-diagonal (X_matrix), imaginary off-diagonal (Y_matrix) and diagonal (Z_matrix). The off-diagonal matrices follow the order given by triangular_indices.

Arguments

• d: The dimension of the Gell-Mann matrices.
• include_identity: A boolean flag indicating whether to include the identity matrix in the set. If this is true, the identity is the first element of the basis

Returns

• A 3D array of Gell-Mann matrices. The last dimension is the index of the matrix in the basis.

Examples

julia> gell_mann_matrices(2,include_identity=false)
2×2×3 Array{ComplexF32, 3}:
[:, :, 1] =
      0.0+0.0im  0.707107+0.0im
 0.707107+0.0im       0.0+0.0im

[:, :, 2] =
 0.0+0.0im       0.0-0.707107im
 0.0+0.707107im  0.0+0.0im

[:, :, 3] =
 0.707107+0.0im        0.0+0.0im
      0.0+0.0im  -0.707107+0.0im

linear_combination!(ρ, xs, set)

Calculate the linear combination of the elements of set with the coefficients xs and store the result in ρ.

linear_combination(xs, set)

Calculate the linear combination of the elements of set with the coefficients xs.

log_likelihood(outcomes, povm, x, ∇ℓπ, cache1, cache2)

Returns the log-likelihood of the outcomes given the povm and the state x. The gradient of the log-likelihood is stored in ∇ℓπ.

maximally_mixed_state(d, ::Type{T}) where {T}

Returns the maximally mixed state of dimension d, represented as a vector of projections in the generalized Gell-Mann basis.

The maximally mixed state is defined as ρ = I / d.

Se also gell_mann_matrices.

orthogonal_projection(ρ, set)

Calculate the orthogonal projection of ρ onto set.

set is an array with one more dimension than ρ.

prediction(outcomes, method::LinearInversion)

Predict the quantum state from the outcomes of a tomography experiment using the LinearInversion method.

prediction(outcomes, method::BayesianInference{T};
    verbose=false,
    σ=T(1e-2),
    log_prior=x -> zero(T),
    x₀=maximally_mixed_state(Int(√size(method.povm, 2)), T),
    nsamples=10^4,
    nwarm=10^3,
    chain=nothing) where {T}

Perform a Bayesian inference on the given outcomes using the BayesianInference method.

Arguments

• outcomes: The outcomes of the experiment.
• method::BayesianInference{T}: The Bayesian inference method.
• verbose=false: Print information about the run.
• σ=T(1e-2): The initial standard deviation of the proposal distribution.
• log_prior=x -> zero(T): The log-prior function.
• x₀=maximally_mixed_state(Int(√size(method.povm, 2)), T): The initial state of the chain.
• nsamples=10^4: The number of samples to take.
• nwarm=10^3: The number of warm-up samples to take.
• chain=nothing: If not nothing, store the chain in this matrix.

Returns

A tuple with the mean state, its projection in method.basis and the covariance matrix. The mean state is already returned in matrix form.

project2density(ρ)

Project a Hermitian matrix ρ to a density matrix by setting the negative eigenvalues to zero and normalizing the trace to 1.

project2pure(ρ)

Project a Hermitian matrix ρ to a pure state by returning the eigenvector corresponding to the largest eigenvalue.

function proposal!(x, x₀, ∇ℓπ₀, σ)

Propose a new state x given the current state x₀.

The proposal is done by sampling a random vector x from a normal distribution with mean x₀ + σ^2 * ∇ℓπ₀ / 2 and covariance matrix σ^2I.

proposal_ratio(x, x₀, ∇ℓπ, ∇ℓπ₀, σ)

Returns the ratio of the transition probability of x₀ given x and the x given x₀.

Used in the acceptance step of the MALA algorithm.

real_orthogonal_projection(ρ, set)

Calculate the real part of the orthogonal projection of ρ onto set.

set is an array with one more dimension than ρ.

This function is useful when the projection is expected to be real, but numerical errors may introduce small imaginary parts.

reduced_representation(povm, outcomes)

Returns a reduced representation of both the povm and the outcomes.

One determines the nonzero elements of outcomes and then selects the corresponding columns of the povm.

This function is used in the Bayesian inference to reduce the size of the problem by ignoring unobserved outcomes.

reduced_representation(history::History)

Create a reduced representation of the given history.

Arguments

• history::History: A History object which contains a history of events.

Return a matrix where each column is a pair (event, count). The event is the unique event from the history and count is the number of times the event has occurred.

Example

julia> h = History([1,1,1,2,1])
History{Int64}([1, 1, 1, 2, 1])

julia> reduced_representation(h)
2×2 Matrix{Int64}:
 2  1
 1  4

reduced_representation(outcomes::Array{T,N}) where {T,N}

Converts a multi-dimensional array of outcomes into a 2D matrix in which the first row contains the indices of non-zero elements and the second row contains the corresponding non-zero values.

outcomes is a multi-dimensional array of outcomes where the outcomes[n] is the number of times the n-th outcome was observed.

The output is a matrix where the first row contains the indices of non-zero elements from the outcomes array and the second row contains the corresponding non-zero values.

This function has an inverse complete_representation.

Examples

julia> outcomes = [0, 1, 0, 2, 0, 3]
6-element Vector{Int64}:
 0
 1
 0
 2
 0
 3

julia> reduced_representation(outcomes)
2×3 Matrix{Int64}:
 2  4  6
 1  2  3

sample(type, n_samples)
sample(type)

Sample n_samples from type.

If n_samples is not provided, a single sample is returned.

Possible values for type are HaarUnitary, HaarVector, Simplex, ProductMeasure, and GinibreEnsamble.

sample_markov_chain(ℓπ, x₀::Vector{T}, nsamples, nwarm, basis;
    verbose=false,
    σ=oftype(T, 1e-2),
    target=0.574,
    minimum=1e-8,
    maximum=100) where {T<:Real}

Sample a Markov chain to sample the posterior of a quantum state tomography experiment using the MALA algorithm.

simulate_outcomes!(probs, N; atol=1e-3)

Simulate the N outcomes of a probability specified by the probs array. The results are stored in the probs array.

atol is the absolute tolerance for the probabilities to be considered non-negative and to sum to 1.

simulate_outcomes(ψ::AbstractVector, povm, N; atol=1e-3)
simulate_outcomes(ρ::AbstractMatrix, povm, N; atol=1e-3)
simulate_outcomes(probs, N; atol=1e-3)

Simulate the N outcomes of a quantum measurement represented by a povm on a quantum state.

The state can be pure or mixed, and it is represented by a vector ψ or a density matrix ρ, respectively. Alternativelly, one can directly provide the probabilities of the outcomes in the probs array.

atol is the absolute tolerance for the probabilities to be considered non-negative and to sum to 1.

step!(x₀, x, ℓπ₀, ∇ℓπ₀, ∇ℓπ, ℓπ_function, parameters, ρ, basis, stats, n, target, min, max)

Perform a step of the MALA algorithm.

triangular_indices(d)

Generate a vector of tuples representing the indices of the lower triangular part of a square matrix of dimension d.

unitary_transform!(povm, unitary)

Apply a unitary transformation to each element of a given POVM (Positive Operator-Valued Measure), modifing it in place.

Arguments

• povm: The POVM to be transformed. It is represented as an array of matrices.
• unitary: The unitary matrix representing the transformation to be applied.

Example

bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
unitary_transform!(bs_povm, half_wave_plate)

unitary_transform(povm, unitary)

Apply a unitary transformation to each element of a given POVM (Positive Operator-Valued Measure).

Arguments

• povm: The POVM to be transformed. It is represented as an array of matrices.
• unitary: The unitary matrix representing the transformation to be applied.

Returns

• A new POVM that is the result of applying the unitary transformation to the input POVM.

Example

bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
unitary_transform!(bs_povm, half_wave_plate)

update_σ!(parameters, n, target, min, max)

Update the parameter σ = parameters[1] of the MALA algorithm given the current iteration n and the acceptance rate parameters[2] / n. The target acceptance rate is target and the minimum and maximum values of σ are min and max, respectively.

cond(povm::Union{AbstractArray{T},AbstractMatrix{T}}, p::Real=2) where {T<:AbstractMatrix}

Calculate the condition number of the linear transformation associated with the povm.

isposdef!(ρ, xs, set)

Calculate the linear combination of the elements of set with the coefficients xs and check if the result is a positive definite matrix.