1. Basics¶

1.1. Basic Representation¶

States and operators in quimb are simply dense numpy arrays or sparse scipy matrices. All functions should directly work with these but the class qarray is also provided as a very thin subclass of numpy.ndarray with a few helpful methods and attributes. The quimbify() function (aliased to qu()) can convert between the various representations.

from quimb import *
data = [1, 2j, -3]

Kets are column vectors, i.e. with shape (d, 1):

qu(data, qtype='ket')

[[ 1.+0.j]
 [ 0.+2.j]
 [-3.+0.j]]

The normalized=True option can be used to ensure a normalized output.

Bras are row vectors, i.e. with shape (1, d):

qu(data, qtype='bra')  # also conjugates the data

[[ 1.-0.j  0.-2.j -3.-0.j]]

And operators are square matrices, i.e. have shape (d, d):

qu(data, qtype='dop')

[[ 1.+0.j  0.-2.j -3.+0.j]
 [ 0.+2.j  4.+0.j  0.-6.j]
 [-3.+0.j  0.+6.j  9.+0.j]]

Which can also be sparse:

qu(data, qtype='dop', sparse=True)

<3x3 sparse matrix of type '<class 'numpy.complex128'>'
	with 9 stored elements in Compressed Sparse Row format>

The sparse format can be specified with the stype keyword. The partial function versions of each of the above are also available:

Note

If a simple 1d-list is supplied and no qtype is given, 'ket' is assumed.

1.2. Basic Operations¶

The ‘dagger’, or hermitian conjugate, operation is performed with the .H attribute:

psi = 1.0j * bell_state('psi-')
psi

[[ 0.+0.j      ]
 [ 0.+0.707107j]
 [-0.-0.707107j]
 [ 0.+0.j      ]]

psi.H

[[ 0.-0.j        0.-0.707107j -0.+0.707107j  0.-0.j      ]]

This is just the combination of .conj() and .T, but only available for scipy.sparse matrices and qarray s (not numpy.ndarray s).

The product of two quantum objects is the dot or matrix product, which, since python 3.5, has been overloaded with the @ symbol. Using it is recommended:

psi = up()
psi

[[1.+0.j]
 [0.+0.j]]

psi.H @ psi  # inner product

[[1.+0.j]]

X = pauli('X')
X @ psi  # act as gate

[[0.+0.j]
 [1.+0.j]]

(psi.H @ X @ psi)  # operator expectation

[[0.+0.j]]

Scalar expectation values might best be computed using the expectation() function (aliased to expec()) which dispatches to accelerated methods:

expec(psi, psi)

1.0

expec(psi, X)

0j

Here’s an example for a much larger (20 qubit), sparse operator expecation, which will be automatically parallelized:

psi = rand_ket(2**20)
A = rand_herm(2**20, sparse=True) + speye(2**20)
A

<1048576x1048576 sparse matrix of type '<class 'numpy.complex128'>'
	with 11534268 stored elements in Compressed Sparse Row format>

expec(A, psi)  # should be ~ 1

0.9998770131258168

%%timeit
expec(A, psi)

59 ms ± 3.85 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

1.3. Combining Objects - Tensor Products¶

There are a number of ways to combine states and operators, i.e. tensoring them together.

Functional form using kron():

>>> kron(psi1, psi2, psi3, ...)
...

This can also be done using the & overload on qarray and scipy matrices:

>>> psi1 & psi2 & psi3
...

Warning

When quimb is imported, it monkey patches the otherwise unused method of &/__and__ of scipy sparse matrices to kron().

Often one wants to sandwich an operator with many identities, ikron() can be used for this:

dims = [2] * 10  # overall space of 10 qubits
X = pauli('X')
IIIXXIIIII = ikron(X, dims, inds=[3, 4])  # act on 4th and 5th spin only
IIIXXIIIII.shape

(1024, 1024)

For more advanced tensor constructions, such as reversing and interleaving identities within operators pkron() can be used:

dims = [2] * 3
XZ = pauli('X') & pauli('Z')
ZIX = pkron(XZ, dims, inds=[2, 0])
ZIX.real.astype(int)

[[ 0  1  0  0  0  0  0  0]
 [ 1  0  0  0  0  0  0  0]
 [ 0  0  0  1  0  0  0  0]
 [ 0  0  1  0  0  0  0  0]
 [ 0  0  0  0  0 -1  0  0]
 [ 0  0  0  0 -1  0  0  0]
 [ 0  0  0  0  0  0  0 -1]
 [ 0  0  0  0  0  0 -1  0]]

ZIX would then act with Z on first spin, and X on 3rd.

1.4. Removing Objects - Partial Trace¶

To remove, or ignore, certain parts of a quantum state the partial trace function partial_trace() (aliased to ptr()) is used. Here, the internal dimensions of a state must be supplied as well as the indicies of which of these subsystems to keep.

For example, if we have a random system of 10 qubits (hilbert space of dimension 2**10), and we want just the reduced density matrix describing the first and last spins:

dims = [2] * 10
D = prod(dims)
psi = rand_ket(D)
rho_ab = ptr(psi, dims, [0, 9])
rho_ab.round(3)  # probably pretty close to identity

[[ 0.256+0.j    -0.013-0.007j  0.02 +0.008j -0.014+0.007j]
 [-0.013+0.007j  0.263+0.j    -0.011+0.008j  0.016+0.004j]
 [ 0.02 -0.008j -0.011-0.008j  0.231+0.j     0.01 +0.008j]
 [-0.014-0.007j  0.016-0.004j  0.01 -0.008j  0.25 +0.j   ]]

partial_trace() accepts dense or sparse, operators or vectors.