Random¶
The submodule numqi.random is designed to generate various random distribution, like random Hermitian matrix, random state, random unitary matrix etc.
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
psi = numqi.random.rand_haar_state(3)
print('psi: ', np.round(psi, 3))
print('norm(psi):', np.linalg.norm(psi))
psi: [-0.591+0.589j -0.112-0.454j -0.104+0.274j] 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$.
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.116+0.431j 0.249+0.384j -0.069+0.117j -0.121+0.042j -0.101-0.453j -0.578-0.083j] norm(psiAB): 0.9999999999999999 psi_A_B_2d: [[-0.116+0.431j 0.249+0.384j -0.069+0.117j] [-0.121+0.042j -0.101-0.453j -0.578-0.083j]]
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.
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.193+0.j -0.119+0.06j -0.052+0.145j] [-0.119-0.06j 0.267+0.j -0.185-0.119j] [-0.052-0.145j -0.185+0.119j 0.54 +0.j ]] trace(rho): 1.0000000000000002 np.linalg.eigvalsh(rho): [0.00170315 0.29627569 0.70202116]
Separable density matrix can also be generated by numqi.random module.
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.
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: [[-2.371+0.j -1.001+0.719j 0.405-0.113j] [-1.001-0.719j -3.335+0.j 0.446-1.058j] [ 0.405+0.113j 0.446+1.058j -0.97 +0.j ]] mat_symmetric: [[-0.758 -2.132 -1.914] [-2.132 2.24 -1.659] [-1.914 -1.659 1.094]]
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).
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.223+0.j -0.182-0.041j -0.024-0.091j] [-0.182+0.041j 0.402+0.j 0.083+0.03j ] [-0.024+0.091j 0.083-0.03j 0.08 -0.j ]] np.linalg.eigvalsh(mat_positive): [0.01800032 0.13780532 0.54985224]
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.
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.262+0.266j 0.874+0.045j 0.292+0.097j] [ 0.596+0.082j 0.315+0.358j -0.364-0.527j] [ 0.461+0.535j -0.062+0.052j -0.128+0.692j]] 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.
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.402 -0.519 0.754] [ 0.86 0.498 -0.115] [-0.315 0.695 0.647]] 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$.
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)$.
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.
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: [1 0 1 1] non_zero=True: [1 1] non_zero=True: [1 1] non_zero=True: [0 1] non_zero=True: [0 1] non_zero=True: [1 0]
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.
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 1 1 1] [0 0 1 0] [0 1 0 0] [1 1 1 0]] 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]]
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')
