From 364debc8dc3e05dd7ab1f2a813ab5d779dd5249d Mon Sep 17 00:00:00 2001
From: "Kevin J. Sung" <kevjsung@umich.edu>
Date: Fri, 22 Sep 2023 06:10:59 -0400
Subject: [PATCH] fix random statevector distribution (#10866)

* fix random statevector distribution

* add citation and release note

* mention Haar measure in doc and use link in release note

* cite better reference

* Format reference

---------

Co-authored-by: Julien Gacon <gaconju@gmail.com>
---
 qiskit/quantum_info/states/random.py          | 19 ++++++++++---------
 ...x-random-statevector-a8dbf991cbbdee8e.yaml |  4 ++++
 2 files changed, 14 insertions(+), 9 deletions(-)
 create mode 100644 releasenotes/notes/fix-random-statevector-a8dbf991cbbdee8e.yaml

diff --git a/qiskit/quantum_info/states/random.py b/qiskit/quantum_info/states/random.py
index 5065b33b54d4..76dbd4e4760e 100644
--- a/qiskit/quantum_info/states/random.py
+++ b/qiskit/quantum_info/states/random.py
@@ -31,7 +31,8 @@ def random_statevector(
 ) -> Statevector:
     """Generator a random Statevector.
 
-    The statevector is sampled from the uniform (Haar) measure.
+    The statevector is sampled from the uniform distribution. This is the measure
+    induced by the Haar measure on unitary matrices.
 
     Args:
         dims (int or tuple): the dimensions of the state.
@@ -40,6 +41,10 @@ def random_statevector(
 
     Returns:
         Statevector: the random statevector.
+
+    Reference:
+        K. Zyczkowski and H. Sommers (2001), "Induced measures in the space of mixed quantum states",
+        `J. Phys. A: Math. Gen. 34 7111 <https://arxiv.org/abs/quant-ph/0012101>`__.
     """
     if seed is None:
         rng = np.random.default_rng()
@@ -49,14 +54,10 @@ def random_statevector(
         rng = default_rng(seed)
 
     dim = np.prod(dims)
-
-    # Random array over interval (0, 1]
-    x = rng.random(dim)
-    x += x == 0
-    x = -np.log(x)
-    sumx = sum(x)
-    phases = rng.random(dim) * 2.0 * np.pi
-    return Statevector(np.sqrt(x / sumx) * np.exp(1j * phases), dims=dims)
+    vec = rng.standard_normal(dim).astype(complex)
+    vec += 1j * rng.standard_normal(dim)
+    vec /= np.linalg.norm(vec)
+    return Statevector(vec, dims=dims)
 
 
 def random_density_matrix(
diff --git a/releasenotes/notes/fix-random-statevector-a8dbf991cbbdee8e.yaml b/releasenotes/notes/fix-random-statevector-a8dbf991cbbdee8e.yaml
new file mode 100644
index 000000000000..75e78790cefc
--- /dev/null
+++ b/releasenotes/notes/fix-random-statevector-a8dbf991cbbdee8e.yaml
@@ -0,0 +1,4 @@
+---
+fixes:
+  - |
+    Fixes the implementation of :func:`.random_statevector` so that it samples from the uniform distribution.