Maximizing Mutual Information and Accessible Information with Python
Suppose Alice sends Bob one of several non-orthogonal quantum states. Bob wants to learn as much as possible about which one she sent. Which measurement should he perform?
This is the problem of accessible information. In this article we solve it numerically for a classic example, the trine ensemble, and check the result against the known analytic answer. Along the way we look at the geometry of the optimal measurement on the Bloch sphere.
1. Problem Setting
An ensemble is a set of states $\rho_i$ with prior probabilities $p_i$:
$$
\mathcal{E}={p_i,\rho_i}_{i=1}^{N}.
$$
Bob measures with a POVM ${E_k}$, where
$$
E_k\ge 0,\qquad \sum_k E_k=\mathbb{1}.
$$
The conditional probability of outcome $k$ given state $i$ is
$$
p(k|i)=\mathrm{Tr}(\rho_i E_k),
$$
and the outcome distribution is
$$
q_k=\sum_i p_i,p(k|i).
$$
The mutual information between the sent index $i$ and the outcome $k$ is
$$
I(X!:!Y)=\sum_{i,k}p_i,p(k|i)\log_2\frac{p(k|i)}{q_k}.
$$
The accessible information is the maximum over all POVMs:
$$
I_{\mathrm{acc}}=\max_I(X!:!Y).
$$
It is bounded above by the Holevo quantity:
$$
I_{\mathrm{acc}}\le\chi=S!\left(\sum_i p_i\rho_i\right)-\sum_i p_i S(\rho_i),
$$
where $S(\rho)=-\mathrm{Tr},\rho\log_2\rho$ is the von Neumann entropy.
Two facts shape the numerical strategy:
- By Davies’ theorem, an optimal POVM can always be chosen with rank-one elements, and at most $d^2$ of them. For a qubit ($d=2$), four outcomes are enough.
- The objective is non-convex in the POVM, so a single local optimization can get stuck. We use many random restarts.
2. The Example: the Trine Ensemble
The trine ensemble consists of three pure qubit states whose Bloch vectors lie in the $x$–$z$ plane, $120^\circ$ apart, each with probability $1/3$:
$$
\vec r_i=\bigl(\sin\phi_i,;0,;\cos\phi_i\bigr),\qquad \phi_i=\frac{2\pi i}{3},\quad i=0,1,2,
$$
$$
\rho_i=\frac12\bigl(\mathbb{1}+\vec r_i\cdot\vec\sigma\bigr).
$$
Because $\sum_i\vec r_i=0$, the average state is maximally mixed, $\bar\rho=\mathbb{1}/2$. Since the states are pure,
$$
\chi=S(\bar\rho)-0=1\ \text{bit}.
$$
Three non-orthogonal states cannot be perfectly distinguished, so the accessible information must be strictly smaller than $1$ bit.
A closed-form candidate
Consider the anti-trine POVM
$$
E_k=\frac23\cdot\frac12\bigl(\mathbb{1}-\vec r_k\cdot\vec\sigma\bigr),\qquad k=0,1,2.
$$
Each element projects onto the state orthogonal to $\rho_k$, and $\sum_k E_k=\mathbb{1}$ because $\sum_k\vec r_k=0$. The conditional probabilities are
$$
p(k|i)=\frac13\bigl(1-\vec r_i\cdot\vec r_k\bigr)=
\begin{cases}
0 & k=i,\[2pt]
\tfrac12 & k\neq i,
\end{cases}
$$
so $q_k=1/3$ and
$$
I=H(Y)-H(Y|X)=\log_2 3-1=\log_2\frac32\approx 0.585\ \text{bits}.
$$
This is the known value of the accessible information of the trine. In the rest of the article we let the computer rediscover it, without telling the optimizer anything about the structure.
3. Strategy for the Numerical Optimization
Parametrizing a valid POVM
Optimizers work best without constraints, so we build a POVM that is valid by construction. Take $K$ arbitrary complex vectors $|v_k\rangle\in\mathbb{C}^2$, form
$$
S=\sum_k|v_k\rangle\langle v_k|,
$$
and define
$$
E_k=S^{-1/2},|v_k\rangle\langle v_k|,S^{-1/2}.
$$
Then $E_k\ge0$ and $\sum_kE_k=S^{-1/2}SS^{-1/2}=\mathbb{1}$ automatically. The $4K$ real numbers (real and imaginary parts of the vectors) can be optimized freely with BFGS.
Bloch representation of a POVM element
For a qubit, every rank-one element can be written
$$
E_k=a_k\bigl(\mathbb{1}+\vec m_k\cdot\vec\sigma\bigr),\qquad |\vec m_k|=1,
$$
with the completeness conditions
$$
\sum_k a_k=1,\qquad \sum_k a_k\vec m_k=\vec 0.
$$
The pair $(a_k,\vec m_k)$ is easy to visualize: $\vec m_k$ is a direction on the Bloch sphere and $a_k$ is its weight.
A two-parameter family for a 3D landscape
To see the optimum as a peak, we define a rotated and “softened” anti-trine family:
$$
E_k(\alpha,\eta)=\frac23\cdot\frac12\bigl(\mathbb{1}+\eta,R_y(\alpha)(-\vec r_k)\cdot\vec\sigma\bigr),\qquad 0\le\eta\le1.
$$
Here $\alpha$ rotates the measurement within the $x$–$z$ plane, and $\eta$ controls the sharpness ($\eta=0$ is a trivial measurement, $\eta=1$ is rank-one). Every member is a valid POVM, so $I(\alpha,\eta)$ is a legitimate surface to plot.
4. The Complete Source Code
1 | import numpy as np |
5. Code Walkthrough
Section 1: Ensemble construction
bloch_to_rho converts a Bloch vector $\vec r$ into $\rho=\tfrac12(\mathbb 1+\vec r\cdot\vec\sigma)$ with a single tensordot. The three trine vectors are generated from the angles $0,,2\pi/3,,4\pi/3$, and rhos is an array of shape (3, 2, 2). entropy_bits computes the von Neumann entropy from the eigenvalues, discarding values below $10^{-12}$ so that $0\log 0$ never produces nan. The Holevo quantity $\chi$ is computed directly from its definition and serves as the upper bound.
Section 2: Mutual information
mutual_information is the heart of the program. The line
1 | cond = np.real(np.einsum("kab,iba->ik", povm, rhos)) |
computes all conditional probabilities $p(k|i)=\mathrm{Tr}(\rho_iE_k)$ in one shot. The index pattern kab,iba->ik contracts $E_k^{ab}\rho_i^{ba}$, which is exactly the trace of a matrix product. The remaining lines compute $q_k$ and the sum $\sum p_i,p(k|i)\log_2\frac{p(k|i)}{q_k}$. The np.where guard enforces the convention $0\log 0=0$, so zero-probability entries (which occur at the optimum!) do not generate warnings or nan.
Section 3: POVM parametrization and optimization
params_to_povm implements $E_k=S^{-1/2}|v_k\rangle\langle v_k|S^{-1/2}$. The matrix $S^{-1/2}$ is built from the eigendecomposition of the Hermitian matrix $S$, with eigenvalues floored at $10^{-12}$ for numerical safety. Because completeness holds by construction, the optimizer needs no constraints and BFGS can be used directly.
optimize_povm runs BFGS from n_restarts random starting points, records the objective after every iteration through the callback, and keeps the best run. Restarts matter here: the landscape has stationary points (for example one corresponding to projective measurements), and some runs may settle there instead of at the global optimum.
We run it for $K=3$ and $K=4$ outcomes. Since Davies’ theorem says four elements suffice for a qubit, the $K=4$ run is the general case, and $K=3$ confirms that the structure of the optimum is simpler than the upper limit suggests.
Section 4: Merging duplicate outcomes
When $K=4$, the optimizer may split one physical outcome into two elements pointing in the same direction. Merging them does not change the mutual information (coarse-graining outcomes with identical likelihood ratios loses nothing). povm_to_bloch extracts the weight $a_k=\mathrm{Tr}E_k/2$ and the direction $\vec m_k$ from each element, and merge_povm groups elements whose directions agree within a tolerance and sums their weights. The merged POVM is rebuilt and its information is recomputed as a consistency check.
Section 5: Reference measurements and the landscapes
Three baselines are computed:
- Anti-trine POVM: the analytic candidate from Section 2 above.
- Trine-aligned POVM: the same construction with $+\vec r_k$ instead of $-\vec r_k$, a natural but inferior guess.
- Best projective measurement: for a two-outcome projective measurement along axis $\hat n$, the likelihoods are $p(\pm|i)=\tfrac12(1\pm\hat n\cdot\vec r_i)$. This is evaluated on a $61\times121$ grid over the whole sphere.
Speed-up note. The grid evaluation is fully vectorized with NumPy broadcasting; the only Python-level loop runs over the three states, not over the 7,381 grid points. The POVM objective uses einsum instead of nested loops, and the unconstrained parametrization lets us use a fast quasi-Newton method instead of a constrained solver. As a result, the whole program, including 40 optimization runs, a $41\times121$ surface, and 2,000 random POVMs, finishes within a few seconds to tens of seconds.
The family rotated_povm(alpha, eta) rotates the anti-trine directions by $\alpha$ in the $x$–$z$ plane and shrinks them by $\eta$, which gives the surface $I(\alpha,\eta)$. Finally, 2,000 random 3-outcome POVMs provide a “no optimization” baseline.
Section 6: Console output
All numerical results are printed: the Holevo bound, the theoretical value $\log_2(3/2)$, the optimizer results, the reference values, the merged POVM in Bloch form, a completeness check $\sum_kE_k=\mathbb 1$, and the matrix of conditional probabilities $p(k|i)$.
Section 7: The figure
All six panels are drawn into a single fig and displayed once with plt.show(). Panels (a), (b), and (c) are 3D.
6. Execution Results
Figure output
Console output
=== Ensemble: trine states (equal priors) === Holevo bound chi : 1.000000 bits Theoretical optimum log2(3/2) : 0.584963 bits === Numerical optimization === 3-outcome POVM (best of 20) : 0.584963 bits 4-outcome POVM (best of 20) : 0.584963 bits 4-outcome POVM after merging : 0.584963 bits (3 distinct outcomes) Spread over restarts (3 outcomes) : min 0.459148 / max 0.584963 === Reference measurements === Best projective measurement : 0.459148 bits Trine-aligned POVM : 0.333333 bits Anti-trine POVM : 0.584963 bits Best of 2000 random 3-outcome POVMs : 0.532563 bits === Optimal POVM (merged): weights and Bloch directions === E_0: weight a = 0.3333, direction m = [-0. -0. -1.] E_1: weight a = 0.3333, direction m = [-0.866 0. 0.5 ] E_2: weight a = 0.3333, direction m = [0.866 0. 0.5 ] Sum of weights = 1.000000 === Completeness check === Sum of E_k = [[1.+0.j 0.+0.j] [0.-0.j 1.+0.j]] === Conditional probabilities p(k|i) of the optimal POVM === [[0. 0.5 0.5] [0.5 0. 0.5] [0.5 0.5 0. ]]
7. Interpreting the Results
Console output
The optimizer should report
$$
I_{\mathrm{acc}}\approx 0.584963\ \text{bits}=\log_2\frac32,
$$
for both the 3-outcome and the 4-outcome searches. The numerical optimum therefore reproduces the analytic value, even though the optimizer started from random vectors and knew nothing about the anti-trine structure.
Several facts stand out:
- The gap to the Holevo bound is large. The Holevo quantity is $\chi=1$ bit, but only about $0.585$ bits are accessible. The Holevo bound is not achievable for this ensemble, because the states are non-orthogonal and no measurement can fully resolve them.
- The 4-outcome optimum collapses to 3 outcomes. After merging, only three distinct directions remain. Allowing the maximum number of outcomes allowed by Davies’ theorem brings no extra information here.
- The optimal conditional probabilities are $0$ and $\tfrac12$. Each outcome rules out exactly one state: outcome $k$ never occurs when state $k$ was sent, and otherwise occurs with probability $1/2$. This matrix is $p(k|i)=\tfrac12(1-\delta_{ik})$ up to a permutation of outcomes.
- The optimal weights are $a_k=1/3$ with directions equal to the anti-trine vectors $-\vec r_k$ (in some order), so $\sum_ka_k=1$ and $\sum_ka_k\vec m_k=\vec0$.
- Restarts matter. The spread over restarts shows that some runs end below the optimum, typically at the plateau of the best projective measurement (about $0.459$ bits). Without restarts, a single unlucky run could have been mistaken for the answer.
Panel (a): projective measurements over the Bloch sphere
The color shows the mutual information obtained by a two-outcome projective measurement along each axis $\hat n$. The map is symmetric under $\hat n\to-\hat n$ and has a pattern with the three-fold symmetry of the trine. The bright regions lie in the $x$–$z$ plane near the signal directions, and the dark regions lie near the $\pm y$ poles: measuring along $y$ is orthogonal to every signal Bloch vector, so it returns pure noise and $I=0$. The maximum over all projective measurements, marked by the red dot, is only about $0.459$ bits, clearly below $\log_2(3/2)$. Projective measurements are not enough: a genuinely generalized measurement beats the best von Neumann measurement.
Panel (b): signal states and the optimal POVM
The blue arrows are the three trine states and the red arrows are the directions of the optimal POVM elements. All six vectors lie in the $x$–$z$ plane, and each red arrow points opposite to one blue arrow. Geometrically, the optimal measurement asks, for each state, the question “is the state not this one?”. The three red vectors form an equilateral triangle, so $\sum_ka_k\vec m_k=0$ is satisfied and the elements add up to the identity.
Panel (c): the 3D landscape $I(\alpha,\eta)$
This surface is the most informative picture of the optimization problem.
- At $\eta=0$ the surface is flat at $I=0$: a POVM with $E_k=\mathbb 1/3$ ignores the state entirely.
- Along $\eta$, information grows monotonically; sharper measurements always help in this family.
- Along $\alpha$, the surface oscillates with period $2\pi/3$, reflecting the three-fold symmetry of the ensemble.
- At $\eta=1$ the peaks (value $\log_2(3/2)\approx0.585$) sit at $\alpha=0,,2\pi/3,,4\pi/3$, where the measurement is exactly the anti-trine.
- Midway between the peaks, at $\alpha=\pi/3$ and its translates, the rotated anti-trine coincides with the trine-aligned POVM. There $I=1/3$ bit, the bottom of the ridge. The same physical hardware, rotated by $60^\circ$ in the Bloch sphere, loses almost half of the information.
So a sharp, correctly oriented measurement is a global maximum within this family, and the landscape shows how quickly the information degrades when the measurement is misaligned.
Panel (d): bar chart
The five bars give the ranking: Holevo bound ($1$) $>$ optimal POVM $=$ anti-trine ($\approx0.585$) $>$ best projective ($\approx0.459$) $>$ trine-aligned ($1/3$). The optimal POVM and the anti-trine bars coincide, confirming the analytic solution. The distance between the gray Holevo bar and the red optimum bar is the part of the information that is encoded in the ensemble but cannot be extracted by any measurement.
Panel (e): convergence
Both curves start from a low value determined by the random initial vectors and climb to the dashed line at $\log_2(3/2)$ within a few dozen BFGS iterations. The 4-outcome run converges to the same level as the 3-outcome run, matching the merging result: the extra degree of freedom is not used.
Panel (f): random POVMs
The histogram of 2,000 random 3-outcome POVMs lies entirely to the left of the red line. Random POVMs typically capture a fraction of the available information, and even the best of 2,000 draws stays below the optimum. This is a sanity check that the optimizer is doing real work, and that the value $\log_2(3/2)$ is a ceiling and not just one more point in the distribution.
8. Conclusion
We computed the accessible information of the trine ensemble by optimizing over POVMs, using a parametrization that enforces completeness automatically. The optimum is
$$
I_{\mathrm{acc}}=\log_2\frac32\approx0.585\ \text{bits},
$$
attained by the three-outcome anti-trine POVM
$$
E_k=\frac23\cdot\frac12\bigl(\mathbb 1-\vec r_k\cdot\vec\sigma\bigr),
$$
while the best projective measurement gives only about $0.459$ bits and the Holevo bound of $1$ bit is out of reach. The same pipeline (rank-one POVM parametrization, multi-start BFGS, and Bloch-sphere visualization) works for any qubit ensemble: change bloch_states and priors and the program finds the information-maximizing measurement for you. For higher dimensions, the parametrization generalizes directly by using $d$-dimensional complex vectors and up to $d^2$ outcomes.
. The dark-to-bright scale makes solver-level errors visible.
first and are then pushed further by the logarithmic growth of $x^{*}$. Expensive resources appear mainly at the ridge peaks, and cheap resources also cover the quiet hours.
