Numerical range¶
Given a complex matrix $A\in\mathbb{C}^{d\times d}$, its numerical range $W(A)\subset \mathbb{C}$ is defined as the set
$$ W(A)=\left\{ x^{\dagger}Ax:x\in\mathbb{C}^{d},\lVert x\rVert_{2}=1\right\} $$
and Hausdorff–Toeplitz theorem proved that the numerical range is convex and compact.
Some properties of the numerical range are listed below
- $\forall \alpha,\beta\in\mathbb{C}, W(\alpha I + \beta F)=\alpha+\beta W(F)$
- $\forall UU^\dagger=U^\dagger U=I, W(F)=W(U^\dagger F U)$
- $W(F^*)=\{\lambda^*:\lambda \in W(F)\}$
Check "Numerical Ranges and Applications in Quantum Information" link for a detailed introduction.
import numpy as np
import matplotlib.pyplot as plt
try:
import numqi
except ImportError:
%pip install numqi
import numqi
np_rng = np.random.default_rng()
hf_randc = lambda *x: np_rng.normal(size=x) + 1j*np_rng.normal(size=x)
Numerical range of normal matrix¶
We said a matrix $A$ is normal iff $AA^{\dagger}=A^{\dagger}A$. And normal matrix is diagonalizable by unitary matrix
$$ A=U\Sigma U^\dagger $$
where $U$ is some unitary matrix and $\Sigma$ is a diagonal matrix. Two typical kinds of normal matrices are Hermitian matrix and unitary matrix. For the Hermitian matrix, its eigenvalues are real and for the unitary matrix, its eigenvalues are complex with unit modulus.
A good property for normal matrix is that its numerical range is the convex hull of its eigenvalues. Below we show this numerically for a random normal matrix.
N0 = 4
matU = numqi.random.rand_special_orthogonal_matrix(N0, tag_complex=True)
EVL = hf_randc(N0)
matA = (matU * EVL) @ matU.T.conj()
boundary = numqi.matrix_space.get_matrix_numerical_range(matA, num_point=200)
fig,ax = plt.subplots()
ax.plot(boundary.real, boundary.imag)
ind0 = slice(None,None,5)
ax.plot(boundary.real[ind0], boundary.imag[ind0], '.', color='r', markersize=10)
ax.plot(EVL.real, EVL.imag, '+')
for ind0,x in enumerate(EVL):
ax.text(x.real, x.imag, str(ind0))
The numerical range denotes the area enclosed by the blue curve and the eigenvalues are the red dots labeled by $0,1,2,3$. One can see that boundary of the numerical range are segments of lines connecting the eigenvalues.
Task: think about how should the numerical range of a Hermitian (unitary) matrix look like.
Numerical range of abnormal matrix¶
For an abnormal matrix ($AA^\dagger\ne A^\dagger A$), its numerical range can be quite complicated, but still we can compute it numerically.
algorithms: for $\theta\in[0,2\pi]$, compute the minimum (maximum) eigenvalue of
$$ e^{i\theta}A+e^{-i\theta}A^\dagger $$
and the associated eigenvector $v$. Then $v^\dagger Av$ gives one point on the boundary and the numerical range is the convex hull of these points.
Below, we show the numerical range of a random matrix $A$ which is the direct sum of a normal part and abnormal part.
seed = np_rng.integers(0, 2**32-1)
np_rng = np.random.default_rng(seed=seed) #1457230309
hf_randc = lambda *x: np_rng.normal(size=x) + 1j*np_rng.normal(size=x)
N0 = 5
N1 = 2
matU = numqi.random.rand_special_orthogonal_matrix(N0, tag_complex=True, seed=np_rng)
EVL = hf_randc(N1)
matIN = hf_randc(N0-N1, N0-N1) #abnormal matrix
tmp0 = (matU[:,:N1]*EVL) @ matU[:,:N1].T.conj()
tmp1 = matU[:,N1:] @ matIN @ matU[:,N1:].T.conj()
matA0 = tmp0*0 + tmp1
matA1 = tmp0 + tmp1
boundary0 = numqi.matrix_space.get_matrix_numerical_range(matA0, num_point=200)
boundary1 = numqi.matrix_space.get_matrix_numerical_range(matA1, num_point=200)
fig,ax = plt.subplots()
ax.plot(boundary0.real, boundary0.imag, label='NR(abnormal-part of $A$)')
ax.plot(boundary1.real, boundary1.imag, label='NR($A$)')
ind0 = slice(None,None)
ax.plot(boundary1.real[ind0], boundary1.imag[ind0], '.', markersize=10)
ax.plot(EVL.real, EVL.imag, '+', label='eigenvalue($A$)')
for ind0,x in enumerate(EVL):
ax.text(x.real, x.imag, str(ind0))
ax.legend()
<matplotlib.legend.Legend at 0x7fcfb589f7a0>
One could still see some pattern of the numerical range of the normal part.
Task: run the code above again (every time you run it, the matrix $A$ is different) and observe various numerical range.
A typical example is the following $2\times 2$ abnormal matrix from thesis (page-10) "Numerical Ranges and Applications in Quantum Information" link
$$ A=\begin{bmatrix} i&1\\0&-i \end{bmatrix} $$
Task: can you find out which unit vector (see definition of numerical range) gives the maximum imaginary value of $v^\dagger Av$?
matA = np.array([[1j,1], [0,-1j]])
# matA = np.array([[0,2], [0,0]])
boundary = numqi.matrix_space.get_matrix_numerical_range(matA, num_point=50)
fig,ax = plt.subplots()
ax.plot(boundary.real, boundary.imag)
ind0 = slice(None,None)
ax.plot(boundary.real[ind0], boundary.imag[ind0], '.', color='r', markersize=10)
ax.set_aspect('equal')
Joint algebraic numerical range¶
For a series of matrices $A_1,A_2,\cdots,A_m\in\mathbb{C}^{d\times d}$, their joint numerical range $W(A_1,\cdots,A_m)$ is defined as
$$ W(A_{1},A_{2},\cdots,A_{r})=\left\{ a\in\mathbb{C}^{m}:x\in\mathbb{C}^{d},\lVert x\rVert_{2}=1,a_{i}=x^{\dagger}A_{i}x\right\} $$
and their joint algebraic numerical range (JANR) $L(A_1,\cdots,A_m)$ is defined as
$$ L(A_{1},A_{2},\cdots,A_{r})=\left\{ a\in\mathbb{C}^{r}:\rho\in\mathbb{C}^{d\times d},\rho\succeq0,\mathrm{Tr}[\rho]=1,a_{i}=\mathrm{Tr}[A_{i}\rho]\right\} $$
Joint numerical range is usually non-convex and is star-shaped when $d$ is large enough. While joint algebraic numerical range is convex and compact.
One may notice that $x$ in above equations is simply pure state and $\rho$ is a density matrix. When the operator $A_i$ all are hermitian, one can interpret $W(A_1,\cdots,A_m)$ as the set of all possible expectation values for pure state and $L(A_1,\cdots,A_m)$ as the set of all possible expectation values for mixed state (density matrix).
One can solve the boundary of the joint algebraic numerical range by semidefinite programming
$$ \max\;\beta $$
$$ s.t.\;\begin{cases} \rho\succeq0\\ \mathrm{Tr}[\rho]=1\\ \mathrm{Tr}[\rho A_{i}]=\beta\hat{n}_{i} & i=1,\cdots,m \end{cases} $$
where $\beta$ is the distance to the boundary from the origin.
A special case is that for two Hermitian operators $m=2$, the joint numerical range coincide with the joint algebraic numerical range, and also with the numerical range of the operator $A_1+iA_2$.
Below we are going to reproduce results of the paper "Discontinuity of maximum entropy inference and quantum phase transitions" doi-link
theta_list = np.linspace(0, 2*np.pi, 200)
direction = np.stack([np.cos(theta_list), np.sin(theta_list)], axis=1)
op0 = np.diag([1,1,-1])
op1_list = [
np.array([[1,0,1], [0,1,1], [1,1,-1]]),
np.array([[1,0,1], [0,0,1], [1,1,-1]]),
np.diag([1,0,-1]),
]
fig,ax_list = plt.subplots(1, 3, figsize=(9,4))
for ax,op1 in zip(ax_list,op1_list):
beta_list,op_nr_list,dual_value_list = numqi.matrix_space.get_joint_algebraic_numerical_range([op0,op1], direction, return_info=True)
tmp0 = slice(None, None, 3)
tmp1 = np.angle(dual_value_list[:,0] + 1j*dual_value_list[:,1])
numqi.matrix_space.draw_line_list(ax, op_nr_list[tmp0], tmp1[tmp0], kind='tangent', color='#CCCCCC', radius=1)
ax.plot(op_nr_list[:,0], op_nr_list[:,1])
ax.set_aspect('equal')
We also plot the supporting plane which is obtained from the dual variable in the semidefinite programming.
fig,ax_list = plt.subplots(1, 3, figsize=(9,4))
for ax,op1 in zip(ax_list,op1_list):
boundary0 = numqi.matrix_space.get_matrix_numerical_range(op0+1j*op1, num_point=200)
ax.plot(boundary0.real, boundary0.imag)
ax.set_aspect('equal')
Obviously, we can get the same boundary from the numerical range of $A_1+iA_2$.
Furthermore, we will use this technique to do the entanglement detection.