A DensityMatrix encodes both pure and mixed quantum states as a matrix ρ satisfying Tr(ρ) = 1. It lives in the github.com/itsubaki/q/quantum/density package.
import (
"github.com/itsubaki/q/quantum/density"
"github.com/itsubaki/q/quantum/qubit"
)
// Pure state from a single qubit
qb := qubit.Zero()
rho := density.From(qb)
fmt.Println(rho.IsPure()) // true
fmt.Println(rho.Purity()) // 1.0
// Mixed state from a classical ensemble
states := []density.WeightedState{
{Probability: 0.5, Qubit: qubit.Zero()},
{Probability: 0.5, Qubit: qubit.One()},
}
mixed := density.FromStates(states)
fmt.Println(mixed.IsMixed()) // true
WeightedState
WeightedState pairs a qubit with a classical probability, forming one term in a mixed ensemble.
type WeightedState struct {
Probability float64
Qubit *qubit.Qubit
}
Classical probability of this state in the ensemble. Must be non-negative. The sum across all states in the ensemble must equal 1.0.
The quantum state associated with this weight. All qubits in an ensemble must have the same dimension.
Constructors
From
func From(qb *qubit.Qubit) *DensityMatrix
ρ = |ψ⟩⟨ψ| from a single qubit. Equivalent to calling FromStates with a single entry at probability 1.0.
rho := density.From(qubit.Zero())
FromStates
func FromStates(states []WeightedState) *DensityMatrix
Normalize first if they do not already sum to 1.
states := []density.WeightedState{
{Probability: 0.7, Qubit: qubit.Zero()},
{Probability: 0.3, Qubit: qubit.One()},
}
rho := density.FromStates(states)
Normalize
func Normalize(states []WeightedState) []WeightedState
WeightedState with probabilities rescaled so they sum to 1. The original slice is not modified.
raw := []density.WeightedState{
{Probability: 2, Qubit: qubit.Zero()},
{Probability: 8, Qubit: qubit.One()},
}
normalized := density.Normalize(raw) // probabilities become 0.2, 0.8
IsValid
func IsValid(states []WeightedState, tol ...float64) bool
WeightedState is suitable for constructing a density matrix. Returns false if any of the following hold:
- The slice is empty.
- Qubits have inconsistent dimensions.
- Any probability is negative.
- Probabilities do not sum to 1 (within optional tolerance
tol).
ok := density.IsValid(states) // default tolerance
ok = density.IsValid(states, 1e-10) // custom tolerance
Properties
func (m *DensityMatrix) At(i, j int) complex128
ρᵢⱼ at row i, column j (0-based).
Matrix
func (m *DensityMatrix) Matrix() *matrix.Matrix
*matrix.Matrix. Useful when you need to pass the raw matrix to other parts of the library (e.g., gate.TensorProduct).
Dim
func (m *DensityMatrix) Dim() (rows, cols int)
n-qubit system the matrix is 2ⁿ × 2ⁿ.
NumQubits
func (m *DensityMatrix) NumQubits() int
n such that Dim() == (2ⁿ, 2ⁿ).
Trace
func (m *DensityMatrix) Trace() float64
Tr(ρ). For a valid density matrix this is always 1.0.
Purity
func (m *DensityMatrix) Purity() float64
Tr(ρ²). A value of 1 indicates a pure state; values less than 1 indicate a mixed state. For an n-qubit maximally mixed state the purity equals 1 / 2ⁿ.
IsPure
func (m *DensityMatrix) IsPure(tol ...float64) bool
Purity() ≈ 1 within optional tolerance tol.
IsMixed
func (m *DensityMatrix) IsMixed(tol ...float64) bool
IsPure() is false — i.e., whether the state is a non-trivial mixture.
IsHermitian
func (m *DensityMatrix) IsHermitian(tol ...float64) bool
ρ = ρ† within optional tolerance tol. Every valid density matrix is Hermitian.
Probability
func (m *DensityMatrix) Probability(q *qubit.Qubit) float64
Tr(ρ |q⟩⟨q|) — the probability that a measurement in the computational basis yields the state q.
p := rho.Probability(qubit.Zero()) // prob of measuring |0⟩
ExpectedValue
func (m *DensityMatrix) ExpectedValue(u *matrix.Matrix) float64
Tr(ρ U) — the expectation value of the observable U in state ρ.
ev := rho.ExpectedValue(gate.Z()) // ⟨Z⟩ in state ρ
Seq2
func (m *DensityMatrix) Seq2() iter.Seq2[int, []complex128]
complex128 values.
for i, row := range rho.Seq2() {
fmt.Println(i, row)
}
Operations
Apply
func (m *DensityMatrix) Apply(u *matrix.Matrix) *DensityMatrix
ρ' = U ρ U†. Returns a new DensityMatrix; the receiver is not modified.
rho2 := rho.Apply(gate.H()) // apply Hadamard (single-qubit system)
ApplyChannel
func (m *DensityMatrix) ApplyChannel(channel ...*Channel) *DensityMatrix
*Channel values sequentially. Each channel is expanded from its Kraus operators against the current number of qubits. Returns a new DensityMatrix.
ApplyChannelFunc
func (m *DensityMatrix) ApplyChannelFunc(channel ...ChannelFunc) *DensityMatrix
ChannelFunc values. Each function receives the qubit count n and returns a *Channel, which is then applied via ApplyChannel. This is the preferred high-level API for the built-in noise channels.
noisy := rho.ApplyChannelFunc(
density.BitFlip(0.1, 0),
density.PhaseFlip(0.05, 0),
)
ApplyKraus
func (m *DensityMatrix) ApplyKraus(ops ...*matrix.Matrix) *DensityMatrix
ρ' = Σᵢ Kᵢ ρ Kᵢ†. Returns a new DensityMatrix. Use this when you have hand-crafted operator matrices.
// Custom Kraus operators
rho2 := rho.ApplyKraus(k0, k1)
PartialTrace / TraceOut
func (m *DensityMatrix) PartialTrace(qb ...int) *DensityMatrix
func (m *DensityMatrix) TraceOut(qb ...int) *DensityMatrix
TraceOut is an alias for PartialTrace.
The length of qb must be at most NumQubits() - 1; you cannot trace out all qubits.
// 2-qubit Bell state — trace out qubit 1 to get the reduced state of qubit 0
bell := density.From(qubit.From("00").H(0).CX(0, 1))
rho0 := bell.PartialTrace(1) // reduced state for qubit 0
Qubit indices are 0-based, ordered from most significant to least significant bit, matching the rest of the itsubaki/q API.
TensorProduct
func (m *DensityMatrix) TensorProduct(n *DensityMatrix) *DensityMatrix
ρ ⊗ σ — the tensor product of two density matrices, representing a composite system with no entanglement between the two subsystems.
rhoAB := rhoA.TensorProduct(rhoB)
Project
func (m *DensityMatrix) Project(q *qubit.Qubit, tol ...float64) (float64, *DensityMatrix)
|q⟩. Returns:
p float64 — the probability of obtaining outcome q.*DensityMatrix — the post-measurement state (|q⟩⟨q| ρ |q⟩⟨q|) / p.
If p is zero (within tol) the returned matrix is all-zeros.
p, post := rho.Project(qubit.Zero())
fmt.Printf("prob=%.3f\n", p)
fmt.Println(post.IsPure()) // true after projection
Clone
func (m *DensityMatrix) Clone() *DensityMatrix
Convenience channel methods
DensityMatrix exposes thin wrappers around ApplyChannelFunc for every built-in channel. Each returns a new *DensityMatrix; the receiver is not modified.
| Method | Delegates to |
|---|
BitFlip(p float64, qb int) | density.BitFlip |
PhaseFlip(p float64, qb int) | density.PhaseFlip |
BitPhaseFlip(p float64, qb int) | density.BitPhaseFlip |
Depolarizing(p float64, qb int) | density.Depolarizing |
AmplitudeDamping(gamma float64, qb int) | density.AmplitudeDamping |
PhaseDamping(gamma float64, qb int) | density.PhaseDamping |
Pauli(px, py, pz float64, qb int) | density.Pauli |
Flip(p float64, u *matrix.Matrix, qb int) | density.Flip |
