API

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.

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

Construct a linear inversion method for quantum state tomography.

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.

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)

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)

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])

triangular_indices(d)

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

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

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

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.

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.

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

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

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

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

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.

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.

Simplex(dim::Int)

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

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.

GinibreEnsamble(dim::Int)

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

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.

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.

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.

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.

orthogonal_projection(ρ, set)

Calculate the orthogonal projection of ρ onto set.

set is an array with one more dimension than ρ.

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.

linear_combination(xs, set)

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

linear_combination!(ρ, xs, set)

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

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.

cond(M, p::Real=2)

Condition number of the matrix M, computed using the operator p-norm. Valid values for p are 1, 2 (default), or Inf.

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.

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.