Source code for qutip.random_objects

# -*- coding: utf-8 -*-
# The above line is so that UTF-8 comments won't break Py2.

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
"""
This module is a collection of random state and operator generators.
The sparsity of the ouput Qobj's is controlled by varing the
`density` parameter.
"""

__all__ = [
    'rand_herm', 'rand_unitary', 'rand_ket', 'rand_dm',
    'rand_unitary_haar', 'rand_ket_haar', 'rand_dm_ginibre',
    'rand_dm_hs', 'rand_super_bcsz', 'rand_stochastic', 'rand_super'
]

from numpy import arcsin, sqrt, pi
from scipy.linalg import sqrtm
import numpy as np
import numpy.linalg as la
import scipy.sparse as sp
from qutip.qobj import Qobj
from qutip.operators import create, destroy, jmat
from qutip.states import basis
import qutip.superop_reps as sr


UNITS = np.array([1, 1j])


def rand_jacobi_rotation(A, seed=None):
    """Random Jacobi rotation of a sparse matrix.

    Parameters
    ----------
    A : spmatrix
        Input sparse matrix.

    Returns
    -------
    spmatrix
        Rotated sparse matrix.
    """
    if seed is not None:
        np.random.seed(seed=seed)
    if A.shape[0]!=A.shape[1]:
        raise Exception('Input matrix must be square.')
    n = A.shape[0]
    angle = 2*np.random.random()*np.pi
    a = 1.0/np.sqrt(2)*np.exp(-1j*angle)
    b = 1.0/np.sqrt(2)*np.exp(1j*angle)
    i = np.int(np.floor(np.random.random()*n))
    j = i
    while (i == j):
        j = np.int(np.floor(np.random.random()*n))
    data = np.hstack((np.array([a,-b,a,b], dtype=complex),
                      np.ones(n-2, dtype=complex)))
    diag = np.delete(np.arange(n), [i,j])
    rows = np.hstack(([i,i,j,j], diag))
    cols = np.hstack(([i,j,i,j], diag))
    R = sp.coo_matrix((data,(rows,cols)), shape=[n,n], dtype=complex).tocsr()
    A = R*A*R.conj().transpose()
    return A

def randnz(shape, norm=1 / np.sqrt(2), seed=None):
    # This function is intended for internal use.
    """
    Returns an array of standard normal complex random variates.
    The Ginibre ensemble corresponds to setting ``norm = 1`` [Mis12]_.

    Parameters
    ----------
    shape : tuple
        Shape of the returned array of random variates.
    norm : float
        Scale of the returned random variates, or 'ginibre' to draw
        from the Ginibre ensemble.
    """
    if seed is not None:
        np.random.seed(seed=seed)
    if norm == 'ginibre':
        norm = 1
    return np.sum(np.random.randn(*(shape + (2,))) * UNITS, axis=-1) * norm


[docs]def rand_herm(N, density=0.75, dims=None, pos_def=False, seed=None): """Creates a random NxN sparse Hermitian quantum object. If 'N' is an integer, uses :math:`H=0.5*(X+X^{+})` where :math:`X` is a randomly generated quantum operator with a given `density`. Else uses complex Jacobi rotations when 'N' is given by an array. Parameters ---------- N : int, list/ndarray If int, then shape of output operator. If list/ndarray then eigenvalues of generated operator. density : float Density between [0,1] of output Hermitian operator. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. pos_def : bool (default=False) Return a positive semi-definite matrix (by diagonal dominance). seed : int seed for the random number generator Returns ------- oper : qobj NxN Hermitian quantum operator. Note ---- If given a list/ndarray as input 'N', this function returns a random Hermitian object with eigenvalues given in the list/ndarray. This is accomplished via complex Jacobi rotations. While this method is ~50% faster than the corresponding (real only) Matlab code, it should not be repeatedly used for generating matrices larger than ~1000x1000. """ if seed is not None: np.random.seed(seed=seed) if isinstance(N, (np.ndarray,list)): M = sp.diags(N,0, dtype=complex, format='csr') N = len(N) if dims: _check_dims(dims, N, N) nvals = max([N**2 * density, 1]) M = rand_jacobi_rotation(M) while M.nnz < 0.95 * nvals: M = rand_jacobi_rotation(M) M.sort_indices() elif isinstance(N, (int, np.int32, np.int64)): if dims: _check_dims(dims, N, N) if density < 0.5: M = _rand_herm_sparse(N, density, pos_def) else: M = _rand_herm_dense(N, density, pos_def) else: raise TypeError('Input N must be an integer or array_like.') if dims: return Qobj(M, dims=dims) else: return Qobj(M)
def _rand_herm_sparse(N, density, pos_def): target = (1-(1-density)**0.5) num_elems = (N**2 - 0.666 * N) * target + 0.666 * N * density num_elems = max([num_elems, 1]) num_elems = np.int(num_elems) data = (2 * np.random.rand(num_elems) - 1) + \ (2 * np.random.rand(num_elems) - 1) * 1j row_idx, col_idx = zip(*[divmod(index, N) for index in np.random.choice(N*N, num_elems, replace=False)]) M = sp.coo_matrix((data, (row_idx,col_idx)), dtype=complex, shape=(N,N)).tocsr() M = 0.5 * (M + M.conj().transpose()) if pos_def: M.setdiag(np.abs(M.diagonal())+np.sqrt(2)*N) M.sort_indices() return M def _rand_herm_dense(N, density, pos_def): M = (2 * np.random.rand(N,N) - 1) + \ (2 * np.random.rand(N,N) - 1) * 1j M = 0.5 * (M + M.conj().transpose()) target = (1-(density)**0.5) num_remove = N * (N - 0.666) * target + 0.666 * N * (1 - density) num_remove = max([num_remove, 1]) num_remove = np.int(num_remove) for row, col in [divmod(index, N) for index in np.random.choice(N*N, num_remove, replace=False)]: M[col, row] = 0 M[row, col] = 0 if pos_def: as_vec = M.ravel() M[::N+1] = (np.abs(M.diagonal())+np.sqrt(2)*N) return M
[docs]def rand_unitary(N, density=0.75, dims=None, seed=None): """Creates a random NxN sparse unitary quantum object. Uses :math:`\exp(-iH)` where H is a randomly generated Hermitian operator. Parameters ---------- N : int Shape of output quantum operator. density : float Density between [0,1] of output Unitary operator. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- oper : qobj NxN Unitary quantum operator. """ if dims: _check_dims(dims, N, N) U = (-1.0j * rand_herm(N, density, seed=seed)).expm() U.data.sort_indices() if dims: return Qobj(U, dims=dims, shape=[N, N]) else: return Qobj(U)
[docs]def rand_unitary_haar(N=2, dims=None, seed=None): """ Returns a Haar random unitary matrix of dimension ``dim``, using the algorithm of [Mez07]_. Parameters ---------- N : int Dimension of the unitary to be returned. dims : list of lists of int, or None Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- U : Qobj Unitary of dims ``[[dim], [dim]]`` drawn from the Haar measure. """ if dims is not None: _check_dims(dims, N, N) else: dims = [[N], [N]] # Mez01 STEP 1: Generate an N × N matrix Z of complex standard # normal random variates. Z = randnz((N, N), seed=seed) # Mez01 STEP 2: Find a QR decomposition Z = Q · R. Q, R = la.qr(Z) # Mez01 STEP 3: Create a diagonal matrix Lambda by rescaling # the diagonal elements of R. Lambda = np.diag(R).copy() Lambda /= np.abs(Lambda) # Mez01 STEP 4: Note that R' := Λ¯¹ · R has real and # strictly positive elements, such that # Q' = Q · Λ is Haar distributed. # NOTE: Λ is a diagonal matrix, represented as a vector # of the diagonal entries. Thus, the matrix dot product # is represented nicely by the NumPy broadcasting of # the *scalar* multiplication. In particular, # Q · Λ = Q_ij Λ_jk = Q_ij δ_jk λ_k = Q_ij λ_j. # As NumPy arrays, Q has shape (N, N) and # Lambda has shape (N, ), such that the broadcasting # represents precisely Q_ij λ_j. U = Qobj(Q * Lambda) U.dims = dims return U
[docs]def rand_ket(N=0, density=1, dims=None, seed=None): """Creates a random Nx1 sparse ket vector. Parameters ---------- N : int Number of rows for output quantum operator. If None or 0, N is deduced from dims. density : float Density between [0,1] of output ket state. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[1]]. Returns ------- oper : qobj Nx1 ket state quantum operator. """ if seed is not None: np.random.seed(seed=seed) if N and dims: _check_dims(dims, N, 1) elif dims: N = prod(dims[0]) _check_dims(dims, N, 1) else: dims = [[N],[1]] X = sp.rand(N, 1, density, format='csr') while X.nnz == 0: # ensure that the ket is not all zeros. X = sp.rand(N, 1, density+1/N, format='csr') X.data = X.data - 0.5 Y = X.copy() Y.data = 1.0j * (np.random.random(len(X.data)) - 0.5) X = X + Y X.sort_indices() X = Qobj(X) return Qobj(X / X.norm(), dims=dims)
[docs]def rand_ket_haar(N=2, dims=None, seed=None): """ Returns a Haar random pure state of dimension ``dim`` by applying a Haar random unitary to a fixed pure state. Parameters ---------- N : int Dimension of the state vector to be returned. If None or 0, N is deduced from dims. dims : list of ints, or None Dimensions of the resultant quantum object. If None, [[N],[1]] is used. Returns ------- psi : Qobj A random state vector drawn from the Haar measure. """ if N and dims: _check_dims(dims, N, 1) elif dims: N = prod(dims[0]) _check_dims(dims, N, 1) else: dims = [[N],[1]] psi = rand_unitary_haar(N, seed=seed) * basis(N, 0) psi.dims = dims return psi
[docs]def rand_dm(N, density=0.75, pure=False, dims=None, seed=None): """Creates a random NxN density matrix. Parameters ---------- N : int, ndarray, list If int, then shape of output operator. If list/ndarray then eigenvalues of generated density matrix. density : float Density between [0,1] of output density matrix. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- oper : qobj NxN density matrix quantum operator. Notes ----- For small density matrices., choosing a low density will result in an error as no diagonal elements will be generated such that :math:`Tr(\\rho)=1`. """ if isinstance(N, (np.ndarray, list)): if np.abs(np.sum(N)-1.0) > 1e-15: raise ValueError('Eigenvalues of a density matrix ' 'must sum to one.') H = sp.diags(N, 0, dtype=complex, format='csr') N = len(N) if dims: _check_dims(dims, N, N) nvals = N**2*density H = rand_jacobi_rotation(H, seed=seed) while H.nnz < 0.95*nvals: H = rand_jacobi_rotation(H) H.sort_indices() elif isinstance(N, (int, np.int32, np.int64)): if dims: _check_dims(dims, N, N) if pure: dm_density = sqrt(density) psi = rand_ket(N, dm_density) H = psi * psi.dag() H.data.sort_indices() else: non_zero = 0 tries = 0 while non_zero == 0 and tries < 10: H = rand_herm(N, density, seed=seed) H = H.dag() * H non_zero = H.tr() tries += 1 if tries >= 10: raise ValueError( "Requested density is too low to generate density matrix.") H = H / H.tr() H.data.sort_indices() else: raise TypeError('Input N must be an integer or array_like.') if dims: return Qobj(H, dims=dims) else: return Qobj(H)
[docs]def rand_dm_ginibre(N=2, rank=None, dims=None, seed=None): """ Returns a Ginibre random density operator of dimension ``dim`` and rank ``rank`` by using the algorithm of [BCSZ08]_. If ``rank`` is `None`, a full-rank (Hilbert-Schmidt ensemble) random density operator will be returned. Parameters ---------- N : int Dimension of the density operator to be returned. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. rank : int or None Rank of the sampled density operator. If None, a full-rank density operator is generated. Returns ------- rho : Qobj An N × N density operator sampled from the Ginibre or Hilbert-Schmidt distribution. """ if rank is None: rank = N if rank > N: raise ValueError("Rank cannot exceed dimension.") X = randnz((N, rank), norm='ginibre', seed=seed) rho = np.dot(X, X.T.conj()) rho /= np.trace(rho) return Qobj(rho, dims=dims)
[docs]def rand_dm_hs(N=2, dims=None, seed=None): """ Returns a Hilbert-Schmidt random density operator of dimension ``dim`` and rank ``rank`` by using the algorithm of [BCSZ08]_. Parameters ---------- N : int Dimension of the density operator to be returned. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- rho : Qobj A dim × dim density operator sampled from the Ginibre or Hilbert-Schmidt distribution. """ return rand_dm_ginibre(N, rank=None, dims=dims, seed=seed)
def rand_kraus_map(N, dims=None, seed=None): """ Creates a random CPTP map on an N-dimensional Hilbert space in Kraus form. Parameters ---------- N : int Length of input/output density matrix. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- oper_list : list of qobj N^2 x N x N qobj operators. """ if dims: _check_dims(dims, N, N) # Random unitary (Stinespring Dilation) big_unitary = rand_unitary(N ** 3, seed=seed).data.todense() orthog_cols = np.array(big_unitary[:, :N]) oper_list = np.reshape(orthog_cols, (N ** 2, N, N)) return list(map(lambda x: Qobj(inpt=x, dims=dims), oper_list))
[docs]def rand_super(N=5, dims=None, seed=None): """ Returns a randomly drawn superoperator acting on operators acting on N dimensions. Parameters ---------- N : int Square root of the dimension of the superoperator to be returned. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[[N],[N]], [[N],[N]]]. """ if dims is not None: # TODO: check! pass else: dims = [[[N],[N]], [[N],[N]]] H = rand_herm(N, seed=seed) S = propagator(H, np.random.rand(), [ create(N), destroy(N), jmat(float(N - 1) / 2.0, 'z') ]) S.dims = dims return S
[docs]def rand_super_bcsz(N=2, enforce_tp=True, rank=None, dims=None, seed=None): """ Returns a random superoperator drawn from the Bruzda et al ensemble for CPTP maps [BCSZ08]_. Note that due to finite numerical precision, for ranks less than full-rank, zero eigenvalues may become slightly negative, such that the returned operator is not actually completely positive. Parameters ---------- N : int Square root of the dimension of the superoperator to be returned. enforce_tp : bool If True, the trace-preserving condition of [BCSZ08]_ is enforced; otherwise only complete positivity is enforced. rank : int or None Rank of the sampled superoperator. If None, a full-rank superoperator is generated. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[[N],[N]], [[N],[N]]]. Returns ------- rho : Qobj A superoperator acting on vectorized dim × dim density operators, sampled from the BCSZ distribution. """ if dims is not None: # TODO: check! pass else: dims = [[[N],[N]], [[N],[N]]] if rank is None: rank = N**2 if rank > N**2: raise ValueError("Rank cannot exceed superoperator dimension.") # We use mainly dense matrices here for speed in low # dimensions. In the future, it would likely be better to switch off # between sparse and dense matrices as the dimension grows. # We start with a Ginibre uniform matrix X of the appropriate rank, # and use it to construct a positive semidefinite matrix X X⁺. X = randnz((N**2, rank), norm='ginibre', seed=seed) # Precompute X X⁺, as we'll need it in two different places. XXdag = np.dot(X, X.T.conj()) if enforce_tp: # We do the partial trace over the first index by using dense reshape # operations, so that we can avoid bouncing to a sparse representation # and back. Y = np.einsum('ijik->jk', XXdag.reshape((N, N, N, N))) # Now we have the matrix 𝟙 ⊗ Y^{-1/2}, which we can find by doing # the square root and the inverse separately. As a possible # improvement, iterative methods exist to find inverse square root # matrices directly, as this is important in statistics. Z = np.kron( np.eye(N), sqrtm(la.inv(Y)) ) # Finally, we dot everything together and pack it into a Qobj, # marking the dimensions as that of a type=super (that is, # with left and right compound indices, each representing # left and right indices on the underlying Hilbert space). D = Qobj(np.dot(Z, np.dot(XXdag, Z))) else: D = N * Qobj(XXdag / np.trace(XXdag)) D.dims = [ # Left dims [[N], [N]], # Right dims [[N], [N]] ] # Since [BCSZ08] gives a row-stacking Choi matrix, but QuTiP # expects a column-stacking Choi matrix, we must permute the indices. D = D.permute([[1], [0]]) D.dims = dims # Mark that we've made a Choi matrix. D.superrep = 'choi' return sr.to_super(D)
[docs]def rand_stochastic(N, density=0.75, kind='left', dims=None, seed=None): """Generates a random stochastic matrix. Parameters ---------- N : int Dimension of matrix. density : float Density between [0,1] of output density matrix. kind : str (Default = 'left') Generate 'left' or 'right' stochastic matrix. dims : list Dimensions of quantum object. Used for specifying tensor structure. Default is dims=[[N],[N]]. Returns ------- oper : qobj Quantum operator form of stochastic matrix. """ if seed is not None: np.random.seed(seed=seed) if dims: _check_dims(dims, N, N) num_elems = max([np.int(np.ceil(N*(N+1)*density)/2), N]) data = np.random.rand(num_elems) # Ensure an element on every row and column row_idx = np.hstack([np.random.permutation(N), np.random.choice(N, num_elems-N)]) col_idx = np.hstack([np.random.permutation(N), np.random.choice(N, num_elems-N)]) M = sp.coo_matrix((data, (row_idx, col_idx)), dtype=float, shape=(N,N)).tocsr() M = 0.5 * (M + M.conj().transpose()) num_rows = M.indptr.shape[0]-1 for row in range(num_rows): row_start = M.indptr[row] row_end = M.indptr[row+1] row_sum = np.sum(M.data[row_start:row_end]) M.data[row_start:row_end] /= row_sum if kind=='left': M = M.transpose() if dims: return Qobj(M, dims=dims, shape=[N, N]) else: return Qobj(M)
def _check_dims(dims, N1, N2): if len(dims) != 2: raise Exception("Qobj dimensions must be list of length 2.") if (not isinstance(dims[0], list)) or (not isinstance(dims[1], list)): raise TypeError( "Qobj dimension components must be lists. i.e. dims=[[N],[N]]") if np.prod(dims[0]) != N1 or np.prod(dims[1]) != N2: raise ValueError("Qobj dimensions must match matrix shape.") if len(dims[0]) != len(dims[1]): raise TypeError("Qobj dimension components must have same length.") # TRAILING IMPORTS # qutip.propagator depends on rand_dm, so we need to put this import last. from qutip.propagator import propagator