Random¶

The submodule numqi.random is designed to generate various random distribution, like random Hermitian matrix, random state, random unitary matrix etc.

In [1]:

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

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

Quantum state¶

Quantum states are represented by 1-dimensional complex vectors with unit norm.

TODO: explain Haar measure

In [2]:

psi = numqi.random.rand_haar_state(3)
print('psi: ', np.round(psi, 3))
print('norm(psi):', np.linalg.norm(psi))

psi = numqi.random.rand_haar_state(3) print('psi: ', np.round(psi, 3)) print('norm(psi):', np.linalg.norm(psi))

psi:  [0.151-0.858j 0.157+0.227j 0.349-0.205j]
norm(psi): 1.0

Sometimes, we need to generate bi-partite quantum states in the Hilbert space $\mathcal{H}_A \otimes \mathcal{H}_B$. numqi still represents these state as 1-dimensional complex vectors, but the dimension of the vector is the product of the dimensions of $\mathcal{H}_A$ and $\mathcal{H}_B$. The basis is ordered as

$$ \{|a_1\rangle \otimes |b_1\rangle, |a_1\rangle \otimes |b_2\rangle, \cdots, |a_2\rangle \otimes |b_1\rangle, |a_2\rangle \otimes |b_2\rangle, \cdots\}. $$

For such a ordering, one can call reshape() function to transform the state into a 2-dimensional array with the first dimension for $A$ and the second dimension for $B$.

In [3]:

dimA = 2
dimB = 3
psiAB = numqi.random.rand_bipartite_state(dimA, dimB)
print('psiAB:', np.round(psiAB, 3), sep='\n')
print('norm(psiAB):', np.linalg.norm(psiAB))

psi_A_B_2d = psiAB.reshape(dimA, dimB)
print('psi_A_B_2d:', np.round(psi_A_B_2d, 3), sep='\n')

dimA = 2 dimB = 3 psiAB = numqi.random.rand_bipartite_state(dimA, dimB) print('psiAB:', np.round(psiAB, 3), sep='\n') print('norm(psiAB):', np.linalg.norm(psiAB)) psi_A_B_2d = psiAB.reshape(dimA, dimB) print('psi_A_B_2d:', np.round(psi_A_B_2d, 3), sep='\n')

psiAB:
[ 0.151+0.038j  0.193-0.073j -0.78 +0.223j -0.299+0.159j -0.214+0.158j
  0.188-0.233j]
norm(psiAB): 0.9999999999999999
psi_A_B_2d:
[[ 0.151+0.038j  0.193-0.073j -0.78 +0.223j]
 [-0.299+0.159j -0.214+0.158j  0.188-0.233j]]

Density matrix¶

Density matrix is a positive semi-defintie matrix with unit trace. numqi represents density matrix as a 2-dimensional complex array. The diagonal elements are real and non-negative, and the off-diagonal elements are complex conjugate to each other. The trace of the density matrix is 1.

In [4]:

rho = numqi.random.rand_density_matrix(3)

print('rho:', np.round(rho, 3), sep='\n')
print('trace(rho):', np.trace(rho).real)
# positive semi-definite matrix has non-negative eigenvalues
print('np.linalg.eigvalsh(rho):', np.linalg.eigvalsh(rho))

rho = numqi.random.rand_density_matrix(3) print('rho:', np.round(rho, 3), sep='\n') print('trace(rho):', np.trace(rho).real) # positive semi-definite matrix has non-negative eigenvalues print('np.linalg.eigvalsh(rho):', np.linalg.eigvalsh(rho))

rho:
[[0.261+0.j    0.129+0.251j 0.053+0.108j]
 [0.129-0.251j 0.318+0.j    0.177+0.034j]
 [0.053-0.108j 0.177-0.034j 0.421+0.j   ]]
trace(rho): 1.0
np.linalg.eigvalsh(rho): [0.00240135 0.27339592 0.72420273]

Separable density matrix can also be generated by numqi.random module.

In [5]:

dimA = 3
dimB = 3
rho = numqi.random.rand_separable_dm(dimA, dimB)
print('shape: ', rho.shape)

dimA = 3 dimB = 3 rho = numqi.random.rand_separable_dm(dimA, dimB) print('shape: ', rho.shape)

shape:  (9, 9)

Symmetric and Hermitian matrix¶

Hermitian matrix is a complex matrix that is equal to its conjugate transpose. Symmetric matrix is a special case of Hermitian matrix, where the matrix is real.

In [6]:

mat_hermitian = numqi.random.rand_hermitian_matrix(3)
print('mat_hermitian:', np.round(mat_hermitian, 3), sep='\n')

mat_symmetric = numqi.random.rand_hermitian_matrix(3, tag_complex=False)
print('mat_symmetric:', np.round(mat_symmetric, 3), sep='\n')

mat_hermitian = numqi.random.rand_hermitian_matrix(3) print('mat_hermitian:', np.round(mat_hermitian, 3), sep='\n') mat_symmetric = numqi.random.rand_hermitian_matrix(3, tag_complex=False) print('mat_symmetric:', np.round(mat_symmetric, 3), sep='\n')

mat_hermitian:
[[-1.891+0.j    -0.316-0.441j  0.296+1.905j]
 [-0.316+0.441j -2.222+0.j    -1.195+0.389j]
 [ 0.296-1.905j -1.195-0.389j -1.542+0.j   ]]
mat_symmetric:
[[ 2.324 -0.443 -1.211]
 [-0.443 -1.049 -0.145]
 [-1.211 -0.145 -0.725]]

The parameter eig=(min_eig,max_eig) is to control the range of eigenvalues of the generated Hermitian matrix. The default value is eig=None which means the eigenvalues are not constrained. With this way, one can also generate positive semi-definite matrix by setting eig=(0,1).

In [7]:

mat_positive = numqi.random.rand_hermitian_matrix(3, eig=(0,1))
print('mat_positive:', np.round(mat_positive, 3), sep='\n')
print('np.linalg.eigvalsh(mat_positive):', np.linalg.eigvalsh(mat_positive))

mat_positive = numqi.random.rand_hermitian_matrix(3, eig=(0,1)) print('mat_positive:', np.round(mat_positive, 3), sep='\n') print('np.linalg.eigvalsh(mat_positive):', np.linalg.eigvalsh(mat_positive))

mat_positive:
[[ 0.535+0.j     0.015-0.02j  -0.211+0.252j]
 [ 0.015+0.02j   0.881+0.j     0.026+0.155j]
 [-0.211-0.252j  0.026-0.155j  0.486+0.j   ]]
np.linalg.eigvalsh(mat_positive): [0.16080777 0.76312641 0.97788752]

Unitary matrix¶

Unitary matrix is a complex matrix whose conjugate transpose is its inverse. numqi.random.rand_haar_unitary generates unitary matrix with Haar measure.

In [8]:

mat_unitary = numqi.random.rand_haar_unitary(3)
print('mat_unitary:', np.round(mat_unitary, 3), sep='\n')
print('mat_unitary @ mat_unitary.T.conj():', np.round(mat_unitary @ mat_unitary.T.conj(), 3), sep='\n')

mat_unitary = numqi.random.rand_haar_unitary(3) print('mat_unitary:', np.round(mat_unitary, 3), sep='\n') print('mat_unitary @ mat_unitary.T.conj():', np.round(mat_unitary @ mat_unitary.T.conj(), 3), sep='\n')

mat_unitary:
[[ 0.162-0.375j -0.255-0.745j  0.325-0.328j]
 [-0.824-0.287j  0.444-0.195j -0.048-0.028j]
 [-0.154+0.218j  0.065+0.375j  0.722-0.513j]]
mat_unitary @ mat_unitary.T.conj():
[[1.+0.j 0.-0.j 0.+0.j]
 [0.+0.j 1.+0.j 0.-0.j]
 [0.-0.j 0.+0.j 1.+0.j]]

numqi.random.rand_special_orthogonal_matrix provides more options over numqi.random.rand_haar_unitary but it's may not from Haar measure. Orthogonal matrix is a special case of unitary matrix, where the matrix is real. One can generate orthogonal matrix by setting tag_complex=False (default to False) in numqi.random.rand_special_orthogonal_matrix.

In [9]:

mat_orthonormal = numqi.random.rand_special_orthogonal_matrix(3, tag_complex=False)
print('mat_orthonormal:', np.round(mat_orthonormal, 3), sep='\n')
print('mat_orthonormal @ mat_orthonormal.T:', np.round(mat_orthonormal @ mat_orthonormal.T, 3), sep='\n')

mat_orthonormal = numqi.random.rand_special_orthogonal_matrix(3, tag_complex=False) print('mat_orthonormal:', np.round(mat_orthonormal, 3), sep='\n') print('mat_orthonormal @ mat_orthonormal.T:', np.round(mat_orthonormal @ mat_orthonormal.T, 3), sep='\n')

mat_orthonormal:
[[ 0.969 -0.221 -0.109]
 [ 0.151  0.881 -0.448]
 [ 0.195  0.418  0.887]]
mat_orthonormal @ mat_orthonormal.T:
[[ 1. -0. -0.]
 [-0.  1. -0.]
 [-0. -0.  1.]]

Quantum channel¶

Kraus operator is a set of matrices that satisfies the completeness relation $\sum_i K_i^\dagger K_i = I$.

In [10]:

kop = numqi.random.rand_kraus_op(num_term=4, dim_in=3, dim_out=5)
print('kop.shape:', kop.shape)
tmp0 = np.einsum(kop, [0,1,2], kop.conj(), [0,1,3], [2,3], optimize=True)
print('sum_i K_i K_i^dagger:', np.round(tmp0, 3), sep='\n')

kop = numqi.random.rand_kraus_op(num_term=4, dim_in=3, dim_out=5) print('kop.shape:', kop.shape) tmp0 = np.einsum(kop, [0,1,2], kop.conj(), [0,1,3], [2,3], optimize=True) print('sum_i K_i K_i^dagger:', np.round(tmp0, 3), sep='\n')

kop.shape: (4, 5, 3)
sum_i K_i K_i^dagger:
[[1.+0.j 0.+0.j 0.-0.j]
 [0.+0.j 1.+0.j 0.+0.j]
 [0.+0.j 0.-0.j 1.+0.j]]

Choi operator is a matrix that is related to the Kraus operator. For the channel mapping $\mathcal{H}_{d_i}\to\mathcal{H}_{d_o}$, the shape of the Choi operator is $(d_i,d_o,d_i,d_o)$.

In [11]:

dim_in = 3
dim_out = 5
choi_op = numqi.random.rand_choi_op(dim_in, dim_out)
print('choi_op.shape:', choi_op.shape)
tmp0 = np.trace(choi_op.reshape(dim_in, dim_out, dim_in, dim_out), axis1=1, axis2=3)
print('partial-trace(choi_op):', np.round(tmp0, 3), sep='\n')

dim_in = 3 dim_out = 5 choi_op = numqi.random.rand_choi_op(dim_in, dim_out) print('choi_op.shape:', choi_op.shape) tmp0 = np.trace(choi_op.reshape(dim_in, dim_out, dim_in, dim_out), axis1=1, axis2=3) print('partial-trace(choi_op):', np.round(tmp0, 3), sep='\n')

choi_op.shape: (15, 15)
partial-trace(choi_op):
[[1.-0.j 0.+0.j 0.+0.j]
 [0.-0.j 1.+0.j 0.-0.j]
 [0.-0.j 0.+0.j 1.+0.j]]

Finite Field 2¶

Pauli Group and Clifford group are closely related to algebra over finite field 2. See tutorial/Clifford-group for more details. numqi.random provides functions to generate random elements in these groups. Please also check paper "How to efficiently select an arbitrary clifford group element" doi-link for the algorithm details.

In [12]:

size = 4
bitstr = numqi.random.rand_F2(size)
print('bitstr:', bitstr)

# set not_zero=True to avoid generate zero vector
for _ in range(5):
    print('non_zero=True:', numqi.random.rand_F2(2, not_zero=True))

size = 4 bitstr = numqi.random.rand_F2(size) print('bitstr:', bitstr) # set not_zero=True to avoid generate zero vector for _ in range(5): print('non_zero=True:', numqi.random.rand_F2(2, not_zero=True))

bitstr: [0 0 1 1]
non_zero=True: [1 0]
non_zero=True: [0 1]
non_zero=True: [1 0]
non_zero=True: [0 1]
non_zero=True: [0 1]

Symplectic matrix is a matrix that preserves the symplectic form. The matrix $A$ is called symplectic if it satisfies

$$ A \Lambda A^T= \Lambda \pmod{2}$$

where $\Lambda=\sigma_x\otimes I_n$ and every matrix element of $A$ is either 0 or 1. All symplectic matrices form a group called symplectic group $Sp(n, F_2)$ which is shorted as spf2 in numqi, e.g. numqi.group.spf2 module.

In [13]:

n = 2
sigma_x = np.array([[0, 1], [1, 0]], dtype=np.uint8)
mat_lambda = np.kron(sigma_x, np.eye(n, dtype=np.uint8))
matA = numqi.random.rand_SpF2(n)

print('matA:', matA, sep='\n')
print('Lambda:', mat_lambda, sep='\n')
tmp0 = (matA @ mat_lambda @ matA.T) % 2
print('matA @ Lambda @ matA.T:', tmp0, sep='\n')

n = 2 sigma_x = np.array([[0, 1], [1, 0]], dtype=np.uint8) mat_lambda = np.kron(sigma_x, np.eye(n, dtype=np.uint8)) matA = numqi.random.rand_SpF2(n) print('matA:', matA, sep='\n') print('Lambda:', mat_lambda, sep='\n') tmp0 = (matA @ mat_lambda @ matA.T) % 2 print('matA @ Lambda @ matA.T:', tmp0, sep='\n')

matA:
[[0 0 1 0]
 [0 0 1 1]
 [1 1 1 0]
 [0 1 1 1]]
Lambda:
[[0 0 1 0]
 [0 0 0 1]
 [1 0 0 0]
 [0 1 0 0]]
matA @ Lambda @ matA.T:
[[0 0 1 0]
 [0 0 0 1]
 [1 0 0 0]
 [0 1 0 0]]

Dirichlet distribution¶

wiki/Dirichlet-distribution

related paper "Improving application performance with biased distributions of quantum states" doi-link

In [14]:

alpha = np.array([0.3, 0.3, 0.3])
N0 = 5000
np0 = np_rng.dirichlet(alpha, size=N0)
assert np.abs(np0.sum(axis=1)-1).max()<1e-7

tmp0 = [0, np.pi*2/3, np.pi*4/3]
tmp1 = np0 @ np.array([[np.cos(x),np.sin(x)] for x in tmp0])
fig,ax = plt.subplots()
ax.scatter(tmp1[:,0], tmp1[:,1], s=3, alpha=0.3)
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_aspect('equal')

alpha = np.array([0.3, 0.3, 0.3]) N0 = 5000 np0 = np_rng.dirichlet(alpha, size=N0) assert np.abs(np0.sum(axis=1)-1).max()<1e-7 tmp0 = [0, np.pi*2/3, np.pi*4/3] tmp1 = np0 @ np.array([[np.cos(x),np.sin(x)] for x in tmp0]) fig,ax = plt.subplots() ax.scatter(tmp1[:,0], tmp1[:,1], s=3, alpha=0.3) ax.set_xlim(-1, 1) ax.set_ylim(-1, 1) ax.set_aspect('equal')