Compute the Matsumoto fidelity of two density matrices 1.
Calculate the Matsumoto fidelity between the two density matrices rho and sigma, defined by:
\[
\mathrm{tr}(\rho\#\sigma),
\]
where \(\#\) denotes the matrix geometric mean, which for invertible states is
\[
\rho\#\sigma = \rho^{1/2}\sqrt{\rho^{-1/2}\sigma\rho^{-1/2}}\rho^{1/2}.
\]
For singular states it is defined by the limit
\[
\rho\#\sigma = \lim_{\epsilon\to0}(\rho+\epsilon\mathbb{I})\#(\sigma+\epsilon\mathbb{I}).
\]
The return is a value between \(0\) and \(1\), with \(0\) corresponding to matrices rho and
sigma with orthogonal support, and \(1\) corresponding to the case rho = sigma. The Matsumoto
fidelity is a lower bound for the fidelity.
Parameters:
-
rho
(ndarray)
–
-
sigma
(ndarray)
–
Returns:
-
float | floating
–
The Matsumoto fidelity between rho and sigma.
Raises:
-
ValueError
–
If matrices are not of equal dimension.
Examples:
Consider the following Bell state
\[
u = \frac{1}{\sqrt{2}} \left( |00 \rangle + |11 \rangle \right) \in \mathcal{X}.
\]
The corresponding density matrix of \(u\) may be calculated by:
\[
\rho = u u^* = \frac{1}{2} \begin{pmatrix}
1 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 1
\end{pmatrix} \in \text{D}(\mathcal{X}).
\]
In the event where we calculate the Matsumoto fidelity between states that are identical, we should obtain
the value
of \(1\). This can be observed in |toqito⟩ as follows.
import numpy as np
from toqito.state_metrics import matsumoto_fidelity
rho = 1 / 2 * np.array(
[[1, 0, 0, 1],
[0, 0, 0, 0],
[0, 0, 0, 0],
[1, 0, 0, 1]]
)
sigma = rho
print(np.around(matsumoto_fidelity(rho, sigma), decimals=2))
References
1 Matsumoto, Keiji. Reverse test and quantum analogue of classical fidelity and generalized fidelity. (2010).
Source code in toqito/state_metrics/matsumoto_fidelity.py
| def matsumoto_fidelity(rho: np.ndarray, sigma: np.ndarray) -> float | np.floating:
r"""Compute the Matsumoto fidelity of two density matrices [@matsumoto2010reverse].
Calculate the Matsumoto fidelity between the two density matrices `rho` and `sigma`, defined by:
\[
\mathrm{tr}(\rho\#\sigma),
\]
where \(\#\) denotes the matrix geometric mean, which for invertible states is
\[
\rho\#\sigma = \rho^{1/2}\sqrt{\rho^{-1/2}\sigma\rho^{-1/2}}\rho^{1/2}.
\]
For singular states it is defined by the limit
\[
\rho\#\sigma = \lim_{\epsilon\to0}(\rho+\epsilon\mathbb{I})\#(\sigma+\epsilon\mathbb{I}).
\]
The return is a value between \(0\) and \(1\), with \(0\) corresponding to matrices `rho` and
`sigma` with orthogonal support, and \(1\) corresponding to the case `rho = sigma`. The Matsumoto
fidelity is a lower bound for the fidelity.
Args:
rho: Density operator.
sigma: Density operator.
Returns:
The Matsumoto fidelity between `rho` and `sigma`.
Raises:
ValueError: If matrices are not of equal dimension.
Examples:
Consider the following Bell state
\[
u = \frac{1}{\sqrt{2}} \left( |00 \rangle + |11 \rangle \right) \in \mathcal{X}.
\]
The corresponding density matrix of \(u\) may be calculated by:
\[
\rho = u u^* = \frac{1}{2} \begin{pmatrix}
1 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 1
\end{pmatrix} \in \text{D}(\mathcal{X}).
\]
In the event where we calculate the Matsumoto fidelity between states that are identical, we should obtain
the value
of \(1\). This can be observed in `|toqito⟩` as follows.
```python exec="1" source="above" result="text"
import numpy as np
from toqito.state_metrics import matsumoto_fidelity
rho = 1 / 2 * np.array(
[[1, 0, 0, 1],
[0, 0, 0, 0],
[0, 0, 0, 0],
[1, 0, 0, 1]]
)
sigma = rho
print(np.around(matsumoto_fidelity(rho, sigma), decimals=2))
```
"""
if not np.all(rho.shape == sigma.shape):
raise ValueError("InvalidDim: `rho` and `sigma` must be matrices of the same size.")
# If `rho` or `sigma` is a cvxpy variable then compute Matsumoto fidelity via
# semidefinite programming, so that this function can be used in the
# objective function or constraints of other cvxpy optimization problems.
if isinstance(rho, cvxpy.atoms.affine.vstack.Vstack) or isinstance(sigma, cvxpy.atoms.affine.vstack.Vstack):
w_var = cvxpy.Variable(rho.shape, hermitian=True)
objective = cvxpy.Maximize(cvxpy.real(cvxpy.trace(w_var)))
constraints = [cvxpy.bmat([[rho, w_var], [w_var, sigma]]) >> 0]
problem = cvxpy.Problem(objective, constraints)
return problem.solve()
if not is_density(rho) or not is_density(sigma):
raise ValueError("Matsumoto fidelity is only defined for density operators.")
# If `rho` or `sigma` are *not* cvxpy variables, compute Matsumoto fidelity directly.
# For numerical stability, invert the matrix with larger determinant.
if np.abs(scipy.linalg.det(sigma)) > np.abs(scipy.linalg.det(rho)):
rho, sigma = sigma, rho
eigvals, eigvecs = np.linalg.eigh(rho)
eigvals = np.maximum(eigvals, 0)
tol = np.finfo(float).eps * rho.shape[0] * max(1.0, np.max(eigvals))
# The closed-form geometric mean rho^{1/2} sqrt(rho^{-1/2} sigma rho^{-1/2}) rho^{1/2} is exact
# whenever the inverted state `rho` (the larger-determinant one after the swap above) is
# invertible, even if `sigma` is singular. When `rho` is singular as well -- i.e. both states are
# rank deficient -- the pseudoinverse discards the part of `sigma` outside the support of `rho`
# and no longer reproduces the limit definition (for example it returns 1/sqrt(2) instead of 0
# for |0><0| and |+><+|, whose supports intersect trivially). Fall back to the exact geometric-mean
# characterization tr(rho # sigma) = max{ tr(W) : [[rho, W], [W, sigma]] >> 0 } in that case, which
# matches the semidefinite-programming branch above and handles non-commuting supports correctly.
if np.min(eigvals) <= tol:
w_var = cvxpy.Variable(rho.shape, hermitian=True)
objective = cvxpy.Maximize(cvxpy.real(cvxpy.trace(w_var)))
constraints = [cvxpy.bmat([[rho, w_var], [w_var, sigma]]) >> 0]
return cvxpy.Problem(objective, constraints).solve(solver=cvxpy.CLARABEL)
sqrt_eigvals = np.sqrt(eigvals)
sq_rho = (eigvecs * sqrt_eigvals) @ eigvecs.conj().T
sqinv_rho = (eigvecs / sqrt_eigvals) @ eigvecs.conj().T
# The inner matrix `M = rho^{-1/2} sigma rho^{-1/2}` is Hermitian positive semidefinite, so its
# principal square root follows from an eigendecomposition. This is cheaper and more stable than
# a general `scipy.linalg.sqrtm`, and avoids the LinAlgWarning it raises on rank-deficient `sigma`.
m_mat = sqinv_rho @ sigma @ sqinv_rho
m_mat = (m_mat + m_mat.conj().T) / 2
m_eigvals, m_eigvecs = np.linalg.eigh(m_mat)
sqrt_m = (m_eigvecs * np.sqrt(np.maximum(m_eigvals, 0))) @ m_eigvecs.conj().T
sq_mfid = sq_rho @ sqrt_m @ sq_rho
return np.real(np.trace(sq_mfid))
|