{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# This cell is added by sphinx-gallery\n# It can be customized to whatever you like\n%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Understanding the Haar Measure\n==============================\n\n::: {.meta}\n:property=\\\"og:description\\\": Learn all about the Haar measure and how\nto randomly sample quantum states.\n\n:property=\\\"og:image\\\":\n\n:::\n\n::: {.related}\ntutorial\\_unitary\\_designs Unitary designs quantum\\_volume Quantum\nvolume qsim\\_beyond\\_classical Beyond classical computing with qsim\ntutorial\\_barren\\_plateaus Barren plateaus\n:::\n\n*Author: PennyLane dev team. Posted: 22 March 2021. Last updated: 22\nMarch 2021.*\n\nIf you\\'ve ever dug into the literature about random quantum circuits,\nvariational ansatz structure, or anything related to the structure and\nproperties of unitary operations, you\\'ve likely come across a statement\nlike the following: \\\"Assume that $U$ is sampled uniformly at random\nfrom the Haar measure\\\". In this demo, we\\'re going to unravel this\ncryptic statement and take an in-depth look at what it means. You\\'ll\ngain an understanding of the general concept of *measure*, the Haar\nmeasure and its special properties, and you\\'ll learn how to sample from\nit using tools available in PennyLane and other scientific computing\nframeworks. By the end of this demo, you\\'ll be able to include that\nimportant statement in your own work with confidence!\n\n::: {.note}\n::: {.admonition-title}\nNote\n:::\n\nTo get the most out of this demo, it is helpful if you are familiar with\n[integration of multi-dimensional\nfunctions](https://en.wikipedia.org/wiki/Multiple_integral), the [Bloch\nsphere](https://en.wikipedia.org/wiki/Bloch_sphere), and the conceptual\nideas behind\n[decompositions](https://en.wikipedia.org/wiki/Matrix_decomposition) and\nfactorizations of unitary matrices (see, e.g., 4.5.1 and 4.5.2 of).\n:::\n\nMeasure\n-------\n\n[Measure theory](https://en.wikipedia.org/wiki/Measure_(mathematics)) is\na branch of mathematics that studies things that are\nmeasurable\\-\\--think length, area, or volume, but generalized to\nmathematical spaces and even higher dimensions. Loosely, the measure\ntells you about how \\\"stuff\\\" is distributed and concentrated in a\nmathematical set or space. An intuitive way to understand measure is to\nthink about a sphere. An arbitrary point on a sphere can be parametrized\nby three numbers\\-\\--depending on what you\\'re doing, you may use\nCartesian coordinates $(x, y, z)$, or it may be more convenient to use\nspherical coordinates $(\\rho, \\phi, \\theta)$.\n\nSuppose you wanted to compute the volume of a solid sphere with radius\n$r$. This can be done by integrating over the three coordinates\n$\\rho, \\phi$, and $\\theta$. Your first thought here may be to simply\nintegrate each parameter over its full range, like so:\n\n$$V = \\int_0^{r} \\int_0^{2\\pi} \\int_0^{\\pi} d\\rho~ d\\phi~ d\\theta = 2\\pi^2 r$$\n\nBut, we know that the volume of a sphere of radius $r$ is\n$\\frac{4}{3}\\pi r^3$, so what we got from this integral is clearly\nwrong! Taking the integral naively like this doesn\\'t take into account\nthe structure of the sphere with respect to the parameters. For example,\nconsider two small, infinitesimal elements of area with the same\ndifference in $\\theta$ and $\\phi$, but at different values of $\\theta$:\n\n![\\|](/demonstrations/haar_measure/spherical_int_dtheta.png){.align-center\nwidth=\"50.0%\"}\n\nEven though the differences $d\\theta$ and $d\\phi$ themselves are the\nsame, there is way more \\\"stuff\\\" near the equator of the sphere than\nthere is near the poles. We must take into account the value of $\\theta$\nwhen computing the integral! Specifically, we multiply by the function\n$\\sin\\theta$\\-\\--the properties of the $\\sin$ function mean that the\nmost weight will occur around the equator where $\\theta=\\pi/2$, and the\nleast weight near the poles where $\\theta=0$ and $\\theta=\\pi$.\n\nSimilar care must be taken for $\\rho$. The contribution to volume of\nparts of the sphere with a large $\\rho$ is far more than for a small\n$\\rho$\\-\\--we should expect the contribution to be proportional to\n$\\rho^2$, given that the surface area of a sphere of radius $r$ is\n$4\\pi r^2$.\n\nOn the other hand, for a fixed $\\rho$ and $\\theta$, the length of the\n$d\\phi$ is the same all around the circle. If put all these facts\ntogether, we find that the actual expression for the integral should\nlook like this:\n\n$$V = \\int_0^r \\int_0^{2\\pi} \\int_0^{\\pi} \\rho^2 \\sin \\theta~ d\\rho~ d\\phi~\nd\\theta = \\frac{4}{3}\\pi r^3$$\n\nThese extra terms that we had to add to the integral, $\\rho^2 \\sin\n\\theta$, constitute the *measure*. The measure weights portions of the\nsphere differently depending on where they are in the space. While we\nneed to know the measure to properly integrate over the sphere,\nknowledge of the measure also gives us the means to perform another\nimportant task, that of sampling points in the space uniformly at\nrandom. We can\\'t simply sample each parameter from the uniform\ndistribution over its domain\\-\\--as we experienced already, this\ndoesn\\'t take into account how the sphere is spread out over space. The\nmeasure describes the distribution of each parameter and gives a recipe\nfor sampling them in order to obtain something properly uniform.\n\nThe Haar measure\n----------------\n\nOperations in quantum computing are described by unitary matrices.\nUnitary matrices, like points on a sphere, can be expressed in terms of\na fixed set of coordinates, or parameters. For example, the most general\nsingle-qubit rotation implemented in PennyLane\n(`~.pennylane.Rot`{.interpreted-text role=\"class\"}) is expressed in\nterms of three parameters like so,\n\n$$\\begin{aligned}\nU(\\phi, \\theta, \\omega) = \\begin{pmatrix} e^{-i(\\phi + \\omega)/2}\n \\cos(\\theta/2) & -e^{i(\\phi - \\omega)/2} \\sin(\\theta/2)\n \\\\ e^{-i(\\phi - \\omega)/2} \\sin(\\theta/2) & e^{i(\\phi +\n \\omega)/2} \\cos(\\theta/2) \\end{pmatrix}.\n\\end{aligned}$$\n\nFor every dimension $N$, the unitary matrices of size $N \\times N$\nconstitute the *unitary group* $U(N)$. We can perform operations on\nelements of this group, such as apply functions to them, integrate over\nthem, or sample uniformly over them, just as we can do to points on a\nsphere. When we do such tasks with respect to the sphere, we have to add\nthe measure in order to properly weight the different regions of space.\nThe *Haar measure* provides the analogous terms we need for working with\nthe unitary group.\n\nFor an $N$-dimensional system, the Haar measure, often denoted by\n$\\mu_N$, tells us how to weight the elements of $U(N)$. For example,\nsuppose $f$ is a function that acts on elements of $U(N)$, and we would\nlike to take its integral over the group. We must write this integral\nwith respect to the Haar measure, like so:\n\n$$\\int_{V \\in U(N)} f(V) d\\mu_N(V).$$\n\nAs with the measure term of the sphere, $d\\mu_N$ itself can be broken\ndown into components depending on individual parameters. While the Haar\nmeasure can be defined for every dimension $N$, the mathematical form\ngets quite hairy for larger dimensions\\-\\--in general, an\n$N$-dimensional unitary requires at least $N^2 - 1$ parameters, which is\na lot to keep track of! Therefore we\\'ll start with the case of a single\nqubit $(N=2)$, then show how things generalize.\n\n### Single-qubit Haar measure\n\nThe single-qubit case provides a particularly nice entry point because\nwe can continue our comparison to spheres by visualizing single-qubit\nstates on the Bloch sphere. As expressed above, the measure provides a\nrecipe for sampling elements of the unitary group in a properly uniform\nmanner, given the structure of the group. One useful consequence of this\nis that it provides a method to sample quantum *states* uniformly at\nrandom\\-\\--we simply generate Haar-random unitaries, and apply them to a\nfixed basis state such as $\\vert 0\\rangle$.\n\nWe\\'ll see how this works in good time. First, we\\'ll take a look at\nwhat happens when we ignore the measure and do things *wrong*. Suppose\nwe sample quantum states by applying unitaries obtained by the\nparametrization above, but sample the angles $\\omega, \\phi$, and\n$\\theta$ from the flat uniform distribution between $[0, 2\\pi)$ (fun\nfact: there is a measure implicit in this kind of sampling too! It just\nhas a constant value, because each point is equally likely to be\nsampled).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import pennylane as qml\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n# set the random seed\nnp.random.seed(42)\n\n# Use the mixed state simulator to save some steps in plotting later\ndev = qml.device('default.mixed', wires=1)\n\n@qml.qnode(dev)\ndef not_a_haar_random_unitary():\n # Sample all parameters from their flat uniform distribution\n phi, theta, omega = 2 * np.pi * np.random.uniform(size=3)\n qml.Rot(phi, theta, omega, wires=0)\n return qml.state()\n\nnum_samples = 2021\n\nnot_haar_samples = [not_a_haar_random_unitary() for _ in range(num_samples)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to plot these on the Bloch sphere, we\\'ll need to do one more\nstep, and convert the quantum states into Bloch vectors.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X = np.array([[0, 1], [1, 0]])\nY = np.array([[0, -1j], [1j, 0]])\nZ = np.array([[1, 0], [0, -1]])\n\n# Used the mixed state simulator so we could have the density matrix for this part!\ndef convert_to_bloch_vector(rho):\n \"\"\"Convert a density matrix to a Bloch vector.\"\"\"\n ax = np.trace(np.dot(rho, X)).real\n ay = np.trace(np.dot(rho, Y)).real\n az = np.trace(np.dot(rho, Z)).real\n return [ax, ay, az]\n\nnot_haar_bloch_vectors = np.array([convert_to_bloch_vector(s) for s in not_haar_samples])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this done, let\\'s find out where our \\\"uniformly random\\\" states\nended up:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def plot_bloch_sphere(bloch_vectors):\n \"\"\" Helper function to plot vectors on a sphere.\"\"\"\n fig = plt.figure(figsize=(6, 6))\n ax = fig.add_subplot(111, projection='3d')\n fig.subplots_adjust(left=0, right=1, bottom=0, top=1)\n\n ax.grid(False)\n ax.set_axis_off()\n ax.view_init(30, 45)\n ax.dist = 7\n\n # Draw the axes (source: https://github.com/matplotlib/matplotlib/issues/13575)\n x, y, z = np.array([[-1.5,0,0], [0,-1.5,0], [0,0,-1.5]])\n u, v, w = np.array([[3,0,0], [0,3,0], [0,0,3]])\n ax.quiver(x, y, z, u, v, w, arrow_length_ratio=0.05, color=\"black\", linewidth=0.5)\n\n ax.text(0, 0, 1.7, r\"|0\u27e9\", color=\"black\", fontsize=16)\n ax.text(0, 0, -1.9, r\"|1\u27e9\", color=\"black\", fontsize=16)\n ax.text(1.9, 0, 0, r\"|+\u27e9\", color=\"black\", fontsize=16)\n ax.text(-1.7, 0, 0, r\"|\u2013\u27e9\", color=\"black\", fontsize=16)\n ax.text(0, 1.7, 0, r\"|i+\u27e9\", color=\"black\", fontsize=16)\n ax.text(0,-1.9, 0, r\"|i\u2013\u27e9\", color=\"black\", fontsize=16)\n\n ax.scatter(\n bloch_vectors[:,0], bloch_vectors[:,1], bloch_vectors[:, 2], c='#e29d9e', alpha=0.3\n )\n\nplot_bloch_sphere(not_haar_bloch_vectors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see from this plot that even though our parameters were sampled\nfrom a uniform distribution, there is a noticeable amount of clustering\naround the poles of the sphere. Despite the input parameters being\nuniform, the output is very much *not* uniform. Just like the regular\nsphere, the measure is larger near the equator, and if we just sample\nuniformly, we won\\'t end up populating that area as much. To take that\ninto account we will need to sample from the proper Haar measure, and\nweight the different parameters appropriately.\n\nFor a single qubit, the Haar measure looks much like the case of a\nsphere, minus the radial component. Intuitively, all qubit state vectors\nhave length 1, so it makes sense that this wouldn\\'t play a role here.\nThe parameter that we will have to weight differently is $\\theta$, and\nin fact the adjustment in measure is identical to that we had to do with\nthe polar axis of the sphere, i.e., $\\sin \\theta$. In order to sample\nthe $\\theta$ uniformly at random in this context, we must sample from\nthe distribution $\\hbox{Pr}(\\theta) = \\sin \\theta$. We can accomplish\nthis by setting up a custom probability distribution with\n[rv\\_continuous](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html#scipy.stats.rv_continuous)\nin `scipy`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy.stats import rv_continuous\n\nclass sin_prob_dist(rv_continuous):\n def _pdf(self, theta):\n # The 0.5 is so that the distribution is normalized\n return 0.5 * np.sin(theta)\n\n# Samples of theta should be drawn from between 0 and pi\nsin_sampler = sin_prob_dist(a=0, b=np.pi)\n\n@qml.qnode(dev)\ndef haar_random_unitary():\n phi, omega = 2 * np.pi * np.random.uniform(size=2) # Sample phi and omega as normal\n theta = sin_sampler.rvs(size=1) # Sample theta from our new distribution\n qml.Rot(phi, theta, omega, wires=0)\n return qml.state()\n\nhaar_samples = [haar_random_unitary() for _ in range(num_samples)]\nhaar_bloch_vectors = np.array([convert_to_bloch_vector(s) for s in haar_samples])\n\nplot_bloch_sphere(haar_bloch_vectors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that when we use the correct measure, our qubit states are now\nmuch better distributed over the sphere. Putting this information\ntogether, we can now write the explicit form for the single-qubit Haar\nmeasure:\n\n$$d\\mu_2 = \\sin \\theta d\\theta \\cdot d\\omega \\cdot d\\phi.$$\n\nShow me more math!\n==================\n\nWhile we can easily visualize the single-qubit case, this is no longer\npossible when we increase the number of qubits. Regardless, we can still\nobtain a mathematical expression for the Haar measure in arbitrary\ndimensions. In the previous section, we expressed the Haar measure in\nterms of a set of parameters that can be used to specify the unitary\ngroup $U(2)$. Such a parametrization is not unique, and in fact there\nare multiple ways to *factorize*, or decompose an $N$-dimensional\nunitary operation into a set of parameters.\n\nMany of these parametrizations come to us from the study of photonics.\nHere, arbitrary operations are broken down into elementary operations\ninvolving only a few parameters which correspond directly to parameters\nof the physical apparatus used to implement them (beamsplitters and\nphase shifts). Rather than qubits, such operations act on modes, or\n*qumodes*. They are expressed as elements of the $N$-dimensional\n[special unitary\ngroup](https://en.wikipedia.org/wiki/Special_unitary_group). This group,\nwritten as $SU(N)$, is the continuous group consisting of all\n$N \\times N$ unitary operations with determinant 1 (essentially like\n$U(N)$, minus a potential global phase).\n\n::: {.note}\n::: {.admonition-title}\nNote\n:::\n\nElements of $SU(N)$ and $U(N)$ can still be considered as multi-qubit\noperations in the cases where $N$ is a power of 2, but they must be\ntranslated from continuous-variable operations into qubit operations.\n(In PennyLane, this can be done by feeding the unitaries to the\n`~.pennylane.QubitUnitary`{.interpreted-text role=\"class\"} operation\ndirectly. Alternatively, one can use *quantum compilation* to express\nthe operations as a sequence of elementary gates such as Pauli rotations\nand CNOTs.)\n:::\n\n::: {.admonition}\nTip\n\nIf you haven\\'t had many opportunities to work in terms of qumodes, the\n[Strawberry Fields\ndocumentation](https://strawberryfields.ai/photonics/concepts/photonics.html)\nis a good starting point.\n:::\n\nFor example, we saw already above that for $N=2$, we can write\n\n$$\\begin{aligned}\nU(\\phi, \\theta, \\omega) = \\begin{pmatrix} e^{-i(\\phi + \\omega)/2}\n \\cos(\\theta/2) & -e^{i(\\phi - \\omega)/2} \\sin(\\theta/2)\n \\\\ e^{-i(\\phi - \\omega)/2} \\sin(\\theta/2) & e^{i(\\phi +\n \\omega)/2} \\cos(\\theta/2) \\end{pmatrix}.\n\\end{aligned}$$\n\nThis unitary can be factorized as follows:\n\n$$\\begin{aligned}\nU(\\phi, \\theta, \\omega) =\n \\begin{pmatrix}\n e^{-i\\omega/2} & 0 \\\\ 0 & e^{i\\omega/2}\n \\end{pmatrix}\n \\begin{pmatrix}\n \\cos(\\theta/2) & -\\sin(\\theta/2) \\\\ \\sin(\\theta/2) & \\cos(\\theta/2)\n \\end{pmatrix}\n \\begin{pmatrix}\n e^{-i\\phi/2} & 0 \\\\ 0 & e^{i\\phi/2}\n \\end{pmatrix}\n\\end{aligned}$$\n\nThe middle operation is a beamsplitter; the other two operations are\nphase shifts. We saw earlier that for $N=2$, $d\\mu_2 = \\sin\\theta\nd\\theta d\\omega d\\phi$\\-\\--note how the parameter in the beamsplitter\ncontributes to the measure in a different way than those of the phase\nshifts. As mentioned above, for larger values of $N$ there are multiple\nways to decompose the unitary. Such decompositions rewrite elements in\n$SU(N)$ acting on $N$ modes as a sequence of operations acting only on 2\nmodes, $SU(2)$, and single-mode phase shifts. Shown below are three\nexamples,,:\n\n![](/demonstrations/haar_measure/unitaries.png){.align-center\nwidth=\"95.0%\"}\n\nIn these graphics, every wire is a different mode. Every box represents\nan operation on one or more modes, and the number in the box indicates\nthe number of parameters. The boxes containing a `1` are simply phase\nshifts on individual modes. The blocks containing a `3` are $SU(2)$\ntransforms with 3 parameters, such as the $U(\\phi, \\theta, \\omega)$\nabove. Those containing a `2` are $SU(2)$ transforms on pairs of modes\nwith 2 parameters, similar to the 3-parameter ones but with\n$\\omega = \\phi$.\n\nAlthough the decompositions all produce the same set of operations,\ntheir structure and parametrization may have consequences in practice.\nThe first has a particularly convenient form that leads to a recursive\ndefinition of the Haar measure. The decomposition is formulated\nrecursively such that an $SU(N)$ operation can be implemented by\nsandwiching an $SU(2)$ transformation between two $SU(N-1)$\ntransformations, like so:\n\n| \n\n![](/demonstrations/haar_measure/sun.svg){.align-center width=\"80.0%\"}\n\n| \n\nThe Haar measure is then constructed recursively as a product of 3\nterms. The first term depends on the parameters in the first $SU(N-1)$\ntransformation; the second depends on the parameters in the lone $SU(2)$\ntransformation; and the third term depends on the parameters in the\nother $SU(N-1)$ transformation.\n\n$SU(2)$ is the \\\"base case\\\" of the recursion\\-\\--we simply have the\nHaar measure as expressed above.\n\n| \n\n![](/demonstrations/haar_measure/su2_haar.svg){.align-center\nwidth=\"25.0%\"}\n\n| \n\nMoving on up, we can write elements of $SU(3)$ as a sequence of three\n$SU(2)$ transformations. The Haar measure $d\\mu_3$ then consists of two\ncopies of $d\\mu_2$, with an extra term in between to take into account\nthe middle transformation.\n\n| \n\n![](/demonstrations/haar_measure/su3_haar.svg){.align-center\nwidth=\"80.0%\"}\n\n| \n\nFor $SU(4)$ and upwards, the form changes slightly, but still follows\nthe pattern of two copies of $d\\mu_{N-1}$ with a term in between.\n\n| \n\n![](/demonstrations/haar_measure/su4_premerge.svg){.align-center\nwidth=\"90.0%\"}\n\n| \n\nFor larger systems, however, the recursive composition allows for some\nof the $SU(2)$ transformations on the lower modes to be grouped. We can\ntake advantage of this and aggregate some of the parameters:\n\n| \n\n![](/demonstrations/haar_measure/su4_triangle_merge.svg){.align-center\nwidth=\"100.0%\"}\n\n| \n\nThis leads to one copy of $d\\mu_{N-1}$, which we\\'ll denote as\n$d\\mu_{N-1}^\\prime$, containing only a portion of the full set of terms\n(as detailed in, this is called a *coset measure*).\n\n| \n\n![](/demonstrations/haar_measure/su4_haar.svg){.align-center\nwidth=\"100.0%\"}\n\n| \n\nPutting everything together, we have that\n\n$$d\\mu_N = d\\mu_{N-1}^\\prime \\times \\sin \\theta_{N-1}\n\\sin^{2(N-2)}\\left(\\frac{\\theta_{N-1}}{2}\\right) d\\theta_{N-1} d\\omega_{N-1} \\times d\\mu_{N-1}$$\n\nThe middle portion depends on the value of $N$, and the parameters\n$\\theta_{N-1}$ and $\\omega_{N-1}$ contained in the $(N-1)$\\'th $SU(N)$\ntransformation. This is thus a convenient, systematic way to construct\nthe $N$-dimensional Haar measure for the unitary group. As a final note,\neven though unitaries can be parametrized in different ways, the\nunderlying Haar measure is *unique*. This is a consequence of it being\nan invariant measure, as will be shown later.\n\nHaar-random matrices from the $QR$ decomposition\n================================================\n\nNice-looking math aside, sometimes you just need to generate a large\nnumber of high-dimensional Haar-random matrices. It would be very\ncumbersome to sample and keep track of the distributions of so many\nparameters; furthermore, the measure above requires you to parametrize\nyour operations in a fixed way. There is a much quicker way to perform\nthe sampling by taking a (slightly modified) [QR\ndecomposition](https://en.wikipedia.org/wiki/QR_decomposition) of\ncomplex-valued matrices. This algorithm is detailed in, and consists of\nthe following steps:\n\n1. Generate an $N \\times N$ matrix $Z$ with complex numbers $a+bi$\n where both $a$ and $b$ are normally distributed with mean 0 and\n variance 1 (this is sampling from the distribution known as the\n *Ginibre ensemble*).\n2. Compute a QR decomposition $Z = QR$.\n3. Compute the diagonal matrix\n $\\Lambda = \\hbox{diag}(R_{ii}/|R_{ii}|)$.\n4. Compute $Q^\\prime = Q \\Lambda$, which will be Haar-random.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from numpy.linalg import qr\n\ndef qr_haar(N):\n \"\"\"Generate a Haar-random matrix using the QR decomposition.\"\"\"\n # Step 1\n A, B = np.random.normal(size=(N, N)), np.random.normal(size=(N, N))\n Z = A + 1j * B\n\n # Step 2\n Q, R = qr(Z)\n\n # Step 3\n Lambda = np.diag([R[i, i] / np.abs(R[i, i]) for i in range(N)])\n\n # Step 4\n return np.dot(Q, Lambda)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let\\'s check that this method actually generates Haar-random unitaries\nby trying it out for $N=2$ and plotting on the Bloch sphere.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"@qml.qnode(dev)\ndef qr_haar_random_unitary():\n qml.QubitUnitary(qr_haar(2), wires=0)\n return qml.state()\n\nqr_haar_samples = [qr_haar_random_unitary() for _ in range(num_samples)]\nqr_haar_bloch_vectors = np.array([convert_to_bloch_vector(s) for s in qr_haar_samples])\nplot_bloch_sphere(qr_haar_bloch_vectors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected, we find our qubit states are distributed uniformly over the\nsphere. This particular method is what\\'s implemented in packages like\n`scipy`\\'s\n[unitary\\_group](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.unitary_group.html)\nfunction.\n\nNow, it\\'s clear that this method works, but it is also important to\nunderstand *why* it works. Step 1 is fairly straightforward\\-\\--the base\nof our samples is a matrix full of complex values chosen from a typical\ndistribution. This isn\\'t enough by itself, since unitary matrices also\nhave constraints\\-\\--their rows and columns must be orthonormal. These\nconstraints are where step 2 comes in\\-\\--the outcome of a generic QR\ndecomposition consists of an *orthonormal* matrix $Q$, and and upper\ntriangular matrix $R$. Since our original matrix was complex-valued, we\nend up with a $Q$ that is in fact already unitary. But why not stop\nthere? Why do we then perform steps 3 and 4?\n\nSteps 3 and 4 are needed because, while the QR decomposition yields a\nunitary, it is not a unitary that is properly Haar-random. In, it is\nexplained that a uniform distribution over unitary matrices should also\nyield a uniform distribution over the *eigenvalues* of those matrices,\ni.e., every eigenvalue should be equally likely. Just using the QR\ndecomposition out of the box produces an *uneven* distribution of\neigenvalues of the unitaries! This discrepancy stems from the fact that\nthe QR decomposition is not unique. We can take any unitary diagonal\nmatrix $\\Lambda$, and re-express the decomposition as\n$QR = Q\\Lambda \\Lambda^\\dagger R = Q^\\prime R^\\prime$. Step 3 removes\nthis redundancy by fixing a $\\Lambda$ that depends on $R$, leading to a\nunique value of $Q^\\prime = Q \\Lambda$, and a uniform distribution of\neigenvalues.\n\n::: {.admonition}\nTry it!\n\nUse the `qr_haar` function above to generate random unitaries and\nconstruct a distribution of their eigenvalues. Then, comment out the\nlines for steps 3 and 4 and do the same\\-\\--you\\'ll find that the\ndistribution is no longer uniform. Check out reference for additional\ndetails and examples.\n:::\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fun (and not-so-fun) facts\n==========================\n\nWe\\'ve now learned what the Haar measure is, and both an analytical and\nnumerical means of sampling quantum states and unitary operations\nuniformly at random. The Haar measure also has many neat properties that\nplay a role in quantum computing.\n\nInvariance of the Haar measure\n------------------------------\n\nEarlier, we showed that the Haar measure is used when integrating\nfunctions over the unitary group:\n\n$$\\int_{V \\in U(N)} f(V) d\\mu_N(V).$$\n\nOne of the defining features of the Haar measure is that it is both left\nand right *invariant* under unitary transformations. That is,\n\n$$\\int_{V \\in U(N)} f(\\color{red}{W}V) d\\mu_N(V) = \\int_{V \\in U(N)} f(V\\color{red}{W}) d\\mu_N(V) = \\int_{V \\in U(N)} f(V) d\\mu_N(V).$$\n\nThis holds true for *any* other $N\\times N$ unitary $W$! A consequence\nof such invariance is that if $V$ is Haar-random, then so is $V^T,$\n$V^\\dagger,$ and any product of another unitary matrix and $V$ (where\nthe product may be taken on either side).\n\nAnother consequence of this invariance has to do with the structure of\nthe entries themselves: they must all come from the same distribution.\nThis is because the measure remains invariant under permutations, since\npermutations are unitary\\-\\--the whole thing still has to be Haar random\nno matter how the entries are ordered, so all distributions must be the\nsame. The specific distribution is complex numbers $a+bi$ where both $a$\nand $b$ has mean 0 and variance $1/N$[^1] (so, much like Ginibre\nensemble we used in the QR decomposition above, but with a different\nvariance and constraints due to orthonormality).\n\nConcentration of measure\n------------------------\n\nAn unfortunate (although interesting) property of the Haar measure is\nthat it suffers from the phenomenon of [concentration of\nmeasure](https://en.wikipedia.org/wiki/Concentration_of_measure). Most\nof the \\\"stuff\\\" in the space concentrates around a certain area, and\nthis gets worse as the size of the system increases. You can see the\nbeginnings of by looking at the sphere. For the 3-dimensional sphere, we\nsaw graphically how there is concentration around the equator, and how\nthe measure takes that into account with the additional factor of\n$\\sin \\theta$. This property becomes increasingly prominent for\n[higher-dimensional spheres](https://en.wikipedia.org/wiki/N-sphere).\n\n::: {.important}\n::: {.admonition-title}\nImportant\n:::\n\nThe concentration described here is not referring to what we witnessed\nearlier on, when we sampled quantum states (points on the Bloch sphere)\nincorrectly and found that they clustered around the poles. However,\nthat is not unrelated. Concentration of measure refers to where the\nmeasure itself is concentrated, and which parts of the space should be\nmore heavily weighted. For the case of the sphere, it is the equatorial\narea, and when we didn\\'t sample properly and take that concentration\ninto account, we obtained an uneven distribution.\n:::\n\nLet\\'s consider an $N$-dimensional unit sphere. Points on the sphere, or\nvectors in this space, are parametrized by $N-1$ real coordinates.\nSuppose we have some function $f$ that maps points on that sphere to\nreal numbers. Sample a point $x$ on that sphere from the uniform\nmeasure, and compute the value of $f(x)$. How close do you think the\nresult will be to the mean value of the function, $E[f]$, over the\nentire sphere?\n\nA result called [Levy\\'s\nlemma](https://en.wikipedia.org/wiki/Concentration_of_measure#Concentration_on_the_sphere)\n[^2],[^3] expresses how likely it is that $f(x)$ is a specific distance\naway from the mean. It states that, for an $x$ selected uniformly at\nrandom, the probability that $f(x)$ deviates from $E[f]$ by some amount\n$\\epsilon$ is bounded by:\n\n$$\\hbox{Pr}(|f(x) - E[f]| \\ge \\epsilon) \\leq 2 \\exp\\left[-\\frac{N\\epsilon^2}{9\\pi^3 \\eta^2}\\right].$$\n\nA constraint on the function $f$ is that it must be [Lipschitz\ncontinuous](https://en.wikipedia.org/wiki/Lipschitz_continuity), where\n$\\eta$ is the *Lipschitz constant* of the function. The important aspect\nhere is the likelihood of deviating significantly from the mean by an\namount $\\epsilon$ decreases exponentially with $\\epsilon.$ Furthermore,\nincreasing the dimension $N$ also makes the deviation exponentially less\nlikely.\n\nNow, this result seems unrelated to quantum states\\-\\--it concerns\nhigher-dimensional spheres. However, recall that a quantum state vector\nis a complex vector whose squared values sum to 1, similar to vectors on\na sphere. If you \\\"unroll\\\" a quantum state vector of dimension\n$N = 2^n$ by stacking its real and complex parts, you end with a vector\nof length $2 \\cdot 2^{n}$ which ends up behaving just like a unit vector\non the sphere in this dimension. Given that measure concentrates on\nspheres, and quantum state vectors can be converted to vectors on\nspheres, functions on random quantum states will also demonstrate\nconcentration.\n\nThis is bad news! To do useful things in quantum computing, we need a\nlot of qubits. But the more qubits we have, the more our randomly\nsampled states will look the same (specifically, random states will\nconcentrate around the maximally entangled state[^4]). This has\nimportant consequences for near-term algorithms (as detailed in the next\nsection), and any algorithm that involves uniform sampling of quantum\nstates and operations.\n\nHaar measure and barren plateaus\n--------------------------------\n\nSuppose you are venturing out to solve a new problem using an algorithm\nsuch as the\n`variational quantum eigensolver `{.interpreted-text\nrole=\"doc\"}. A critical component of such methods is the choice of\n`variational ansatz\n`{.interpreted-text role=\"doc\"}. Having now\nlearned a bit about the properties of the Haar measure, you may think it\nwould make sense to use this for the parametrization. Variational\nansaetze are, after all, parametrized quantum circuits, so why not\nchoose an ansatz that corresponds directly to a parametrization for\nHaar-random unitaries? The initial parameter selection will give you a\nstate in the Hilbert space uniformly at random. Then, since this ansatz\nspans the entire Hilbert space, you\\'re guaranteed to be able to\nrepresent the target ground state with your ansatz, and it should be\nable to find it with no issue \\... right?\n\nUnfortunately, while such an ansatz is extremely *expressive* (i.e., it\nis capable of representing any possible state), these ansaetze actually\nsuffer the most from the barren plateau problem[^5],[^6].\n`Barren plateaus `{.interpreted-text\nrole=\"doc\"} are regions in the cost landscape of a parametrized circuit\nwhere both the gradient and its variance approach 0, leading the\noptimizer to get stuck in a local minimum. This was explored recently in\nthe work of[^7], wherein closeness to the Haar measure was actually used\nas a metric for expressivity. The closer things are to the Haar measure,\nthe more expressive they are, but they are also more prone to exhibiting\nbarren plateaus.\n\n![Image source:[^8]. A highly expressive ansatz that can access much of\nthe space of possible unitaries (i.e., an ansatz capable of producing\nunitaries in something close to a Haar-random manner) is very likely to\nhave flat cost landscapes and suffer from the barren plateau\nproblem.](/demonstrations/haar_measure/holmes-costlandscapes.png){.align-center\nwidth=\"50.0%\"}\n\nIt turns out that the types of ansaetze know as *hardware-efficient\nansaetze* also suffer from this problem if they are \\\"random enough\\\"\n(this notion will be formalized in a future demo). It was shown in[^9]\nthat this is a consequence of the concentration of measure phenomenon\ndescribed above. The values of gradients and variances can be computed\nfor classes of circuits on average by integrating with respect to the\nHaar measure, and it is shown that these values decrease exponentially\nin the number of qubits, and thus huge swaths of the cost landscape are\nsimply and fundamentally flat.\n\nConclusion\n==========\n\nThe Haar measure plays an important role in quantum\ncomputing\\-\\--anywhere you might be dealing with sampling random\ncircuits, or averaging over all possible unitary operations, you\\'ll\nwant to do so with respect to the Haar measure.\n\nThere are two important aspects of this that we have yet to touch upon,\nhowever. The first is whether it is efficient to sample from the Haar\nmeasure\\-\\--given that the number of parameters to keep track of is\nexponential in the number of qubits, certainly not. But a more\ninteresting question is do we *need* to always sample from the full Haar\nmeasure? The answer to this is \\\"no\\\" in a very interesting way.\nDepending on the task at hand, you may be able to take a shortcut using\nsomething called a *unitary design*. In an upcoming demo, we will\nexplore the amazing world of unitary designs and their applications!\n\nReferences\n==========\n\n[^1]: E. Meckes (2019) [\\\"The Random Matrix Theory of the Classical\n Compact\n Groups\\\"](https://case.edu/artsci/math/esmeckes/Haar_book.pdf),\n Cambridge University Press.\n\n[^2]: M. Gerken (2013) \\\"Measure concentration: Levy\\'s Lemma\\\"\n ([lecture\n notes](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.679.2560)).\n\n[^3]: P. Hayden, D. W. Leung, and A. Winter (2006) \\\"Aspects of generic\n entanglement\\\", [Comm. Math. Phys. Vol. 265, No. 1, pp.\n 95-117](https://link.springer.com/article/10.1007%2Fs00220-006-1535-6).\n ([arXiv](https://arxiv.org/abs/quant-ph/0407049))\n\n[^4]: P. Hayden, D. W. Leung, and A. Winter (2006) \\\"Aspects of generic\n entanglement\\\", [Comm. Math. Phys. Vol. 265, No. 1, pp.\n 95-117](https://link.springer.com/article/10.1007%2Fs00220-006-1535-6).\n ([arXiv](https://arxiv.org/abs/quant-ph/0407049))\n\n[^5]: J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H.\n Neven (2018) \\\"Barren plateaus in quantum neural network training\n landscapes\\\", [Nature Communications,\n 9(1)](http://dx.doi.org/10.1038/s41467-018-07090-4).\n ([arXiv](https://arxiv.org/abs/1803.11173))\n\n[^6]: Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles (2021)\n \\\"Connecting ansatz expressibility to gradient magnitudes and barren\n plateaus\\\". ([arXiv](https://arxiv.org/abs/2101.02138))\n\n[^7]: Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles (2021)\n \\\"Connecting ansatz expressibility to gradient magnitudes and barren\n plateaus\\\". ([arXiv](https://arxiv.org/abs/2101.02138))\n\n[^8]: Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles (2021)\n \\\"Connecting ansatz expressibility to gradient magnitudes and barren\n plateaus\\\". ([arXiv](https://arxiv.org/abs/2101.02138))\n\n[^9]: J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H.\n Neven (2018) \\\"Barren plateaus in quantum neural network training\n landscapes\\\", [Nature Communications,\n 9(1)](http://dx.doi.org/10.1038/s41467-018-07090-4).\n ([arXiv](https://arxiv.org/abs/1803.11173))\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 0
}