Adding random unitary (#1857)

* Adding random unitary that does not use spicy

* Changhegog

CHANGELOG.rst

@@ -40,6 +40,7 @@ Added
 Changed
 -------
 
+- Change random_state to take in dim over number of qubits (#1857)
 - The ``Exception`` subclasses have been moved to an ``.exceptions`` module
   within each package (for example, ``qiskit.exceptions.QiskitError``) (#1600).
 - The ``QiskitTestCase`` and testing utilities are now included as part of
@@ -80,6 +81,7 @@ Deprecated
 Fixed
 -----
 
+- Fixed #829 by removing dependence on scipy unitary_group (#1857).
 - Fixed a bug with measurement sampling optimization in BasicAer
   qasm_simulator (#1624).
 - Fixed a bug where barriers didn't plot over all qubits when using matplotlib (#1718).

qiskit/quantum_info/states/_states.py

@@ -38,13 +38,13 @@ def basis_state(str_state, num):
         raise QiskitError('size of bitstring is greater than num.')
 
 
-def random_state(num, seed=None):
+def random_state(dim, seed=None):
     """
     Return a random quantum state from the uniform (Haar) measure on
     state space.
 
     Args:
-        num (int): the number of qubits
+        dim (int): the dim of the state spaxe
         seed (int): Optional. To set a random seed.
 
     Returns:
@@ -53,11 +53,11 @@ def random_state(num, seed=None):
     if seed is not None:
         np.random.seed(seed)
     # Random array over interval (0, 1]
-    x = np.random.random(1 << num)
+    x = np.random.random(dim)
     x += x == 0
     x = -np.log(x)
     sumx = sum(x)
-    phases = np.random.random(1 << num)*2.0*np.pi
+    phases = np.random.random(dim)*2.0*np.pi
     return np.sqrt(x/sumx)*np.exp(1j*phases)
 
 

qiskit/tools/qi/qi.py

@@ -20,12 +20,12 @@ import warnings
 
 import numpy as np
 import scipy.linalg as la
-from scipy.stats import unitary_group
 
 
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info import pauli_group
 from qiskit.quantum_info import purity as new_purity
+from qiskit.quantum_info import random_state
 
 
 ###############################################################
@@ -336,19 +336,32 @@ def outer(vector1, vector2=None):
 # Random Matrices.
 ###############################################################
 
-def random_unitary_matrix(length, seed=None):
+def random_unitary_matrix(dim, seed=None):
     """
-    Return a random unitary ndarray.
+    Return a random unitary ndarray over num qubits.
 
     Args:
-        length (int): the length of the returned unitary.
+        dim (int): the dim of the state space.
         seed (int): Optional. To set a random seed.
     Returns:
-        ndarray: U (length, length) unitary ndarray.
+        ndarray: U (2**num, 2**num) unitary ndarray.
     """
-    if seed is not None:
-        np.random.seed(seed)
-    return unitary_group.rvs(length)
+    unitary = np.zeros([dim, dim], dtype=complex)
+    for j in range(dim):
+        if j == 0:
+            a = random_state(dim, seed)
+        else:
+            a = random_state(dim)
+        unitary[:, j] = np.copy(a)
+        # Grahm-Schmidt Orthogonalize
+        i = j-1
+        while i >= 0:
+            dc = np.vdot(unitary[:, i], a)
+            unitary[:, j] = unitary[:, j]-dc*unitary[:, i]
+            i = i - 1
+        # normalize
+        unitary[:, j] = unitary[:, j] * (1.0 / np.sqrt(np.vdot(unitary[:, j], unitary[:, j])))
+    return unitary
 
 
 def random_density_matrix(length, rank=None, method='Hilbert-Schmidt', seed=None):

requirements.txt

@@ -6,5 +6,5 @@ numpy>=1.13
 pillow>=4.2.1
 ply>=3.10
 psutil>=5
-scipy>=0.19,!=0.19.1
+scipy>=0.19
 sympy>=1.3

test/python/quantum_info/test_states.py

@@ -84,12 +84,12 @@ class TestStates(QiskitTestCase):
         E_P0_last = 0
         for ii in range(number):
             state = basis_state(bin(3)[2:].zfill(3), 3)
-            E_P0 = (E_P0_last*ii)/(ii+1)+state_fidelity(state, random_state(3, seed=ii))/(ii+1)
+            E_P0 = (E_P0_last*ii)/(ii+1)+state_fidelity(state, random_state(2**3, seed=ii))/(ii+1)
             E_P0_last = E_P0
         self.assertAlmostEqual(E_P0, 1/8, places=2)
 
     def test_random_state_circuit(self):
-        state = random_state(3, seed=40)
+        state = random_state(2**3, seed=40)
         q = QuantumRegister(3)
         qc = QuantumCircuit(q)
         qc.initialize(state, [q[0], q[1], q[2]])