Generating Random Quantum States & Operators¶

QuTiP includes a collection of random state, unitary and channel generators for simulations, Monte Carlo evaluation, theorem evaluation, and code testing. Each of these objects can be sampled from one of several different distributions including the default distributions used by QuTiP versions prior to 3.2.0.

For example, a random Hermitian operator can be sampled by calling rand_herm function:

In [1]: rand_herm(5)
Out[1]: 
Quantum object: dims = [[5], [5]], shape = (5, 5), type = oper, isherm = True
Qobj data =
[[ 0.21008894+0.j         -0.46056584-0.47165765j  0.        +0.j
  -0.15248215-0.48974871j  0.        +0.j        ]
 [-0.46056584+0.47165765j  0.18241967+0.j          0.        +0.j
  -0.20212601+0.15804465j  0.01329707-0.42871682j]
 [ 0.        +0.j          0.        +0.j          0.        +0.j
   0.        +0.j         -0.0608447 +0.38781155j]
 [-0.15248215+0.48974871j -0.20212601-0.15804465j  0.        +0.j
   0.        +0.j          0.02299724+0.33505856j]
 [ 0.        +0.j          0.01329707+0.42871682j -0.0608447 -0.38781155j
   0.02299724-0.33505856j  0.        +0.j        ]]

Random Variable Type	Sampling Functions	Dimensions
State vector (`ket`)	rand_ket, rand_ket_haar	\(N \times 1\)
Hermitian operator (`oper`)	rand_herm	\(N \times 1\)
Density operator (`oper`)	rand_dm, rand_dm_hs, rand_dm_ginibre	\(N \times N\)
Unitary operator (`oper`)	rand_unitary, rand_unitary_haar	\(N \times N\)
CPTP channel (`super`)	rand_super, rand_super_bcsz	\((N \times N) \times (N \times N)\)

In all cases, these functions can be called with a single parameter \(N\) that dimension of the relevant Hilbert space. The optional dims keyword argument allows for the dimensions of a random state, unitary or channel to be broken down into subsystems.

In [2]: print(rand_super_bcsz(7).dims)
[[[7], [7]], [[7], [7]]]

In [3]: print(rand_super_bcsz(6, dims=[[[2, 3], [2, 3]], [[2, 3], [2, 3]]]).dims)
[[[2, 3], [2, 3]], [[2, 3], [2, 3]]]

Several of the distributions supported by QuTiP support additional parameters as well, namely density and rank. In particular, the rand_herm and rand_dm functions return quantum objects such that a fraction of the elements are identically equal to zero. The ratio of nonzero elements is passed as the density keyword argument. By contrast, the rand_dm_ginibre and rand_super_bcsz take as an argument the rank of the generated object, such that passing rank=1 returns a random pure state or unitary channel, respectively. Passing rank=None specifies that the generated object should be full-rank for the given dimension.

For example,

In [4]: rand_dm(5, density=0.5)
Out[4]: 
Quantum object: dims = [[5], [5]], shape = (5, 5), type = oper, isherm = True
Qobj data =
[[ 0.48467318+0.j          0.02458384-0.09377379j  0.13215612-0.01977424j
  -0.07476071-0.07190307j  0.0647332 +0.08671008j]
 [ 0.02458384+0.09377379j  0.09549324+0.j          0.05185442+0.12098514j
   0.04983752-0.08919698j  0.        +0.j        ]
 [ 0.13215612+0.01977424j  0.05185442-0.12098514j  0.24597381+0.j
  -0.0859455 -0.11157708j  0.        +0.j        ]
 [-0.07476071+0.07190307j  0.04983752+0.08919698j -0.0859455 +0.11157708j
   0.10932586+0.j          0.        +0.j        ]
 [ 0.0647332 -0.08671008j  0.        +0.j          0.        +0.j
   0.        +0.j          0.06453391+0.j        ]]

In [5]: rand_dm_ginibre(5, rank=2)
Out[5]: 
Quantum object: dims = [[5], [5]], shape = (5, 5), type = oper, isherm = True
Qobj data =
[[ 0.39908994+0.j          0.11093243+0.0421428j   0.07558034-0.11219794j
   0.13747317-0.08786001j -0.02270823-0.01477814j]
 [ 0.11093243-0.0421428j   0.07564242+0.j         -0.03007169-0.15174557j
   0.04905063-0.02944462j  0.00753004-0.04216709j]
 [ 0.07558034+0.11219794j -0.03007169+0.15174557j  0.39803457+0.j
   0.0046961 +0.06889389j  0.0977376 +0.07311324j]
 [ 0.13747317+0.08786001j  0.04905063+0.02944462j  0.0046961 -0.06889389j
   0.07895757+0.j         -0.006409  -0.0338791j ]
 [-0.02270823+0.01477814j  0.00753004+0.04216709j  0.0977376 -0.07311324j
  -0.006409  +0.0338791j   0.04827549+0.j        ]]

See the API documentation: Random Operators and States for details.

Warning

When using the density keyword argument, setting the density too low may result in not enough diagonal elements to satisfy trace constraints.

Random objects with a given eigen spectrum¶

Note

New in QuTiP 3.2

It is also possible to generate random Hamiltonian (rand_herm) and densitiy matrices (rand_dm) with a given eigen spectrum. This is done by passing an array of eigenvalues as the first argument to either function. For example,

In [6]: eigs = np.arange(5)

In [7]: H = rand_herm(eigs, density=0.5)

In [8]: H
Out[8]: 
Quantum object: dims = [[5], [5]], shape = (5, 5), type = oper, isherm = True
Qobj data =
[[ 1.5       -1.11022302e-16j -0.93095066-5.08262600e-01j
   0.        +0.00000000e+00j -0.93095066-5.08262600e-01j
   0.        +0.00000000e+00j]
 [-0.93095066+5.08262600e-01j  1.25      -1.11022302e-16j
   0.        +0.00000000e+00j  0.25      -1.11022302e-16j
   0.        +0.00000000e+00j]
 [ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
   3.        +0.00000000e+00j  0.        +0.00000000e+00j
  -1.        +0.00000000e+00j]
 [-0.93095066+5.08262600e-01j  0.25      -1.11022302e-16j
   0.        +0.00000000e+00j  1.25      -1.11022302e-16j
   0.        +0.00000000e+00j]
 [ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
  -1.        +0.00000000e+00j  0.        +0.00000000e+00j
   3.        +0.00000000e+00j]]

In [9]: H.eigenenergies()
Out[9]: 
array([-2.22044605e-16,  1.00000000e+00,  2.00000000e+00,  3.00000000e+00,
        4.00000000e+00])

In order to generate a random object with a given spectrum QuTiP applies a series of random complex Jacobi rotations. This technique requires many steps to build the desired quantum object, and is thus suitable only for objects with Hilbert dimensionality \(\lesssim 1000\).

Composite random objects¶

In many cases, one is interested in generating random quantum objects that correspond to composite systems generated using the qutip.tensor.tensor function. Specifying the tensor structure of a quantum object is done using the dims keyword argument in the same fashion as one would do for a qutip.Qobj object:

In [10]: rand_dm(4, 0.5, dims=[[2,2], [2,2]])
Out[10]: 
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[ 0.14979702+0.j         -0.09953542-0.01733519j  0.05311938+0.05517837j
   0.12900729-0.11820813j]
 [-0.09953542+0.01733519j  0.28659651+0.j          0.        +0.j
   0.00908027+0.16374315j]
 [ 0.05311938-0.05517837j  0.        +0.j          0.15949409+0.j
   0.1012504 -0.06445248j]
 [ 0.12900729+0.11820813j  0.00908027-0.16374315j  0.1012504 +0.06445248j
   0.40411238+0.j        ]]