scnc
/
qiskit-documentation
mirror of https://github.com/Qiskit/documentation.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
qiskit-documentation/learning/courses/quantum-diagonalization-alg.../krylov.ipynb

2688 lines
132 KiB
Plaintext
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"cells": [
{
"cell_type": "markdown",
"id": "2b1776f4-c259-494e-a33d-1ebb0459426c",
"metadata": {},
"source": [
"{/* cspell:ignore Hndt longrightarrow Bigg Jörg Liesen Zdenek Strakos Yousef Saad MINRES vstar nabla */}\n",
"\n",
"# Krylov quantum diagonalization\n",
"\n",
"In this lesson on Krylov quantum diagonalization (KQD) we will answer the following:\n",
"\n",
"* What is the Krylov method, generally?\n",
"* Why does the Krylov method work and under what conditions?\n",
"* How does quantum computing play a role?\n",
"\n",
"The quantum part of the calculations are based largely on work in Ref [\\[1\\]](#references).\n",
"\n",
"The video below gives an overview of Krylov methods in classical computing, motivates their use, and explains how quantum computing can play a role in that workstream. The subsequent text offers more detail and implements a Krylov method both classically, and using a quantum computer.\n",
"\n",
"<ImageLink title=\"Quantum Krylov diagonalization\" alt=\"QKD video thumbnail\" href=\"https://video.ibm.com/recorded/134360991\" src=\"/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-thumb.avif\"/>"
]
},
{
"cell_type": "markdown",
"id": "a9ba88fb-734e-44d8-ada4-c588e2654921",
"metadata": {},
"source": [
"## 1. Introduction to Krylov methods\n",
"\n",
"A __Krylov subspace method__ can refer to any of several methods built around what is called the __Krylov subspace__. A complete review of these is beyond the scope of this lesson, but Ref [\\[2-4\\]](#references) can all give substantially more background. Here, we will focus on what a Krylov subspace is, how and why it is useful in solving eigenvalue problems, and finally how it can be implemented on a quantum computer."
]
},
{
"cell_type": "markdown",
"id": "9b47a985-ec1d-41e2-8516-6cfcd7713e5f",
"metadata": {},
"source": [
"__Definition:__ Given a symmetric, positive semi-definite $N\\times N$ matrix $A$, the Krylov space $\\mathcal{K}^r$ of order $r$ is the space spanned by vectors obtained by multiplying higher powers of a matrix $A$, up to $r-1\\leq N$, with a reference vector $\\vert v \\rangle$.\n",
"\n",
"$$\n",
"\\mathcal{K}^r = \\text{span}\\left\\{ \\vert v \\rangle, A \\vert v \\rangle, A^2 \\vert v \\rangle, ..., A^{r-1} \\vert v \\rangle \\right\\}\n",
"$$\n",
"\n",
"Although the vectors above span what we call a Krylov subspace, there is no reason to think that they will be orthogonal. One often uses an iterative orthonormalizing process similar to __Gram-Schmidt orthogonalization__. Here the process is slightly different since each new vector is made orthogonal to the others as it is generated. In this context this is called __Arnoldi iteration__. Starting with the initial vector $|v\\rangle$, one generates the next vector $A|v\\rangle$, and then ensures that this second vector is orthogonal to the first by subtracting off its projection on $|v\\rangle$. That is\n",
"\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"|v_0\\rangle &=\\frac{|v\\rangle}{\\left|\\left| |v\\rangle \\right|\\right|}\\\\\n",
"\n",
"|v_1\\rangle &=\\frac{A|v\\rangle-\\langle v_0|A|v\\rangle |v_0\\rangle}{\\left|\\left|A|v\\rangle-\\langle v_0|A|v\\rangle |v_0\\rangle \\right|\\right|}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"It is now easy to see that $|v_0\\rangle \\perp |v_1\\rangle,$ since\n",
"$$\n",
"\\langle v_0 | v_1\\rangle=\\frac{\\langle v_0 | A|v\\rangle-\\langle v_0 |A|v\\rangle\\langle v_0|v_0\\rangle}{\\left|\\left| A|v\\rangle-\\langle A|v\\rangle|v_0\\rangle |v_0\\rangle \\right|\\right|}=0\n",
"$$\n",
"\n",
"We do the same for the next vector, ensuring it is orthogonal to both the previous two:\n",
"$$\n",
"|v_2\\rangle=\\frac{A |v_1\\rangle-\\langle v_0|A |v_1\\rangle |v_0\\rangle-\\langle v_1|A |v_1\\rangle |v_1\\rangle}{\\left|\\left| A |v_1\\rangle-\\langle v_0|A |v_1\\rangle |v_0\\rangle-\\langle v_1|A |v_1\\rangle |v_1\\rangle\\right|\\right|}\n",
"$$\n",
"If we repeat this process for all $r$ vectors, we have a complete orthonormal basis for a Krylov space. Note that the orthogonalization process here will yield zero once $r>m$, since $m$ orthogonal vector necessarily span the full space. The process will also yield zero if any vector is an eigenvector of $A$ since all subsequent vectors will be multiples of that vector."
]
},
{
"cell_type": "markdown",
"id": "5351576c-bc3c-4e21-8839-9c14eb10747f",
"metadata": {},
"source": [
"### 1.1 A simple example: Krylov by hand\n",
"\n",
"Let us step through a generation of a Krylov subspace generation on a trivially small matrix, so that we can see the process. We start with an initial matrix $A$ of interest to us:\n",
"$$\n",
"A=\\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\n",
"$$\n",
"For this small example, we can determine the eigenvectors and eigenvalues easily even by hand. We show the numerical solution here."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49017c87-a971-4d2d-b6dc-b8ec9b64184a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The eigenvalues are [2.58578644 4. 5.41421356]\n",
"The eigenvectors are [[ 5.00000000e-01 -7.07106781e-01 5.00000000e-01]\n",
" [ 7.07106781e-01 1.37464400e-16 -7.07106781e-01]\n",
" [ 5.00000000e-01 7.07106781e-01 5.00000000e-01]]\n"
]
}
],
"source": [
"# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use linalg.eig in this lesson.\n",
"\n",
"import numpy as np\n",
"\n",
"A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])\n",
"eigenvalues, eigenvectors = np.linalg.eig(A)\n",
"print(\"The eigenvalues are \", eigenvalues)\n",
"print(\"The eigenvectors are \", eigenvectors)"
]
},
{
"cell_type": "markdown",
"id": "a0f96bc7-fef9-48c7-858a-5b44772d710f",
"metadata": {},
"source": [
"We record them here for later comparison:\n",
"$$\n",
"\\begin{aligned}\n",
"a_0&=2.59,&|0\\rangle&=&\\begin{pmatrix}1/2\\\\-\\sqrt{2}/2\\\\1/2\\end{pmatrix}\\\\\n",
"\\\\\n",
"a_1&=4,&|1\\rangle&=&\\begin{pmatrix}\\sqrt{2}/2\\\\0\\\\-\\sqrt{2}/2\\end{pmatrix}\\\\\n",
"\\\\\n",
"a_2&=5.41,&|2\\rangle&=&\\begin{pmatrix}1/2\\\\\\sqrt{2}/2\\\\1/2\\end{pmatrix}\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "7f237ba1-d380-4553-b5bc-94422b4f39c1",
"metadata": {},
"source": [
"We would like to study how this process works (or fails) as we increase the dimension of our Krylov subspace, $r$. To this end, we will apply this process:\n",
"\n",
"* Generate a subspace of the full vector space starting with a randomly-chosen vector $|v\\rangle$ (call it $|v_0\\rangle$ if it is already normalized, as above).\n",
"* Project the full matrix $A$ onto that subspace, and find the eigenvalues of that projected matrix $\\tilde{A}$.\n",
"* Increase the size of the subspace by generating more vectors, ensuring that they are orthonormal, using a process similar to Gram-Schmidt orthogonalization.\n",
"* Project $A$ onto the larger subspace and find the eigenvalues of the resulting matrix, $\\tilde{A}$.\n",
"* Repeat this until the eigenvalues converge (or in this toy case, until you have generated vectors spanning the full vector space of the original matrix $A$).\n",
"\n",
"A normal implementation of the Krylov method would not need to solve the eigenvalue problem for the matrix projected on every Krylov subspace as it is built. You could construct the subspace of the desired dimension, project the matrix onto that subspace, and diagonalize the projected matrix. Projecting and diagonalizing at each subspace dimension is only done for checking convergence.\n",
"\n",
"#### Dimension $r=1$:\n",
"\n",
"We choose a random vector, say\n",
"$$\n",
"|v_0\\rangle=\\begin{pmatrix}1\\\\0\\\\0\\end{pmatrix}\n",
"$$\n",
"\n",
"If it is not already normalized, normalize it.\n",
"\n",
"We now project our matrix $A$ onto the subspace of this one vector:\n",
"$$\n",
"\\tilde{A}_0=\\langle v_0| A|v_0\\rangle=\\begin{pmatrix}1&0&0\\end{pmatrix}\\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}1\\\\0\\\\0\\end{pmatrix}=(4)\n",
"$$\n",
"This is our projection of the matrix onto our Krylov subspace when it contains just a single vector, $|v_0\\rangle$. The eigenvalue of this matrix is trivially 4. We can think of this as our zeroth-order estimate of the eigenvalues (in this case just one) of $A$. Although it is a poor estimate, it is the correct order of magnitude."
]
},
{
"cell_type": "markdown",
"id": "a0eeb9b1-86a5-4bc8-8651-66588aea7686",
"metadata": {},
"source": [
"#### Dimension $r=2$:\n",
"\n",
"We now generate the next vector in our subspace through operation with $A$ on the previous vector:\n",
"\n",
"$$\n",
"A|v_0\\rangle=\\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}1\\\\0\\\\0\\end{pmatrix}=\\begin{pmatrix}4\\\\-1\\\\0\\end{pmatrix}\n",
"$$\n",
"\n",
"Now we subtract off the projection of this vector onto our previous vector to ensure orthogonality.\n",
"$$\n",
"|v_1\\rangle=A|v_0\\rangle-\\langle v_0 |A|v_0\\rangle|v_0\\rangle\n",
"$$\n",
"$$\n",
"|v_1\\rangle=\\begin{pmatrix}4\\\\-1\\\\0\\end{pmatrix}-\\begin{pmatrix}1& 0& 0\\end{pmatrix}\\begin{pmatrix}4\\\\-1\\\\0\\end{pmatrix}\\begin{pmatrix}1\\\\0\\\\0\\end{pmatrix}=\\begin{pmatrix}0\\\\-1\\\\0\\end{pmatrix}\n",
"$$\n",
"\n",
"If it is not already normalized, normalize it. In this case, the vector was already normalized, so\n",
"\n",
"$$\n",
"|v_1\\rangle=\\begin{pmatrix}0\\\\-1\\\\0\\end{pmatrix}\n",
"$$\n",
"\n",
"We now project our matrix A onto the subspace of these two vectors:\n",
"$$\n",
"\\tilde{A}_1= \\begin{pmatrix} 1&0&0\\\\0&-1&0 \\end{pmatrix} \\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}1&0\\\\0&-1\\\\0&0\\end{pmatrix}=\\begin{pmatrix}1&0&0\\\\0&-1&0\\end{pmatrix}\\begin{pmatrix}4&1\\\\-1&-4\\\\0&1\\end{pmatrix}=\\begin{pmatrix}4&1\\\\1&4\\end{pmatrix}\n",
"$$\n",
"\n",
"We are still left with the problem of determining the eigenvalues of this matrix. But this matrix is slightly smaller than the full matrix. In problems involving very large matrices, working with this smaller subspace may be highly advantageous.\n",
"\n",
"$$\n",
"\\det(\\tilde{A_1}-\\lambda I)=0\n",
"$$\n",
"$$\n",
"\\begin{vmatrix} 4-\\lambda&1\\\\1&4-\\lambda\\end{vmatrix} =(4-\\lambda)^2-1=0\n",
"$$\n",
"$$\n",
"4-\\lambda=±1→\\lambda=3,5\n",
"$$\n",
"Although this is still not a good estimate, it is better than the zeroth order estimate. We will carry this out for one more iteration, to ensure the process is clear. However, this undercuts the point of the method, since we will end up diagonalizing a 3x3 matrix in the next iteration, meaning we have not saved time or computational power."
]
},
{
"cell_type": "markdown",
"id": "43ef6c6d-ba18-4e8f-94df-83c404d0015a",
"metadata": {},
"source": [
"#### Dimension $r=3$:\n",
"\n",
"We now generate the next vector in our subspace through operation with A on the previous vector:\n",
"$$\n",
"A|v_1\\rangle=\\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}0\\\\-1\\\\0\\end{pmatrix}=\\begin{pmatrix}1\\\\-4\\\\1\\end{pmatrix}\n",
"$$\n",
"Now we subtract off the projection of this vector onto both our previous vectors to ensure orthogonality.\n",
"$$\n",
"\\begin{aligned}\n",
"|v_2\\rangle&=A|v_1\\rangle-\\langle v_0 |A|v_1\\rangle|v_0\\rangle-\\langle v_1 |A|v_1\\rangle|v_1\\rangle\\\\\n",
"|v_2\\rangle&=\\begin{pmatrix}1\\\\-4\\\\1\\end{pmatrix}-\\begin{pmatrix}1& 0& 0 \\end{pmatrix}\\begin{pmatrix}1\\\\-4\\\\1\\end{pmatrix}\\begin{pmatrix}1\\\\0\\\\0\\end{pmatrix}-\\begin{pmatrix}0&-1& 0\\end{pmatrix}\\begin{pmatrix}1\\\\-4\\\\1\\end{pmatrix}\\begin{pmatrix}0\\\\-1\\\\0\\end{pmatrix}=\\begin{pmatrix}0\\\\0\\\\1\\end{pmatrix}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"If it is not already normalized, normalize it. In this case, the vector was already normalized, so\n",
"$$\n",
"|v_2 \\rangle=\\begin{pmatrix}0\\\\0\\\\1\\end{pmatrix}\n",
"$$\n",
"We now project our matrix $A$ onto the subspace of these vectors:\n",
"\n",
"$$\n",
"\\tilde{A}_2=\\begin{pmatrix}1&0&0\\\\0&-1&0\\\\0&0&1\\end{pmatrix}\\begin{pmatrix}4&-1&0\\\\-1&4&-1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}1&0&0\\\\0&-1&0\\\\0&0&1\\end{pmatrix}=\\begin{pmatrix}4&-1&0\\\\1&-4&1\\\\0&-1&4\\end{pmatrix}\\begin{pmatrix}1&0&0\\\\0&-1&0\\\\0&0&1\\end{pmatrix}=\\begin{pmatrix}4&1&0\\\\1&4&1\\\\0&1&4\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "5f19b90b-93fa-4ccb-b8a2-3bbd55e4dab0",
"metadata": {},
"source": [
"We now determine the eigenvalues:\n",
"$$\n",
"\\det(\\tilde{A}_2-\\lambda I)=0\n",
"$$\n",
"$$\n",
"\\begin{vmatrix}4-\\lambda&1&0\\\\1&4-\\lambda&1\\\\0&1&4-\\lambda\\end{vmatrix} = (4-\\lambda)((4-\\lambda)^3-1)-(4-\\lambda)=0\\\\\n",
"$$\n",
"$$\n",
"4-\\lambda=0,4-\\lambda=±2^{1/2}→\\lambda=4-2^{1/2},4,4+2^{1/2}≈2.59,4,5.41\n",
"$$\n",
"These eigenvalues are exactly the eigenvalues of the original matrix $A$. This must be the case, since we have expanded our Krylov subspace to span the entire vector space of the original matrix $A$.\n",
"\n",
"In this example, the Krylov method may not appear particularly easier than direct diagonalization. Indeed, as we will see in later sections, the Krylov method is only advantageous above a certain matrix dimension; this is intended to help us solve eigenvalue/eigenvector problems of extremely large matrices.\n",
"\n",
"![An image showing a very large matrix being projected onto a Krylov subspace, that is, rows of Krylov vectors making a matrix on the left, a Hamiltonian, then columns of Krylov vectors on the right.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig1.avif)\n",
"\n",
"This is the only example we will show worked “by hand”, but section 2 below shows computational examples.\n",
"\n",
"#### Clarification of terms\n",
"\n",
"A common misconception is that there is just a single Krylov subspace for a given problem. But of course, since there are many initial vectors to which our matrix could be applied, there are many possible Krylov subspaces. We will only use the phrase \"__the__ Krylov subspace\" to refer to a specific Krylov subspace already defined for a specific example. For general problem-solving approaches we will refer to \"__a__ Krylov subspace\".\n",
"A final clarification is that it is valid to refer to a \"Krylov __space__\". One often sees it called a \"Krylov __subspace__\" because of its use in the context of projecting matrices from an initial space into a subspace. In keeping with that context, we will mostly refer to it as a subspace here."
]
},
{
"cell_type": "markdown",
"id": "f3117328-61e0-45a2-b8f2-05356d210fae",
"metadata": {},
"source": [
"#### Check-in question\n",
"\n",
"Explain why it is not (a) useful, and (b) possible to extend the dimension of the Krylov subspace $r$ beyond the dimension $N$ of the matrix of interest.\n",
"\n",
"__Answer:__\n",
"\n",
"(a) Since we are orthonormalizing the vectors as we produce them, a set of $N$ such vectors will form a complete basis, meaning a linear combination of them can be used to create any vector in the space. (b) The orthogonalization process consists of subtracting off the projection of a new vector onto all previous vectors. If all previous vectors span the full vector space, then subtracting off projections onto the full subspace will always leave us with a zero vector.\n",
"\n",
"#### Check-in question\n",
"\n",
"Suppose a fellow researcher is demonstrating the Krylov method applied to a small toy matrix for some colleagues. This fellow researcher plans to use the matrix and initial vector:\n",
"$$\n",
"A=\\begin{pmatrix}2&1&3\\\\1&2&3\\\\3&3&5\\end{pmatrix}\n",
"$$\n",
"and\n",
"$$\n",
"|\\psi\\rangle=\\frac{1}{\\sqrt{2}}\\begin{pmatrix}1\\\\-1\\\\0\\end{pmatrix}.\n",
"$$\n",
"Is there something wrong with this? How would you advise your colleague?\n",
"\n",
"__Answer:__\n",
"\n",
"Your colleague has accidentally chosen an eigenvector for his/her initial vector. Acting with the matrix on the initial vector will simply return the same vector back, scaled by the eigenvalue. This will not generate a subspace of increasing dimension. Advise your colleague to select a different initial vector, making sure it is not an eigenvector.\n",
"\n",
"\n",
"#### Check-in question\n",
"\n",
"Apply the Krylov method to the matrix below, selecting an appropriate new initial vector. Write down the estimates of the minimum eigenvalue at the 0th and 1st order of your Krylov subspace.\n",
"$$\n",
"A=\\begin{pmatrix}1&1&0\\\\1&1&1\\\\0&1&1\\end{pmatrix}\n",
"$$\n",
"\n",
"__Answer:__\n",
"\n",
"There are many possible answers depending on the choice of initial vector. We will choose:\n",
"$$\n",
"|v_0\\rangle=\\frac{1}{\\sqrt{3}}\\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix}.\n",
"$$\n",
"To get $|v_1\\rangle$ we apply $A$ once to $|v_0\\rangle$, and then make $|v_1\\rangle$ orthogonal to $|v_0\\rangle.$\n",
"$$\n",
"A|v_0\\rangle=\\begin{pmatrix}1&1&0\\\\1&1&1\\\\0&1&1\\end{pmatrix}\\frac{1}{\\sqrt{3}}\\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix} = \\frac{1}{\\sqrt{3}}\\begin{pmatrix}2\\\\3\\\\2\\end{pmatrix}\n",
"$$\n",
"$$\n",
"A|v_0\\rangle - \\langle v_0|A|v_0\\rangle |v_0\\rangle=\\frac{1}{\\sqrt{3}}\\begin{pmatrix}2\\\\3\\\\2\\end{pmatrix} - \\frac{1}{\\sqrt{3}}\\begin{pmatrix}1&1&1\\end{pmatrix}\\frac{1}{\\sqrt{3}}\\begin{pmatrix}2\\\\3\\\\2\\end{pmatrix}\\frac{1}{\\sqrt{3}}\\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix} = \\frac{1}{\\sqrt{3}}\\begin{pmatrix}2\\\\3\\\\2\\end{pmatrix}-\\frac{7}{3}\\frac{1}{\\sqrt{3}}\\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix}=\\sqrt{\\frac{3}{2}}\\begin{pmatrix}-1/3\\\\2/3\\\\-1/3\\end{pmatrix}\n",
"$$\n",
"At 0th order, the projection onto our Krylov subspace is\n",
"\n",
"$$\n",
"\\langle v_0|A|v_0\\rangle=\\frac{1}{\\sqrt{3}}\\begin{pmatrix}1&1&1\\end{pmatrix} \\begin{pmatrix}1&1&0\\\\1&1&1\\\\0&1&1\\end{pmatrix} \\frac{1}{\\sqrt{3}}\\begin{pmatrix}1\\\\1\\\\1\\end{pmatrix} = \\frac{7}{3}\n",
"$$\n",
"At 1st order, the projection onto this Krylov subspace is\n",
"\n",
"$$\n",
"\\langle V^1|A|V^1\\rangle=\\begin{pmatrix}\\frac{1}{\\sqrt{3}}&\\frac{1}{\\sqrt{3}}&\\frac{1}{\\sqrt{3}}\\\\-\\sqrt{\\frac{1}{6}}&\\sqrt{\\frac{2}{3}}&-\\sqrt{\\frac{1}{6}}\\end{pmatrix} \\begin{pmatrix}1&1&0\\\\1&1&1\\\\0&1&1\\end{pmatrix} \\begin{pmatrix}\\frac{1}{\\sqrt{3}}&-\\sqrt{\\frac{1}{6}}\\\\\\frac{1}{\\sqrt{3}}& \\sqrt{\\frac{2}{3}} \\\\ \\frac{1}{\\sqrt{3}}&-\\sqrt{\\frac{1}{6}}\\end{pmatrix}\n",
"$$\n",
"This can be done by hand, but is most easily done using numpy:\n",
"\n",
"```\n",
"import numpy as np\n",
"vstar = np.array([[1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)],[-1/np.sqrt(6),np.sqrt(2/3),-1/np.sqrt(6)]]\n",
")\n",
"A = np.array([[1, 1, 0],\n",
" [1, 1, 1],\n",
" [0, 1, 1]])\n",
"v = np.array([[1/np.sqrt(3),-1/np.sqrt(6)],[1/np.sqrt(3),np.sqrt(2/3)],[1/np.sqrt(3),-1/np.sqrt(6)]])\n",
"proj = vstar@A@v\n",
"print(proj)\n",
"eigenvalues, eigenvectors = np.linalg.eig(proj)\n",
"print(\"The eigenvalues are \", eigenvalues)\n",
"print(\"The eigenvectors are \", eigenvectors)\n",
"```\n",
"outputs:\n",
"```\n",
"[[ 2.33333333 0.47140452]\n",
" [ 0.47140452 -0.33333333]]\n",
"The eigenvalues are [ 2.41421356 -0.41421356]\n",
"The eigenvectors are [[ 0.98559856 -0.16910198]\n",
" [ 0.16910198 0.98559856]]\n",
"```\n",
"The minimum eigenvalue estimate is -0.414."
]
},
{
"cell_type": "markdown",
"id": "5a4f969e-d346-485a-9e4c-a04dbc9492e9",
"metadata": {},
"source": [
"### 1.2 Types of Krylov methods\n",
"\n",
"\"Krylov subspace methods\" can refer to any of several iterative techniques used to solve large linear systems and eigenvalue problems. What they all have in common is that they construct an approximate solution from a Krylov subspace\n",
"\n",
"$$\\mathcal{K}^r(A,|v\\rangle ) = \\text{span}\\{|v\\rangle, A|v\\rangle, A^2|v\\rangle, ..., A^{r-1}|v\\rangle\\},$$\n",
"\n",
"where $|v\\rangle$ is the initial guess (see Ref [\\[5\\]](#references)). They differ in how they choose the best approximation from this subspace, balancing factors such as convergence rate, memory usage, and overall computational cost. The focus of this lesson is to leverage quantum computing in the context of Krylov subspace methods; an exhaustive discussion of these methods is beyond its scope. The brief definitions below are for context only and include some references for investigating these methods further.\n",
"\n",
"__The conjugate gradient (CG) method__: This method is used for solving symmetric, positive definite linear systems[\\[6\\]](#references). It minimizes the A-norm of the error at each iteration, making it particularly effective for systems arising from discretized elliptic PDEs[\\[7\\]](#references). We will use this approach in the next section to motivate why a Krylov subspace would be an effective subspace in which to probe for improved solutions to linear systems.\n",
"\n",
"__The generalized minimal residual (GMRES) method__: This is designed for solving general nonsymmetric linear systems. It minimizes the residual norm over a Krylov space at each iteration, making it robust but potentially memory-intensive for large systems[\\[7\\]](#references).\n",
"\n",
"__The minimal residual (MINRES) method__: This method is used for solving symmetric indefinite linear systems. It's similar to GMRES but takes advantage of the matrix symmetry to reduce computational cost[\\[8\\]](#references).\n",
"\n",
"Other approaches of note include the __full orthogonalization method (FOM)__, which is closely related to Arnoldi's method for eigenvalue problems, the __bi-conjugate gradient (BiCG) method__, and the __induced dimension reduction (IDR) method__."
]
},
{
"cell_type": "markdown",
"id": "774fa41b-9fbb-46cc-9c64-37f7a008e161",
"metadata": {},
"source": [
"### 1.3 Why the Krylov subspace method works\n",
"\n",
"Here we will motivate that the Krylov subspace method should be an efficient way to approximate matrix eigenvalues via iterative refinement of eigenvector approximations, through the lens of steepest descent. We will argue that given an initial guess of a ground state, the space of successive corrections to that initial guess that yields the fastest convergence is a Krylov subspace. We stop short of a rigorous proof of convergence behavior.\n",
"\n",
"Assume our matrix of interest $A$ is symmetric and positive definite. This makes our argument most relevant to the CG method above. We make no assumptions about sparsity here; nor are we claiming that $A$ must be a Hermitian (which it needs to be if it is a Hamiltonian).\n",
"\n",
"We typically wish to solve a problem of the form\n",
"$$\n",
"A|x\\rangle=|b\\rangle.\n",
"$$\n",
"One might imagine that $|b\\rangle=c|x\\rangle$ where $c$ is some constant, as in an eigenvalue problem. But our problem statement remains more general for now.\n",
"\n",
"We start with a vector $|x_0\\rangle$ that is an approximate solution. Although there are parallels between this guess $|x_0\\rangle$ and $|v_0\\rangle$ in Section 1.1, we are not leveraging these here. Our guess $|x_0\\rangle$ has error, which we call $|e_0\\rangle:$\n",
"$$\n",
"|e_0\\rangle:=|x\\rangle−|x_0\\rangle.\n",
"$$\n",
"We also define the residual $R_0:$\n",
"$$\n",
"|R_0\\rangle=|b\\rangle−A|x_0\\rangle.\n",
"$$\n",
"Here we use capital $R$ to distinguish the residual from the dimension of our Krylov subspace $r$.\n",
"\n",
"![A true eigenvector labeled x, a guess labeled x 0 and a graphical representation of hte error between those two.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig2.avif)\n",
"\n",
"We now want to make a correction step of the form\n",
"$$\n",
"|x_1\\rangle=|x_0\\rangle+|p_0\\rangle,\n",
"$$\n",
"which we hope improves our approximation. Here $|p_0\\rangle$ is some vector yet to be determined. Let $|e_1\\rangle$ be the error after the correction is made. Then\n",
"$$\n",
"|e_1\\rangle=|x\\rangle−|x_1\\rangle=|x\\rangle−(|x_0\\rangle+|p_0\\rangle)=|e_0\\rangle−|p_0\\rangle.\n",
"$$\n",
"\n",
"![A true eigenvector and an update to the initial guess. The updated guess is closer to the true eigenvector.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig3.avif)\n",
"\n",
"We are interested in how our error behaves when transformed by our matrix. So let us calculate the $A$-norm of the error. That is\n",
"$$\n",
"\\begin{aligned}\n",
"∥|e_0\\rangle−|p_0\\rangle∥_A^2&=\\left(\\langle e_0|A−\\langle p_0|A\\right)\\left(|e_0\\rangle−|p_0\\rangle\\right)\\\\\n",
" & = \\langle e_0|A|e_0 \\rangle − \\langle e_0|A|p_0\\rangle − \\langle p_0|A|e_0\\rangle+\\langle p_0|A|p_0\\rangle\\\\\n",
" & = \\langle e_0|A|e_0\\rangle−2\\langle e_0|A|p_0\\rangle+\\langle p_0|A|p_0\\rangle\\\\\n",
" & = d−2\\langle R_0|p_0\\rangle +\\langle p_0|A|p_0\\rangle,\n",
" \\end{aligned}\n",
"$$\n",
"where we have used the symmetry of $A$ and also that $A |e_0\\rangle = |R_0\\rangle.$ Here $d$ is some constant independent of $|p_0\\rangle$. As mentioned in Section 1.2, the $A$-norm of the error is not the only quantity we might choose to minimize, but it is a good one. We want to see how this quantity varies with our choice of correction vectors $|p_0\\rangle.$ So we define the function $f$ by setting\n",
"$$\n",
"f(|p_0\\rangle)=\\langle p_0|A|p_0\\rangle−2\\langle R_0|p_0\\rangle+d.\n",
"$$\n",
"$f$ is just the error $|e_1\\rangle$ as a function of the correction $|p_0\\rangle$ measured in the $A$-norm. Hence, we want to choose $|p_0\\rangle$ such that $f(|p_0\\rangle)$ is as small as possible. For this purpose, we compute the gradient of $f$. Using the symmetry of $A$ we have\n",
"$$\n",
"\\nabla f(|p_0\\rangle) = 2(A|p_0\\rangle−|R_0\\rangle).\n",
"$$\n",
"The gradient points in the direction of steepest ascent, meaning its opposite gives us the direction in which the function decreases the most: the direction of __steepest descent__. At our initial guess $|x_0\\rangle$, where $|p_0\\rangle=0$, we have that\n",
"$\\nabla f(0) = -2|R_0\\rangle.$\n",
"Thus, the function $f$ decreases the most in the direction of the residual $|R_0\\rangle.$ So our initial choice would benefit most by the addition of the vector $|p_0\\rangle=\\alpha_0 |R_0\\rangle$ for some scalar $\\alpha_0$.\n",
"\n",
"In the next step, we choose, again, a vector $|p_1\\rangle$ and add its value to the current approximation. Using the same argument as before we choose $|p_1\\rangle = \\alpha_1 |R_1\\rangle$ for some scalar $\\alpha_1$. We continue in this manner, such that the $k^\\text{th}$ iteration of our vector is\n",
"$$\n",
"|x_{k+1}\\rangle=|x_0\\rangle+\\alpha_0 |R_0\\rangle+\\alpha_1 |R_1\\rangle+⋯+\\alpha_k |R_k\\rangle.\n",
"$$\n",
"Equivalently, we want to build up the space from which we choose our improved estimates by adding $|R_0\\rangle$, $|R_1\\rangle$, and so on, in order. The $k^\\text{th}$ estimated vector lies in\n",
"$$\n",
"|x_{k+1}\\rangle\\in |x_0\\rangle+\\text{span}\\{|R_0\\rangle,|R_1\\rangle,…,|R_k\\rangle \\}.\n",
"$$\n",
"Now, using the relation that\n",
"$$\n",
"|R_{k+1}\\rangle=|b\\rangle−A |x_{k+1}\\rangle=|b\\rangle−A(|x_k\\rangle+\\alpha_k |R_k\\rangle)=|R_k\\rangle−\\alpha_k A |R_k\\rangle,\n",
"$$\n",
"we see that\n",
"$$\n",
"\\text{span} \\{|R_0\\rangle,|R_1\\rangle,…,|R_k\\rangle \\}=\\text{span} \\{|R_0\\rangle,A|R_0\\rangle,…,A^{k−1}|R_0\\rangle \\}.\n",
"$$\n",
"That is, the space we build up that most efficiently approximates the correct solution $|x\\rangle$ is exactly the space built up by successive operation of the matrix $A$ on $|R_0\\rangle.$ A Krylov subspace _is_ the space spanned by the vectors of successive directions of steepest descent.\n",
"\n",
"Finally we reiterate that we have made no numerical claims about the scaling of this approach, nor have we discussed the comparative benefit for sparse matrices. This is only meant to motivate the use of Krylov subspace methods, and add some intuitive sense for them. We will now explore the behavior of these methods numerically."
]
},
{
"cell_type": "markdown",
"id": "2df5e11f-62fe-4237-bcaf-36dcc88f11a7",
"metadata": {},
"source": [
"#### Check-in question\n",
"\n",
"In the workflow above, we proposed minimizing the $A$-norm of the error. What other quantities might one consider minimizing in seeking the ground state and its eigenvalue?\n",
"\n",
"__Answer:__\n",
"\n",
"One could imagine using the residual vector instead of the $A$-norm of the error. There might be cases in which considering the error vector itself is useful."
]
},
{
"cell_type": "markdown",
"id": "56bcf069-b3b5-4c51-b91f-7e9bb62c5c94",
"metadata": {},
"source": [
"## 2. Krylov methods in classical computation\n",
"\n",
"In this section we implement Arnoldi iterations computationally so that we may leverage a Krylov subspace in solving eigenvalue problems. We will apply this first to a small-scale example, then examine how computation time scales as the size of the matrix of interest increases. A key idea here will be that the generation of the vectors spanning the Krylov space will be a large contributor to total computing time required. The memory required will vary between specific Krylov methods. But memory constraints can limit the use of traditional Krylov methods."
]
},
{
"cell_type": "markdown",
"id": "952459b6-75dd-4044-8cec-7004fbc48ddf",
"metadata": {},
"source": [
"### 2.1 Simple small-scale example\n",
"\n",
"In the process of creating a Krylov subspace, we will need to orthonormalize the vectors in our subspace. Let us define a function that takes an established vector from our subspace `vknown` (not assumed to be normalized) and a candidate vector to add to our subspace `vnext` and make `vnext` orthogonal to `vknown` and normalized. Let us further define a function that steps through this process for all established vectors in our Krylov subspace to ensure a fully orthonormal set."
]
},
{
"cell_type": "code",
"execution_count": 100,
"id": "fb4f96ea-906c-4819-a9cb-8e9409cac051",
"metadata": {},
"outputs": [],
"source": [
"# vknown is some established vector in our subspace. vnext is one we wish to add, which must be orthogonal to vknown.\n",
"\n",
"\n",
"def orthog_pair(vknown, vnext):\n",
" vknown = vknown / np.sqrt(vknown.T @ vknown)\n",
" diffvec = vknown.T @ vnext * vknown\n",
" vnext = vnext - diffvec\n",
" return vnext\n",
"\n",
"\n",
"# v is the candidate vector to be added to our subspace. s is the existing subspace.\n",
"\n",
"\n",
"def orthoset(v, s):\n",
" v = v / np.sqrt(v.T @ v)\n",
" temp = v\n",
" for i in range(len(s)):\n",
" temp = orthog_pair(s[i], temp)\n",
" v = temp / np.sqrt(temp.T @ temp)\n",
" return v"
]
},
{
"cell_type": "markdown",
"id": "61a2704e-3985-4bd4-a9a9-253cc96c708b",
"metadata": {},
"source": [
"Now let us define a function that builds an iteratively larger and larger Krylov subspace, until the space of Krylov vectors spans the full space of the original matrix. This will enable us to see how well the eigenvalues obtained using our Krylov subspace method match the exact values, as a function of Krylov subspace dimension. Importantly, our function `krylov_full_build` returns the Krylov vectors, the projected Hamiltonians, the eigenvalues, and the time required."
]
},
{
"cell_type": "code",
"execution_count": 101,
"id": "818cceee-bd88-472a-8af9-0a6c0dec5646",
"metadata": {},
"outputs": [],
"source": [
"# Necessary imports and definitions to track time in microseconds\n",
"import time\n",
"\n",
"\n",
"def time_mus():\n",
" return int(time.time() * 1000000)\n",
"\n",
"\n",
"# This function\n",
"\n",
"\n",
"def krylov_full_build(v0, matrix):\n",
" t0 = time_mus()\n",
" b = v0 / np.sqrt(v0 @ v0.T)\n",
" A = matrix\n",
" ks = []\n",
" ks.append(b)\n",
" Hs = []\n",
" eigs = []\n",
" Hs.append(b.T @ A @ b)\n",
" eigs.append(np.array([b.T @ A @ b]))\n",
" k_tot_times = []\n",
"\n",
" for j in range(len(A) - 1):\n",
" vec = A @ ks[j].T\n",
" ortho = orthoset(vec, ks)\n",
" ks.append(ortho)\n",
" ksarray = np.array(ks)\n",
" Hs.append(ksarray @ A @ ksarray.T)\n",
" eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)\n",
" k_tot_times.append(time_mus() - t0)\n",
"\n",
" # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues, and the total time required.\n",
" return (ks, Hs, eigs, k_tot_times)"
]
},
{
"cell_type": "markdown",
"id": "4cd7bdcc-65fd-4199-b9fd-d0a9f62cba17",
"metadata": {},
"source": [
"We will test this on a matrix that is still quite small, but larger than we might want to do by hand."
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "3e2aa939-b568-4982-b64a-9f2e4b0da2b5",
"metadata": {},
"outputs": [],
"source": [
"# Define our small test matrix\n",
"test_matrix = np.array(\n",
" [\n",
" [4, -1, 0, 1, 0],\n",
" [-1, 4, -1, 2, 1],\n",
" [0, -1, 4, 3, 3],\n",
" [1, 2, 3, 4, 0],\n",
" [0, 1, 3, 0, 4],\n",
" ]\n",
")\n",
"\n",
"# Give the test matrix and an initial guess as arguments in the function defined above. Calculate outputs.\n",
"test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(\n",
" np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix\n",
")"
]
},
{
"cell_type": "markdown",
"id": "2f74fcf0-8959-4862-b33d-c8b1eddda2a8",
"metadata": {},
"source": [
"We can check our functions by ensuring that in the last step (when the Krylov space is the full vector space of the original matrix) the eigenvalues from the Krylov method exactly match those of the exact numerical diagonalization:"
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "704db2f6-5081-404d-9ae9-7264a62ea432",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-1.36956923 8.43756009 2.9040308 5.34436028 4.68361806]\n",
"[-1.36956923 8.43756009 2.9040308 4.68361806 5.34436028]\n"
]
}
],
"source": [
"print(np.linalg.eig(test_matrix).eigenvalues)\n",
"print(test_eigs[len(test_matrix) - 1])"
]
},
{
"cell_type": "markdown",
"id": "f0ddf041-d7e8-4131-bb4c-3310d8351c2c",
"metadata": {},
"source": [
"That was successful. Of course, what really matters is how good our approximation is as a function of the dimension of our Krylov subspace dimension. Because we are often concerned with finding ground states and other minimum eigenvalues (and for other more algebraic reasons explained below), let's look at our estimate of the lowest eigenvalue as a function of Krylov subspace dimension. That is"
]
},
{
"cell_type": "code",
"execution_count": 104,
"id": "48b23ca9-a5ba-4a9c-bbc2-7fba75b36124",
"metadata": {},
"outputs": [],
"source": [
"def errors(matrix, krylov_eigs):\n",
" targ_min = min(np.linalg.eig(matrix).eigenvalues)\n",
" err = []\n",
" for i in range(len(matrix)):\n",
" err.append(min(krylov_eigs[i]) - targ_min)\n",
" return err"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0912e12b-63bd-4026-ac8e-56bc5db8736e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/0912e12b-63bd-4026-ac8e-56bc5db8736e-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"krylov_error = errors(test_matrix, test_eigs)\n",
"\n",
"plt.plot(krylov_error)\n",
"plt.axhline(y=0, color=\"red\", linestyle=\"--\") # Add dashed red line at y=0\n",
"plt.xlabel(\"Order of Krylov subspace\") # Add x-axis label\n",
"plt.ylabel(\"Error in minimum eigenvalue\") # Add y-axis label\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7af1b947-205b-46f8-bb23-d439ae6f076a",
"metadata": {},
"source": [
"We see that the minimum eigenvalue is reached fairly accurately once the Krylov subspace has grown to $\\mathcal{K}^2,$ and is perfect by $\\mathcal{K}^3.$"
]
},
{
"cell_type": "markdown",
"id": "7ee80bf2-e3e8-4124-9afc-ef2ef247f831",
"metadata": {},
"source": [
"### 2.2 Time scaling with matrix dimension\n",
"\n",
"Let us convince ourselves that the Krylov method can be advantageous over exact numerical eigensolvers in the following way:\n",
"* Construct random matrices (not sparse, not the ideal application for KQD)\n",
"* Determine eigenvalues using two methods: directly using NumPy and using a Krylov subspace.\n",
"* We choose a cutoff for how precise our eigenvalues must be, before we accept the Krylov estimates.\n",
"* Compare the wall time required to solve in these two ways.\n",
"\n",
"__Caveats:__ As we will discuss in detail below, Krylov quantum diagonalization is best applied to operators whose matrix representations are sparse and/or can be written using a small number of groups of commuting Pauli operators. The random matrices we are using here do not fit that description. These are only useful in probing the scale at which classical Krylov methods might be useful.\n",
"Secondly, in using the Krylov method we will calculate eigenvalues using many different-sized Krylov subspaces. We will report the time required for the minimum-dimension Krylov subspace that achieves our required accuracy for the ground state eigenvalue. Again, this is a bit different from solving a problem that is intractable for exact eigensolvers, since we are using the exact solution to assess the dimension needed.\n",
"\n",
"We begin by generating our set of random matrices."
]
},
{
"cell_type": "code",
"execution_count": 106,
"id": "a5f8f1f6-b334-43e4-ad40-fd689c07a3e1",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"# Set the random seed\n",
"np.random.seed(42)\n",
"\n",
"# how many random matrices will we make\n",
"num_matrix = 200\n",
"\n",
"matrices = []\n",
"for m in range(1, num_matrix):\n",
" matrices.append(np.random.rand(m, m))"
]
},
{
"cell_type": "markdown",
"id": "d1c31d97-7e76-4577-ab03-ce8f0aebf134",
"metadata": {},
"source": [
"We now diagonalize each matrix directly, using numpy. We calculate the time required for diagonalization for later comparison."
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "d9cb6e36-9ad3-403f-823f-0072e9cf53fe",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/d9cb6e36-9ad3-403f-823f-0072e9cf53fe-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"matrix_numpy_times = []\n",
"matrix_numpy_eigs = []\n",
"for mm in range(num_matrix - 1):\n",
" t0 = time_mus()\n",
" matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))\n",
" matrix_numpy_times.append(time_mus() - t0)\n",
"\n",
"plt.plot(matrix_numpy_times)\n",
"plt.xlabel(\"Dimension of matrix\") # Add x-axis label\n",
"plt.ylabel(\"Time to diagonalize (microsec)\") # Add y-axis label\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "4e5d4af5-78da-4254-a726-4b0a72f98fe6",
"metadata": {},
"source": [
"Note that in the image above, the anomalously high time around a dimension of 125 may be due to the random nature of the matrices or due to implementation on the classical processor used, but it is not reproducible. Re-running the code will yield a different profile with different anomalous peaks.\n",
"\n",
"Now for each matrix we will build up a Krylov subspace and calculate eigenvalues in steps. At each step, we will check to see if the lowest eigenvalue has been obtained to within our specified absolute error. The subspace that first gives us eigenvalues within our specified error is the subspace for which we will record computation times. Executing this cell may take several minutes, depending on processor speed. Feel free to skip evaluation or reduce the maximum dimension of matrices diagonalized. Looking at the pre-calculated results is sufficient."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7fbe57ec-507a-412e-b7fe-09aaea28e103",
"metadata": {},
"outputs": [],
"source": [
"# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size at which that errer is achieved.\n",
"abserr = 0.05\n",
"accept_subspace_size = []\n",
"\n",
"# Lists to store total time spent on the Krylov method, and he the subset of that time spent on diagonalizing the projected matrix.\n",
"matrix_krylov_tot_times = []\n",
"matrix_krylov_dim = []\n",
"\n",
"# Step through all our random matrices\n",
"for mm in range(0, num_matrix - 1):\n",
" test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(\n",
" np.ones(len(matrices[mm])), matrices[mm]\n",
" )\n",
" # We have not yet found a Krylov subspace that produces our minimum eigenvalue to within the required error.\n",
" found = 0\n",
" for j in range(0, len(matrices[mm]) - 1):\n",
" # If we still haven't found the desired subspace...\n",
" if found == 0:\n",
" # ...but if this one satisfies the requirement, then record everything\n",
" if (\n",
" abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])\n",
" < abserr\n",
" ):\n",
" accept_subspace_size.append(j)\n",
" matrix_krylov_tot_times.append(test_k_tot_times[j])\n",
" matrix_krylov_dim.append(mm)\n",
" found = 1"
]
},
{
"cell_type": "markdown",
"id": "5046bfdb-9ff6-4142-9dc6-260ce26fd3f5",
"metadata": {},
"source": [
"Let us plot the times we have obtained for these two methods for comparison:"
]
},
{
"cell_type": "code",
"execution_count": 94,
"id": "ed3302d5-a0a8-479b-9731-755fc5a1b684",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/ed3302d5-a0a8-479b-9731-755fc5a1b684-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(matrix_numpy_times, color=\"blue\")\n",
"plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color=\"green\")\n",
"plt.xlabel(\"Dimension of matrix\") # Add x-axis label\n",
"plt.ylabel(\"Time to diagonalize (microsec)\") # Add y-axis label\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e3c5ff65-964f-4b7c-8ca7-89958b3d9465",
"metadata": {},
"source": [
"These are the actual times required, but for the purposes of discussion, let us smooth these curves by averaging over a few adjacent points / matrix dimensions. This is done below:"
]
},
{
"cell_type": "code",
"execution_count": 95,
"id": "41ac6108-526d-476b-8255-db9bb52a30ca",
"metadata": {},
"outputs": [],
"source": [
"smooth_numpy_times = []\n",
"smooth_krylov_times = []\n",
"\n",
"# Choose the number of adjacent points over which to average forward; the same will be used backward.\n",
"smooth_steps = 10\n",
"\n",
"# We will do this smoothing for all points/matrix dimensions\n",
"for i in range(len(matrix_krylov_tot_times)):\n",
" # Ensure we don't exceed the boundaries of our lists\n",
" start = max(0, i - smooth_steps)\n",
" end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)\n",
"\n",
" # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov and the NumPy calculations.\n",
" smooth_count = 0\n",
" smooth_numpy_sum = 0\n",
" smooth_krylov_sum = 0\n",
"\n",
" for j in range(start, end):\n",
" smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]\n",
" smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]\n",
" smooth_count = smooth_count + 1\n",
"\n",
" # Appending the averaged adjacent values to our new smooth lists\n",
" smooth_numpy_times.append(smooth_numpy_sum / smooth_count)\n",
" smooth_krylov_times.append(smooth_krylov_sum / smooth_count)"
]
},
{
"cell_type": "code",
"execution_count": 96,
"id": "6963ee9b-55a8-42ab-bc4f-d470bb81e543",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/6963ee9b-55a8-42ab-bc4f-d470bb81e543-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(smooth_numpy_times, color=\"blue\")\n",
"plt.plot(smooth_krylov_times, color=\"green\")\n",
"plt.xlabel(\"Dimension of matrix\") # Add x-axis label\n",
"plt.ylabel(\"Time to diagonalize (smoothed, microsec)\") # Add y-axis label\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7e1bbb91-5da7-43ea-8b99-d04ea1196e76",
"metadata": {},
"source": [
"Note that the time required for the building of a Krylov subspace initially exceeds the time required for numpy's full diagonalization. But as the size of the matrix grows, the Krylov method becomes advantageous. This is true even if we lower our acceptable error, but the advantage sets in at a larger matrix size. This is worth picking apart.\n",
"\n",
"The time complexity of numerical diagonalization is $O(n^3)$ (with some variation among algorithms). The time complexity of generating an orthonormal basis of $n$ vectors is also $O(n^3)$. So the advantage of the Krylov method is __not__ related to the use of $\\it{some}$ orthonormal basis, but to the use of a particular orthonormal basis that effectively picks out the eigenvalues of interest. We have already seen this from the sketch of a proof in the first section of this lesson, and this is critical for the convergence guarantees in Krylov methods."
]
},
{
"cell_type": "markdown",
"id": "da603cad-2e8e-4494-a0fd-f082953c35cc",
"metadata": {},
"source": [
"Let us review our progress so far:\n",
"* For very large matrices, the Krylov subspace method may yield approximate eigenvalues within required tolerances faster than traditional diagonalization algorithms.\n",
"* For such very large matrices, the generation of a Krylov subspace is the most time-consuming part of the Krylov subspace method.\n",
"* Thus an efficient way of generating a Krylov subspace would be highly valuable.\n",
"This is finally where quantum computer comes into the picture."
]
},
{
"cell_type": "markdown",
"id": "7782cac1-a86b-4a70-971f-3db5922838d7",
"metadata": {},
"source": [
"#### Check-in question\n",
"\n",
"Refer to the smoothed plot of diagonalization times vs. matrix dimension above.\n",
"(a) At approximately what matrix dimension did the Krylov method become faster, according to this plot?\n",
"(b) What aspects of the calculation could change that dimension at which the Krylov method becomes faster?\n",
"\n",
"__Answer:__\n",
"(a) Answers may vary if you re-run the calculation, but the Krylov method becomes faster at approximately dimension 50.\n",
"(b) There are many possible answers. Some important factors are the precision we require and the sparsity of the matrices being diagonalized."
]
},
{
"cell_type": "markdown",
"id": "03b62a85-4737-41ec-b264-55d22fe814bc",
"metadata": {},
"source": [
"## 3. Krylov via time evolution\n",
"\n",
"Everything that we have described so far can be done classically. So how and when would we use a quantum computer? For very large matrices, the Krylov method can require long computing times and large amounts of memory. The time required for matrix operation of $H$ on $|v\\rangle$ scales like $O(N^2)$ in the worst case. Even multiplying sparse matrices on a vector (the typical case for classical Krylov-type solvers) has a time complexity scaling like $O(N)$. This is done for every vector we want in our subspace. The subspace dimension $r$ is usually not a significant fraction of $N$, and often scales like $\\log(N)$. So generating all vectors scales like $O(N^2 \\log(N))$ in the worst case. Although there are other steps, like orthogonalization, this is the dominant scaling to keep in mind.\n",
"\n",
"Quantum computing allows us to change what attributes of the problem determine the scaling of the time and resources required. Instead of dependence on matrix size $N$ across the board, we will see things like number of shots and number of non-commuting Pauli terms that make up the Hamiltonian. Let’s explore how this works."
]
},
{
"cell_type": "markdown",
"id": "1636f90b-6bdf-4b11-aee3-9be3f2e82d79",
"metadata": {},
"source": [
"### 3.1 Time evolution\n",
"\n",
"Recall that the operator that time-evolves a quantum state is $e^{-iHt/\\hbar}$ (and it is very common, especially in quantum computing to drop the $\\hbar$ from the notation). One way of understanding and even realizing such an exponential function of an operator is to look at its Fourier series expansion. Note that this operation acting on some initial vector $|v\\rangle$ yields a sum of terms with increasing powers of $H$ applied to the initial state. It looks like we can just make our Krylov subspace by time-evolving our initial guess state!\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"e^{-iHt/\\hbar}→e^{-iHt}&≈1-iHt-\\frac{(H^2 t^2)}{2}+⋯\\\\\n",
"e^{-iHt} |v\\rangle &≈ |v\\rangle-iHt|v\\rangle-\\frac{(H^2 t^2)}{2}|v\\rangle+⋯\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"The caveat is in realizing the time evolution on a real quantum computer. Many of the terms in the Hamiltonian will not commute with each other. So while some simple exponential operators like $e^{-iZ}$ correspond to simple circuits, general Hamiltonians do not. And since they contain non-commuting terms, we can’t simply decompose the exponential into a product of simple ones, the way we can with numbers.\n",
"\n",
"$$\n",
"e^{-iHt}=e^{-i(H_1+H_2+⋯+H_n)t}\\neq e^{-iH_1 t} e^{-iH_2 t}... e^{-iH_n t}\n",
"$$\n",
"\n",
"So this is not trivial, but this is a well-studied process in quantum computing. We carry out time-evolution on quantum computers using a process called trotterization, which in itself is a rich subject[\\[10\\]](#references). But at a very high level, by breaking the time evolution into very small steps, say $m$ steps of size $dt$, we limit the effects of the non-commutativity of terms.\n",
"\n",
"$$\n",
"e^{-iHt}=e^{-i(H_1+H_2+⋯+H_n )t} = (e^{-i(H_1+H_2+⋯+H_n )t/m} )^m\n",
"≈(e^{-iH_1 dt} e^{-iH_2 dt} …e^{-iH_n dt} )^m\n",
"$$\n",
"where $dt = t/m$.\n",
"\n",
"Let us call a Krylov subspace of order r that we generated in the classical context using powers of H directly a “power Krylov subspace”.\n",
"\n",
"$$\n",
"\\mathcal{K}_P^r (H,|v\\rangle)=\\text{span}\\{|v\\rangle,H|v\\rangle,H^2 |v\\rangle… H^{r-1} |v\\rangle\\}\n",
"$$\n",
"\n",
"Now we generate a similar space using the unitary time-evolution operator $U \\equiv e^{-iHt}$; we'll refer to this as the “unitary Krylov space” $\\mathcal{K}_U^r$. The power Krylov subspace $\\mathcal{K}_P^r$ that we use classically cannot be generated directly on a quantum computer as $H$ is not a unitary operator. Using the unitary Krylov subspace can be shown to give similar convergence guarantees as the power Krylov subspace, namely, the ground state error converges efficiently as long as the initial state $|v\\rangle$ has overlap with the true ground state that is not exponentially vanishing, and as long as there is a sufficient gap between eigenvalues. See Ref [\\[1\\]](#references) for a more precise discussion of convergence.\n",
"\n",
"Here, powers of $U$ become different time steps (the $k^\\text{th}$ power of $U$ is stepping forward by a time $k \\times dt$). We can label the element of the subspace that is time-evolved for total time $k dt$ as $|\\psi_k\\rangle$."
]
},
{
"cell_type": "markdown",
"id": "619c81f7-be3c-4ffa-87f0-2ab6420381ca",
"metadata": {},
"source": [
"$$\n",
"\\begin{aligned}\n",
"U&=e^{-iHdt}\\\\\n",
"U^k&=e^{-iH(kdt)}\\\\\n",
"\\mathcal{K}_U^r&=\\text{span}\\{|\\psi\\rangle,U|\\psi\\rangle,U^2 |\\psi\\rangle… U^{r-1} |\\psi\\rangle\\}\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "ca6352b4-88c2-4dab-9f82-1750223ec8af",
"metadata": {},
"source": [
"We can project our Hamiltonian H on to the unitary Krylov subspace, $\\mathcal{K}_U^r$. In other words, we calculate each matrix element of $H$ in the $\\mathcal{K}_U^r$ basis. We'll refer to this projected matrix as $\\tilde{H}$.\n",
"\n",
"### 3.2 How to implement on a quantum computer\n",
"\n",
"The matrix elements of $\\tilde{H}$ are given by the expectation values $\\langle \\psi_m |H| \\psi_n\\rangle$, which can be estimated using the quantum computer. Keep in mind that $H$ can be written as a sum of Pauli operators on different qubits, and that not all the Pauli operators can be measured simultaneously. We can sort the Pauli terms into groups of commuting terms, and measure all those at once. But we may need many such groups to cover all the terms. So the number of distinct commuting groups into which the terms can be partitioned, $N_\\text{GCP}$ becomes important.\n",
"\n",
"$$\n",
"H=\\sum_{\\alpha=1}^{N_\\text{GCP}} c_\\alpha P_\\alpha\n",
"$$\n",
"Here $P_\\alpha$ is a Pauli string of the form $P_\\alpha \\sim IZIXII...YZXIX$ or a set of such Pauli strings that commute with one another. Given that we can write $H$ as such a sum of measureable operators, the following expressions for matrix elements of $\\tilde{H}$ can be realized using the Qiskit Runtime primitive Estimator.\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\tilde{H}_{mn}&=\\langle \\psi_m |H| \\psi_n\\rangle\\\\\n",
"&=\\langle \\psi e^{iHt_m} |H| \\psi e^{-iHt_n}\\rangle\\\\\n",
"&=\\langle \\psi e^{iHmdt} |H|\\psi e^{-iHndt}\\rangle\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Where $\\vert \\psi_n \\rangle = e^{-i H t_n} \\vert \\psi \\rangle$ are the vectors of the unitary Krylov space and $t_n = n dt$ are the multiples of the time step $dt$ chosen. On a quantum computer, the calculation of each matrix elements can be done with any algorithm which allows us to obtain overlap between quantum states. In this lesson we'll focus on the Hadamard test. Given that the $\\mathcal{K}_U$ has dimension $r$, the Hamiltonian projected into the subspace will have dimensions $r \\times r$. With $r$ small enough (generally $r<<100$ is sufficient to obtain convergence of estimates of eigenvalues) we can then easily diagonalize the projected Hamiltonian $\\tilde{H},$ classically. However, we cannot directly diagonalize $\\tilde{H}$ because of the non-orthogonality of the Krylov space vectors. We'll have to measure their overlaps and construct a matrix $\\tilde{S}$\n",
"\n",
"$$\n",
"\\tilde{S}_{mn} = \\langle \\psi_m \\vert \\psi_n \\rangle\n",
"$$\n",
"\n",
"This allows us to solve the eigenvalue problem in a non-orthogonal space (also called generalized eigenvalue problem)\n",
"\n",
"$$\n",
"\\tilde{H} \\ \\vec{c} = E \\ \\tilde{S} \\ \\vec{c}\n",
"$$\n",
"\n",
"One can then obtain estimates of the eigenvalues and eigenstates of $H$ by looking at the solutions of this generalized eigenvalue problem. For example, the estimate of the ground state energy is obtained by taking the smallest eigenvalue $E$ and the ground state from the corresponding eigenvector $\\vec{c}$. The coefficients in $\\vec{c}$ determines the contribution of the different vectors that span $\\mathcal{K}_U$."
]
},
{
"cell_type": "markdown",
"id": "7d29b607-3f2c-445c-9630-3cde4edf3299",
"metadata": {},
"source": [
"#### Generalized eigenvalue problem\n",
"\n",
"Why can we not simply diagonalize $\\tilde{H}$? Since $\\tilde{S}$ contains the information about the geometry of the Krylov basis (which is nonorthogonal in all but very special cases), $\\tilde{H}$ on its own does not describe a projection of the full Hamiltonian, so its eigenvalues have no particular relation to those of the full Hamiltonian -- they could be any random values. Solving the generalized eigenvalue problem is required to obtain the approximate eigenvalues and eigenvectors corresponding to the projection of the full Hamiltonian into the Krylov space.."
]
},
{
"cell_type": "markdown",
"id": "2f3ba07d-0c16-47fb-aeed-b4ef0a1f25e3",
"metadata": {},
"source": [
"![A circuit diagram with many layers indicating that the circuit must be used many times with different states to perform the modified Hadamard test.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig4.avif)\n",
"\n",
"The Figure shows a circuit representation of the modified Hadamard test, a method that is used to compute the overlap between different quantum states. For each matrix element $\\tilde{H}_{i,j}$, a Hadamard test between the state $\\vert \\psi_i \\rangle$, $\\vert \\psi_j \\rangle$ is carried out. This is highlighted in the figure by the color scheme for the matrix elements and the corresponding $\\text{Prep} \\; \\psi_i$, $\\text{Prep} \\; \\psi_j$ operations. Thus, a set of Hadamard tests for all the possible combinations of Krylov space vectors is required to compute all the matrix elements of the projected Hamiltonian $\\tilde{H}$. The top wire in the Hadamard test circuit is an ancilla qubit which is measured either in the X or Y basis, its expectation value determines the value of the overlap between the states. The bottom wire represents all the qubits of the system Hamiltonian. The $\\text{Prep} \\; \\psi_i$ operation prepares the system qubit in the state $\\vert \\psi_i \\rangle$ controlled by the state of the ancilla qubit (similarly for $\\text{Prep} \\; \\psi_j$) and the operation $P$ represents Pauli decomposition of the system Hamiltonian $H = \\sum_i P_i$. The implementation of this is on a quantum computer is shown in greater detail below."
]
},
{
"cell_type": "markdown",
"id": "c12c2855-f276-48a2-8901-13baf0384778",
"metadata": {},
"source": [
"## 4. Krylov quantum diagonalization on a quantum computer\n",
"\n",
"We will now implement Krylov quantum diagonalization on a real quantum computer. Let's start by importing some useful packages."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "afa233e1-1e80-4843-a958-0f84cec707ea",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pylab as plt\n",
"from typing import Union, List\n",
"import warnings\n",
"\n",
"from qiskit.quantum_info import SparsePauliOp, Pauli\n",
"from qiskit.circuit import Parameter\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"from qiskit.circuit.library import PauliEvolutionGate\n",
"from qiskit.synthesis import LieTrotter\n",
"\n",
"# from qiskit.providers.fake_provider import Fake20QV1\n",
"from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator\n",
"\n",
"import itertools as it\n",
"\n",
"warnings.filterwarnings(\"ignore\")"
]
},
{
"cell_type": "markdown",
"id": "e4a3fba7-051f-4df4-b78a-7bfeb18caf1e",
"metadata": {},
"source": [
"We define the function below to solve the generalized eigenvalue problem we just explained above."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5446eb74-126f-4db1-b018-d5f4613f79e7",
"metadata": {},
"outputs": [],
"source": [
"def solve_regularized_gen_eig(\n",
" h: np.ndarray,\n",
" s: np.ndarray,\n",
" threshold: float,\n",
" k: int = 1,\n",
" return_dimn: bool = False,\n",
") -> Union[float, List[float]]:\n",
" \"\"\"\n",
" Method for solving the generalized eigenvalue problem with regularization\n",
"\n",
" Args:\n",
" h (numpy.ndarray):\n",
" The effective representation of the matrix in our Krylov subspace\n",
" s (numpy.ndarray):\n",
" The matrix of overlaps between vectors of our Krylov subspace\n",
" threshold (float):\n",
" Cut-off value for the eigenvalue of s\n",
" k (int):\n",
" Number of eigenvalues to return\n",
" return_dimn (bool):\n",
" Whether to return the size of the regularized subspace\n",
"\n",
" Returns:\n",
" lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem\n",
"\n",
"\n",
" \"\"\"\n",
" s_vals, s_vecs = sp.linalg.eigh(s)\n",
" s_vecs = s_vecs.T\n",
" good_vecs = np.array([vec for val, vec in zip(s_vals, s_vecs) if val > threshold])\n",
" h_reg = good_vecs.conj() @ h @ good_vecs.T\n",
" s_reg = good_vecs.conj() @ s @ good_vecs.T\n",
" if k == 1:\n",
" if return_dimn:\n",
" return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)\n",
" else:\n",
" return sp.linalg.eigh(h_reg, s_reg)[0][0]\n",
" else:\n",
" if return_dimn:\n",
" return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)\n",
" else:\n",
" return sp.linalg.eigh(h_reg, s_reg)[0][:k]"
]
},
{
"cell_type": "markdown",
"id": "7dce32ec-33eb-474d-88bf-e07f6563d6a2",
"metadata": {},
"source": [
"At least in initial benchmarking, it is useful to know an exact classical solution to check convergence behavior. The function below calculates the ground state energy of a Hamiltonian, using the Hamiltonian and the number of qubits as arguments."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "f2d43bcc-5210-448d-b3d9-996796f782f6",
"metadata": {},
"outputs": [],
"source": [
"def single_particle_gs(H_op, n_qubits):\n",
" \"\"\"\n",
" Find the ground state of the single particle(excitation) sector\n",
" \"\"\"\n",
" H_x = []\n",
" for p, coeff in H_op.to_list():\n",
" H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))\n",
"\n",
" H_z = []\n",
" for p, coeff in H_op.to_list():\n",
" H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))\n",
"\n",
" H_c = H_op.coeffs\n",
"\n",
" print(\"n_sys_qubits\", n_qubits)\n",
"\n",
" n_exc = 1\n",
" sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))\n",
" print(\"n_exc\", n_exc, \", subspace dimension\", sub_dimn)\n",
"\n",
" few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)\n",
"\n",
" sparse_vecs = [\n",
" set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)\n",
" ] # list all of the possible sets of n_exc indices of 1s in n_exc-particle states\n",
"\n",
" m = 0\n",
" for i, i_set in enumerate(sparse_vecs):\n",
" for j, j_set in enumerate(sparse_vecs):\n",
" m += 1\n",
"\n",
" if len(i_set.symmetric_difference(j_set)) <= 2:\n",
" for p_x, p_z, coeff in zip(H_x, H_z, H_c):\n",
" if i_set.symmetric_difference(j_set) == p_x:\n",
" sgn = ((-1j) ** len(p_x.intersection(p_z))) * (\n",
" (-1) ** len(i_set.intersection(p_z))\n",
" )\n",
" else:\n",
" sgn = 0\n",
"\n",
" few_particle_H[i, j] += sgn * coeff\n",
"\n",
" gs_en = min(np.linalg.eigvalsh(few_particle_H))\n",
" print(\"single particle ground state energy: \", gs_en)\n",
" return gs_en"
]
},
{
"cell_type": "markdown",
"id": "c70998c4-e65c-456e-b3b4-61a289c8dd84",
"metadata": {},
"source": [
"### 4.1 Step 1: Map problem to quantum circuits and operators\n",
"\n",
"Now we will define a Hamiltonian. This is distinct from the function above in that the function above takes a Hamiltonian as an argument and returns only the ground state, and it does so classically. This Hamiltonian we define here determines the energy levels of all energy eigenstates, and this Hamiltonian can be constructed using Pauli operators and implemented on a quantum computer.\n",
"\n",
"We choose a Hamiltonian corresponding to a chain of spins which can have any orientation in space, called a \"Heisenberg chain\". We assume that the $i^\\text{th}$ spin can be influenced by its nearest neighbors (the $(i-1)^\\text{th}$ and $(i+1)^\\text{th}$ spins) but not by more distant neighbors. We also allow for the possibility that the interaction between spins is different when the spins point along different axes. This asymmetry sometimes arises, for example, due to the structure of the crystal lattice into which spins are embedded."
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "8b547f8a-df47-4e56-921b-3955eb7c19a9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[('ZZIIIIIIII', 2), ('IZZIIIIIII', 2), ('IIZZIIIIII', 2), ('IIIZZIIIII', 2), ('IIIIZZIIII', 2), ('IIIIIZZIII', 2), ('IIIIIIZZII', 2), ('IIIIIIIZZI', 2), ('IIIIIIIIZZ', 2), ('XXIIIIIIII', 1), ('IXXIIIIIII', 1), ('IIXXIIIIII', 1), ('IIIXXIIIII', 1), ('IIIIXXIIII', 1), ('IIIIIXXIII', 1), ('IIIIIIXXII', 1), ('IIIIIIIXXI', 1), ('IIIIIIIIXX', 1), ('YYIIIIIIII', 3), ('IYYIIIIIII', 3), ('IIYYIIIIII', 3), ('IIIYYIIIII', 3), ('IIIIYYIIII', 3), ('IIIIIYYIII', 3), ('IIIIIIYYII', 3), ('IIIIIIIYYI', 3), ('IIIIIIIIYY', 3)]\n"
]
}
],
"source": [
"# Define problem Hamiltonian.\n",
"n_qubits = 10\n",
"# coupling strength for XX, YY, and ZZ interactions\n",
"JX = 1\n",
"JY = 3\n",
"JZ = 2\n",
"\n",
"# Define the Hamiltonian:\n",
"H_int = [[\"I\"] * n_qubits for _ in range(3 * (n_qubits - 1))]\n",
"for i in range(n_qubits - 1):\n",
" H_int[i][i] = \"Z\"\n",
" H_int[i][i + 1] = \"Z\"\n",
"for i in range(n_qubits - 1):\n",
" H_int[n_qubits - 1 + i][i] = \"X\"\n",
" H_int[n_qubits - 1 + i][i + 1] = \"X\"\n",
"for i in range(n_qubits - 1):\n",
" H_int[2 * (n_qubits - 1) + i][i] = \"Y\"\n",
" H_int[2 * (n_qubits - 1) + i][i + 1] = \"Y\"\n",
"H_int = [\"\".join(term) for term in H_int]\n",
"H_tot = [\n",
" (term, JZ)\n",
" if term.count(\"Z\") == 2\n",
" else (term, JY)\n",
" if term.count(\"Y\") == 2\n",
" else (term, JX)\n",
" for term in H_int\n",
"]\n",
"\n",
"# Get operator\n",
"H_op = SparsePauliOp.from_list(H_tot)\n",
"print(H_tot)"
]
},
{
"cell_type": "markdown",
"id": "759b8b55-894b-4d3b-93de-5df4f63f9609",
"metadata": {},
"source": [
"The code below restricts the Hamiltonian to single particle states, and uses the spectral norm to set a good size for our time step $dt$. We heuristically choose a value for the time-step `dt` (based on upper bounds on the Hamiltonian norm). Ref [\\[9\\]](#references) showed that a sufficiently small timestep is $\\pi/\\vert \\vert H \\vert \\vert$, and that it is preferable up to a point to underestimate this value rather than overestimate, since overestimating can allow contributions from high-energy states to corrupt even the optimal state in the Krylov space. On the other hand, choosing $dt$ to be too small leads to worse conditioning of the Krylov subspace, since the Krylov basis vectors differ less from timestep to timestep."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "3dd96fbe-47bb-444b-b403-c67c2dcc9d07",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.17453292519943295"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Get Hamiltonian restricted to single-particle states\n",
"single_particle_H = np.zeros((n_qubits, n_qubits))\n",
"for i in range(n_qubits):\n",
" for j in range(i + 1):\n",
" for p, coeff in H_op.to_list():\n",
" p_x = Pauli(p).x\n",
" p_z = Pauli(p).z\n",
" if all(p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)):\n",
" sgn = ((-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))) * (\n",
" (-1) ** p_z[i]\n",
" )\n",
" else:\n",
" sgn = 0\n",
" single_particle_H[i, j] += sgn * coeff\n",
"for i in range(n_qubits):\n",
" for j in range(i + 1, n_qubits):\n",
" single_particle_H[i, j] = np.conj(single_particle_H[j, i])\n",
"\n",
"# Set dt according to spectral norm\n",
"dt = np.pi / np.linalg.norm(single_particle_H, ord=2)\n",
"dt"
]
},
{
"cell_type": "markdown",
"id": "cdb8c08b-caca-4bd5-ae93-065c42485e13",
"metadata": {},
"source": [
"Let us specify the maximum Krylov dimension we wish to use, though we will check for convergence at smaller dimensions. We also specify the number of Trotter steps to use in the time evolution. In this case, we set the Krylov dimension to be equal to the number of Trotter steps plus one. This is consistent with generating an additional Krylov vector with each application of $U$, plus the initial vector. For the sake of this lesson, we will choose a small Krylov dimension of just 5. This is quite limited in the context of realistic applications, but it is sufficient for this example. We will explore methods in later lessons that allow us to scale and project our Hamiltonians onto larger subspaces."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "d95c6f17-7275-474a-b1e5-a41c4539ed83",
"metadata": {},
"outputs": [],
"source": [
"# Set parameters for quantum Krylov algorithm\n",
"krylov_dim = 5 # size of krylov subspace\n",
"num_trotter_steps = 6\n",
"dt_circ = dt / num_trotter_steps"
]
},
{
"cell_type": "markdown",
"id": "29291306-28a3-4c12-94ef-b3ffd5521c76",
"metadata": {},
"source": [
"#### State preparation\n",
"\n",
"Pick a reference state $\\vert \\psi \\rangle$ that has some overlap with the ground state. For this Hamiltonian, We use the a state with an excitation in the middle qubit $\\vert 00..010...00 \\rangle$ as our reference state."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "410192fb-8197-4860-8c3a-2e874e2f9c56",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/410192fb-8197-4860-8c3a-2e874e2f9c56-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qc_state_prep = QuantumCircuit(n_qubits)\n",
"qc_state_prep.x(int(n_qubits / 2) + 1)\n",
"qc_state_prep.draw(\"mpl\", scale=0.5)"
]
},
{
"cell_type": "markdown",
"id": "dee3c383-fc2e-4420-8d04-861f5f4f4576",
"metadata": {},
"source": [
"#### Time evolution\n",
"\n",
"We can realize the time-evolution operator generated by a given Hamiltonian: $U=e^{-iHt}$ via the [Lie-Trotter approximation](https://docs.quantum.ibm.com/api/qiskit/qiskit.synthesis.LieTrotter). For simplicity we use the built-in ```PauliEvolutionGate``` in the time-evolution circuit. The general syntax for this is this."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "4b3dddd0-391a-412a-9663-9e15a153125a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<qiskit.circuit.instructionset.InstructionSet at 0x7f486e895900>"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Parameter(\"t\")\n",
"\n",
"## Create the time-evo op circuit\n",
"evol_gate = PauliEvolutionGate(\n",
" H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)\n",
")\n",
"\n",
"qr = QuantumRegister(n_qubits)\n",
"qc_evol = QuantumCircuit(qr)\n",
"qc_evol.append(evol_gate, qargs=qr)"
]
},
{
"cell_type": "markdown",
"id": "295bf203-2bc0-4c02-912b-423e5ca392bb",
"metadata": {},
"source": [
"We will use a version of this below in the Hadamard test, but stepping forward for times $dt$."
]
},
{
"cell_type": "markdown",
"id": "03f71981-3a22-4ba6-8281-7b20895dbda3",
"metadata": {},
"source": [
"#### Hadamard test\n",
"\n",
"Recall that we wish to calculate the matrix elements of both $\\tilde{H}$ and the Gram matrix $\\tilde{S}$ using the Hadamard test. Let's review how this works in this context, focusing first on the construction of $\\tilde{H}.$ The overall process is depicted graphically below. The layers of colored state preparation blocks $\\text{Prep}|\\psi_i\\rangle$ serve as a reminder that this process is carried out for all combinations of $|\\psi_i\\rangle$ and $|\\psi_j\\rangle$ in our subspace.\n",
"\n",
"![An image of a quantum circuit diagram with many layers indicating that the circuit must be evaluated for many different states in order to perform the Hadamard test.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig5.avif)\n",
"\n",
"The states of the system at the steps indicated are:\n",
"$$\n",
"\\begin{aligned}\n",
" \\text{Step 0:}\\qquad|\\Psi\\rangle & = |0\\rangle|0\\rangle^N \\\\\n",
" \\text{Step 1:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle + |1\\rangle \\Big)|0\\rangle^N \\\\\n",
" \\text{Step 2:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle|0\\rangle^N+|1\\rangle |\\psi_i\\rangle\\Big)\\\\\n",
" \\text{Step 3:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |0\\rangle^N+|1\\rangle P |\\psi_i\\rangle\\Big) \\\\\n",
" \\text{Step 4:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |\\psi_j\\rangle+|1\\rangle P|\\psi_i\\rangle\\Big)\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Here $P$ is a Pauli term in the decomposition of the Hamiltonian (note that it cannot be a linear combination of multiple commuting Pauli terms since that would not be unitary -- grouping is possible using a different construction we will show later) $\\text{Prep} \\; \\psi_i$, $\\text{Prep} \\; \\psi_j$ are controlled operations that prepare $|\\psi_i\\rangle$, $|\\psi_j\\rangle$ vectors of the unitary Krylov space, with $|\\psi_k\\rangle = e^{-i H k dt } \\vert \\psi \\rangle = e^{-i H k dt } U_{\\psi} \\vert 0 \\rangle^N$. Applying measurements of $X$ and $Y$ to this circuit calculates the real and imaginary parts, respectively, of the matrix elements we require.\n",
"\n",
"Starting from Step 4 above, apply the Hadamard gate $H$ to the zeroth qubit.\n",
"\n",
"$$\n",
"\\begin{equation*}\n",
" |\\Psi\\rangle \\longrightarrow\\quad\\frac{1}{2}|0\\rangle\\Big( |\\psi_j\\rangle + P|\\psi_i\\rangle\\Big) + \\frac{1}{2}|1\\rangle\\Big(|\\psi_j\\rangle - P|\\psi_i\\rangle\\Big)\n",
"\\end{equation*}\n",
"$$\n",
"\n",
"Then measure either $X$ or $Y$.\n",
"\n",
"$$\n",
"\\begin{equation*}\n",
"\\begin{split}\n",
" \\Rightarrow\\quad\\langle X\\rangle &= \\frac{1}{4}\\Bigg(\\Big\\|| \\psi_j\\rangle + P|\\psi_i\\rangle \\Big\\|^2-\\Big\\||\\psi_j\\rangle - P|\\psi_i\\rangle\\Big\\|^2\\Bigg) \\\\\n",
" &= \\text{Re}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n",
"\\end{split}\n",
"\\end{equation*}\n",
"$$\n",
"\n",
"From the identity $|a + b\\|^2 = \\langle a + b | a + b \\rangle = \\|a\\|^2 + \\|b\\|^2 + 2\\text{Re}\\langle a | b \\rangle$. Similarly, measuring $Y$ yields\n",
"\n",
"$$\n",
"\\begin{equation*}\n",
" \\langle Y\\rangle = \\text{Im}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n",
"\\end{equation*}\n",
"$$\n",
"\n",
"Adding these steps to the time-evolution we set up previously we write the following."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "ac536079-96f2-435b-a15f-72c207fa24d1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Circuit for calculating the real part of the overlap in S via Hadamard test\n"
]
},
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/ac536079-96f2-435b-a15f-72c207fa24d1-1.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Create the time-evo op circuit\n",
"evol_gate = PauliEvolutionGate(\n",
" H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n",
")\n",
"\n",
"## Create the time-evo op dagger circuit\n",
"evol_gate_d = PauliEvolutionGate(\n",
" H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n",
")\n",
"evol_gate_d = evol_gate_d.inverse()\n",
"\n",
"# Put pieces together\n",
"qc_reg = QuantumRegister(n_qubits)\n",
"qc_temp = QuantumCircuit(qc_reg)\n",
"qc_temp.compose(qc_state_prep, inplace=True)\n",
"for _ in range(num_trotter_steps):\n",
" qc_temp.append(evol_gate, qargs=qc_reg)\n",
"for _ in range(num_trotter_steps):\n",
" qc_temp.append(evol_gate_d, qargs=qc_reg)\n",
"qc_temp.compose(qc_state_prep.inverse(), inplace=True)\n",
"\n",
"# Create controlled version of the circuit\n",
"controlled_U = qc_temp.to_gate().control(1)\n",
"\n",
"# Create hadamard test circuit for real part\n",
"qr = QuantumRegister(n_qubits + 1)\n",
"qc_real = QuantumCircuit(qr)\n",
"qc_real.h(0)\n",
"qc_real.append(controlled_U, list(range(n_qubits + 1)))\n",
"qc_real.h(0)\n",
"\n",
"print(\"Circuit for calculating the real part of the overlap in S via Hadamard test\")\n",
"qc_real.draw(\"mpl\", fold=-1, scale=0.5)"
]
},
{
"cell_type": "markdown",
"id": "b88673a4-859a-44f7-ac0c-38ecf17faf36",
"metadata": {},
"source": [
"We warned already about the depth involved in Trotter circuits. Performing the Hadamard test in these conditions can yield an even deeper circuit, especially once we decompose to native gates. This will increase even more if we account for the topology of the device. So before using any time on the quantum computer, it is a good idea to check the 2-qubit depth of our circuit."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f40d206a-9bd5-4355-8007-cecff21f7fbe",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of layers of 2Q operations 34993\n"
]
}
],
"source": [
"print(\n",
" \"Number of layers of 2Q operations\",\n",
" qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),\n",
")"
]
},
{
"cell_type": "markdown",
"id": "bbd4f53b-8946-4713-89c2-e6bbbfa80609",
"metadata": {},
"source": [
"A circuit of this depth cannot return usable results on modern quantum computers. If we are to construct $\\tilde{H}$ and $\\tilde{S},$ we need a better way. This is the reason for the efficient Hadamard test introduced below."
]
},
{
"cell_type": "markdown",
"id": "6cda7f15-7c76-40aa-bf2e-bfea141c8474",
"metadata": {},
"source": [
"### 4. 2 Step 2. Optimize circuits and operators for target hardware\n",
"\n",
"#### Efficient Hadamard test\n",
"\n",
"We can optimize the deep circuits for the Hadamard test that we have obtained by introducing some approximations and relying on some assumption about the model Hamiltonian. For example, consider the following circuit for the Hadamard test:\n",
"\n",
"![An image of a quantum circuit diagram with many layers indicating that the circuit must be evaluated for many different unitary operators in order to perform the modified, efficient Hadamard test.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig6.avif)\n",
"\n",
"Assume we can classically calculate $E_0$, the eigenvalue of $|0\\rangle^N$ under the Hamiltonian $H$.\n",
"This is satisfied when the Hamiltonian preserves the U(1) symmetry. Although this may seem like a strong assumption, there are many cases where it is safe to assume that there is a vacuum state (in this case it maps to the $|0\\rangle^N$ state) which is unaffected by the action of the Hamiltonian. This is true for example for chemistry Hamiltonians that describe stable molecule (where the number of electrons is conserved).\n",
"Given that the gate $\\text{Prep} \\; \\psi_0$, prepares the desired reference state $\\ket{\\psi_0} = \\text{Prep} \\; \\psi_0 \\ket{0} = e^{-i H 0 dt} U_{\\psi_0} \\ket{0}$, e.g., to prepare the HF state for chemistry $\\text{Prep} \\; \\psi_0$ would be a product of single-qubit NOTs, so controlled-$\\text{Prep} \\; \\psi_0$ is just a product of CNOTs.\n",
"Then the circuit above implements the following state prior to measurement:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\text{Step 0:}\\qquad|\\Psi\\rangle & = \\ket{0} \\ket{0}^{N}\\\\\n",
" \\text{Step 1:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\left(\\ket{0}\\ket{0}^N+ \\ket{1} \\ket{0}^N\\right)\\\\\n",
" \\text{Step 2:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\left(|0\\rangle|0\\rangle^N+|1\\rangle|\\psi_0\\rangle\\right)\\\\\n",
" \\text{Step 3:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\left(e^{i\\phi}\\ket{0}\\ket{0}^N+\\ket{1} U\\ket{\\psi_0}\\right)\\\\\n",
" \\text{Step 4:}\\qquad|\\Psi\\rangle & = \\frac{1}{\\sqrt{2}}\\left(e^{i\\phi}\\ket{0} \\ket{\\psi_0}+\\ket{1} U\\ket{\\psi_0}\\right)\\\\\n",
" & = \\frac{1}{2}\\left(\\ket{+}\\left(e^{i\\phi}\\ket{\\psi_0}+U\\ket{\\psi_0}\\right)+\\ket{-}\\left(e^{i\\phi}\\ket{\\psi_0}-U\\ket{\\psi_0}\\right)\\right)\\\\\n",
" & = \\frac{1}{2}\\left(\\ket{+i}\\left(e^{i\\phi}\\ket{\\psi_0}-iU\\ket{\\psi_0}\\right)+\\ket{-i}\\left(e^{i\\phi}\\ket{\\psi_0}+iU\\ket{\\psi_0}\\right)\\right)\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"where we have used the classical simulable phase shift $ U\\ket{0}^N = e^{i\\phi}\\ket{0}^N$ from step 2 to 3. Therefore the expectation values are\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\langle X\\otimes P\\rangle&=\\frac{1}{4}\n",
" \\Big(\n",
" \\left(e^{-i\\phi}\\bra{\\psi_0}+\\bra{\\psi_0}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi_0}+U\\ket{\\psi_0}\\right)\n",
" \\\\\n",
" &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi_0}-\\bra{\\psi_0}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi_0}-U\\ket{\\psi_0}\\right)\n",
" \\Big)\\\\\n",
" &=\\text{Re}\\left[e^{-i\\phi}\\bra{\\psi_0}PU\\ket{\\psi_0}\\right],\n",
"\\end{aligned}\n",
"\n",
"$$\n",
"\n",
"$$\n",
"\n",
"\\begin{aligned}\n",
" \\langle Y\\otimes P\\rangle&=\\frac{1}{4}\n",
" \\Big(\n",
" \\left(e^{-i\\phi}\\bra{\\psi_0}+i\\bra{\\psi_0}U^\\dagger\\right)P\\left(e^{i\\phi_0}\\ket{\\psi_0}-iU\\ket{\\psi_0}\\right)\n",
" \\\\\n",
" &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi_0}-i\\bra{\\psi_0}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi_0}+iU\\ket{\\psi_0}\\right)\n",
" \\Big)\\\\\n",
" &=\\text{Im}\\left[e^{-i\\phi}\\bra{\\psi_0}PU\\ket{\\psi_0}\\right].\n",
"\\end{aligned}\n",
"\n",
"$$\n",
"\n",
"Using these assumptions we were able to write the expectation values of operators of interest with fewer controlled operations. In fact, we only need to implement the controlled state preparation $\\text{Prep} \\; \\psi_0$ and not controlled time evolutions. Reframing our calculation as above will allow us to greatly reduce the depth of the resulting circuits.\n",
"\n",
"Note that as a bonus, since the Pauli operator now appears as a measurement at the end of the circuit rather than as a controlled gate in the middle, it can be measured alongside other commuting Pauli operators as in the decomposition $$H=\\sum_{\\alpha = 1}^{N_\\text{GCP}}c_\\alpha P_\\alpha $$ given above."
]
},
{
"cell_type": "markdown",
"id": "cd24795e-3ef0-4629-87bf-8e04c471d665",
"metadata": {},
"source": [
"### Decompose time-evolution operator with Trotter decomposition\n",
"\n",
"Instead of implementing the time-evolution operator exactly we can use the Trotter decomposition to implement an approximation of it. Repeating several times a certain order Trotter decomposition gives us further reduction of the error introduced from the approximation. In the following, we directly build the Trotter implementation in the most efficient way for the interaction graph of the Hamiltonian we are considering (nearest neighbor interactions only). In practice we insert Pauli rotations $R_{xx}$, $R_{yy}$, $R_{zz}$ with a parametrized angle $t$ which correspond to the approximate implementation of $e^{-i (XX + YY + ZZ) t}$. Given the difference in definition of the Pauli rotations and the time-evolution that we are trying to implement, we'll have to use the parameter $2*dt$ to achieve a time-evolution of $dt$. Furthermore, we reverse the order of the operations for odd number of repetitions of the Trotter steps, which is functionally equivalent but allows for synthesizing adjacent operations in a single $SU(2)$ unitary. This gives a much shallower circuit than what is obtained using the generic `PauliEvolutionGate()` functionality."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "8d35e6cd-c0f6-4871-83da-86472d66be2c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/8d35e6cd-c0f6-4871-83da-86472d66be2c-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Parameter(\"t\")\n",
"\n",
"# Create instruction for rotation about XX+YY-ZZ:\n",
"Rxyz_circ = QuantumCircuit(2)\n",
"Rxyz_circ.rxx(2 * t, 0, 1)\n",
"Rxyz_circ.ryy(2 * t, 0, 1)\n",
"Rxyz_circ.rzz(-2 * t, 0, 1)\n",
"Rxyz_instr = Rxyz_circ.to_instruction(label=\"RXX+YY-ZZ\")\n",
"\n",
"interaction_list = [\n",
" [[i, i + 1] for i in range(0, n_qubits - 1, 2)],\n",
" [[i, i + 1] for i in range(1, n_qubits - 1, 2)],\n",
"] # linear chain\n",
"\n",
"qr = QuantumRegister(n_qubits)\n",
"trotter_step_circ = QuantumCircuit(qr)\n",
"for i, color in enumerate(interaction_list):\n",
" for interaction in color:\n",
" trotter_step_circ.append(Rxyz_instr, interaction)\n",
" if i < len(interaction_list) - 1:\n",
" trotter_step_circ.barrier()\n",
"reverse_trotter_step_circ = trotter_step_circ.reverse_ops()\n",
"\n",
"qc_evol = QuantumCircuit(qr)\n",
"for step in range(num_trotter_steps):\n",
" if step % 2 == 0:\n",
" qc_evol = qc_evol.compose(trotter_step_circ)\n",
" else:\n",
" qc_evol = qc_evol.compose(reverse_trotter_step_circ)\n",
"\n",
"qc_evol.decompose().draw(\"mpl\", fold=-1, scale=0.5)"
]
},
{
"cell_type": "markdown",
"id": "2a3d8d4a-fc0f-4702-ac69-67f63d0c0e30",
"metadata": {},
"source": [
"We prepare an initial state again for this efficient Hadamard test."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "d1d0b9de-65d4-4a46-975d-6cfaaac05f9a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/d1d0b9de-65d4-4a46-975d-6cfaaac05f9a-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"control = 0\n",
"excitation = int(n_qubits / 2) + 1\n",
"controlled_state_prep = QuantumCircuit(n_qubits + 1)\n",
"controlled_state_prep.cx(control, excitation)\n",
"controlled_state_prep.draw(\"mpl\", fold=-1, scale=0.5)"
]
},
{
"cell_type": "markdown",
"id": "415e94f8-6594-42ab-b55b-80a34852b943",
"metadata": {},
"source": [
"#### Template circuits for calculating matrix elements of $\\tilde{S}$ and $\\tilde{H}$ via Hadamard test\n",
"\n",
"The only difference between the circuits used in the Hadamard test will be the phase in the time-evolution operator and the observables measured. Therefore we can prepare a template circuit which represent the generic circuit for the Hadamard test, with placeholders for the gates that depend on the time-evolution operator."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "27a54efa-affb-41bc-8523-822f8d92d34f",
"metadata": {},
"outputs": [],
"source": [
"# Parameters for the template circuits\n",
"parameters = []\n",
"for idx in range(1, krylov_dim):\n",
" parameters.append(dt_circ * (idx))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "37382668-1999-4475-b50f-2887449d5c93",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/37382668-1999-4475-b50f-2887449d5c93-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Create modified hadamard test circuit\n",
"qr = QuantumRegister(n_qubits + 1)\n",
"qc = QuantumCircuit(qr)\n",
"qc.h(0)\n",
"qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)\n",
"qc.barrier()\n",
"qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)\n",
"qc.barrier()\n",
"qc.x(0)\n",
"qc.compose(controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True)\n",
"qc.x(0)\n",
"\n",
"qc.decompose().draw(\"mpl\", fold=-1)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "4ad71f04-ec89-473b-9b7a-db4d3936c85e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The optimized circuit has 2Q gates depth: 74\n"
]
}
],
"source": [
"print(\n",
" \"The optimized circuit has 2Q gates depth: \",\n",
" qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),\n",
")"
]
},
{
"cell_type": "markdown",
"id": "6f60d587-8ae1-4daa-ab76-a032c2221ab7",
"metadata": {},
"source": [
"This depth is substantially reduced compared to the original Hadamard test. This depth is manageable on modern quantum computers, though it is still quite high. We will need to use state-of-the-art error mitigation to obtain useful results.\n",
"\n",
"Select a backend on which to run our quantum Krylov calculation, so that we can transpile our circuit for running on that quantum computer."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "5683752b-d520-4cad-9fe9-7ebb7f495ba0",
"metadata": {},
"outputs": [],
"source": [
"# from qiskit_ibm_runtime import QiskitRuntimeService, Batch\n",
"\n",
"# service = QiskitRuntimeService(channel = 'ibm_quantum')\n",
"# service = QiskitRuntimeService()\n",
"# backend_id = 'ibm_cusco'\n",
"# backend_id = 'test_eagle_us-east'\n",
"# backend = service.backend(backend_id)\n",
"\n",
"from qiskit_ibm_runtime import Batch\n",
"\n",
"service = QiskitRuntimeService()\n",
"backend_id = \"ibm_torino\"\n",
"backend = service.backend(backend_id)"
]
},
{
"cell_type": "markdown",
"id": "0337f370-edff-4dcb-9570-26818ee51d3a",
"metadata": {},
"source": [
"We now transpile our circuits and operators."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "8a7c9ea1-5bc3-40b7-aa55-040cbb5f8c28",
"metadata": {},
"outputs": [],
"source": [
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"target = backend.target\n",
"basis_gates = list(target.operation_names)\n",
"pm = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, basis_gates=basis_gates\n",
")\n",
"\n",
"qc_trans = pm.run(qc)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "054bcfba-fde3-4213-bbd7-d7a6dd6c03c8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"52\n",
"OrderedDict({'rz': 630, 'sx': 521, 'cz': 228, 'x': 14, 'barrier': 8})\n"
]
},
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/054bcfba-fde3-4213-bbd7-d7a6dd6c03c8-1.avif\" alt=\"Output of the previous code cell\" />"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(qc_trans.depth(lambda x: x[0].num_qubits == 2))\n",
"print(qc_trans.count_ops())\n",
"qc_trans.draw(\"mpl\", fold=-1, idle_wires=False, scale=0.5)"
]
},
{
"cell_type": "markdown",
"id": "5395cc4a-ce59-4f8f-a126-1f8ca5507f83",
"metadata": {},
"source": [
"After optimization, our transpiled two-qubit depth is further reduced.\n",
"\n",
"### 4.3 Step 3. Execute using a Qiskit Runtime primitive\n",
"\n",
"We now create PUBs for execution with Estimator."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "5d949e77-d7af-47aa-91e1-40fcf5940996",
"metadata": {},
"outputs": [],
"source": [
"# Define observables to measure for S\n",
"observable_S_real = \"I\" * (n_qubits) + \"X\"\n",
"observable_S_imag = \"I\" * (n_qubits) + \"Y\"\n",
"\n",
"observable_op_real = SparsePauliOp(\n",
" observable_S_real\n",
") # define a sparse pauli operator for the observable\n",
"observable_op_imag = SparsePauliOp(observable_S_imag)\n",
"\n",
"layout = qc_trans.layout # get layout of transpiled circuit\n",
"observable_op_real = observable_op_real.apply_layout(\n",
" layout\n",
") # apply physical layout to the observable\n",
"observable_op_imag = observable_op_imag.apply_layout(layout)\n",
"observable_S_real = (\n",
" observable_op_real.paulis.to_labels()\n",
") # get the label of the physical observable\n",
"observable_S_imag = observable_op_imag.paulis.to_labels()\n",
"\n",
"observables_S = [[observable_S_real], [observable_S_imag]]\n",
"\n",
"\n",
"# Define observables to measure for H\n",
"# Hamiltonian terms to measure\n",
"observable_list = []\n",
"for pauli, coeff in zip(H_op.paulis, H_op.coeffs):\n",
" # print(pauli)\n",
" observable_H_real = pauli[::-1].to_label() + \"X\"\n",
" observable_H_imag = pauli[::-1].to_label() + \"Y\"\n",
" observable_list.append([observable_H_real])\n",
" observable_list.append([observable_H_imag])\n",
"\n",
"layout = qc_trans.layout\n",
"\n",
"observable_trans_list = []\n",
"for observable in observable_list:\n",
" observable_op = SparsePauliOp(observable)\n",
" observable_op = observable_op.apply_layout(layout)\n",
" observable_trans_list.append([observable_op.paulis.to_labels()])\n",
"\n",
"observables_H = observable_trans_list\n",
"\n",
"\n",
"# Define a sweep over parameter values\n",
"params = np.vstack(parameters).T\n",
"\n",
"\n",
"# Estimate the expectation value for all combinations of\n",
"# observables and parameter values, where the pub result will have\n",
"# shape (# observables, # parameter values).\n",
"pub = (qc_trans, observables_S + observables_H, params)"
]
},
{
"cell_type": "markdown",
"id": "e0f30365-57d2-4ece-ba62-4e534a810de6",
"metadata": {},
"source": [
"Circuits for $t=0$ are classically calculable. We carry this out before moving on to the $t\\neq 0$ case using a quantum computer."
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "298c43a4-e54c-4ac3-95bd-12f535d44790",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(10+0j)\n"
]
}
],
"source": [
"from qiskit.quantum_info import StabilizerState, Pauli\n",
"\n",
"\n",
"qc_cliff = qc.assign_parameters({t: 0})\n",
"\n",
"\n",
"# Get expectation values from experiment\n",
"S_expval_real = StabilizerState(qc_cliff).expectation_value(\n",
" Pauli(\"I\" * (n_qubits) + \"X\")\n",
")\n",
"S_expval_imag = StabilizerState(qc_cliff).expectation_value(\n",
" Pauli(\"I\" * (n_qubits) + \"Y\")\n",
")\n",
"\n",
"# Get expectation values\n",
"S_expval = S_expval_real + 1j * S_expval_imag\n",
"\n",
"H_expval = 0\n",
"for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n",
" # Get expectation values from experiment\n",
" expval_real = StabilizerState(qc_cliff).expectation_value(\n",
" Pauli(pauli[::-1].to_label() + \"X\")\n",
" )\n",
" expval_imag = StabilizerState(qc_cliff).expectation_value(\n",
" Pauli(pauli[::-1].to_label() + \"Y\")\n",
" )\n",
" expval = expval_real + 1j * expval_imag\n",
"\n",
" # Fill-in matrix elements\n",
" H_expval += coeff * expval\n",
"\n",
"\n",
"print(H_expval)"
]
},
{
"cell_type": "markdown",
"id": "9b86af1c-7e37-427d-b4f0-fbca38d45a21",
"metadata": {},
"source": [
"Although we were able to reduce our gate depth by orders of magnitude using the efficient Hadamard test, the depth is still sufficient to require state-of-the-art error mitigation. Below, we specify attributes of the mitigation being used. All of the methods used are important, but it is worth called out [probabilistic error amplification (PEA)](https://docs.quantum.ibm.com/guides/error-mitigation-and-suppression-techniques#probabilistic-error-amplification-pea) specifically. This powerful technique comes with a great deal of quantum overhead. The calculation done here can take 20 minutes or more to run on a real quantum computer. You may wish to play with the parameters below to increase or decrease precision and consequentially overhead. The default settings below yield high-fidelity results."
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "2dff41e1-d417-4af7-aa6c-537a9b0e0c7c",
"metadata": {},
"outputs": [],
"source": [
"# Experiment options\n",
"num_randomizations = 300\n",
"num_randomizations_learning = 20\n",
"max_batch_circuits = 20\n",
"shots_per_randomization = 100\n",
"learning_pair_depths = [0, 4, 24]\n",
"noise_factors = [1, 1.3, 1.6]\n",
"\n",
"# Base option formatting\n",
"options = {\n",
" # Builtin resilience settings for ZNE\n",
" \"resilience\": {\n",
" \"measure_mitigation\": True,\n",
" \"zne_mitigation\": True,\n",
" \"zne\": {\"noise_factors\": noise_factors},\n",
" # TREX noise learning configuration\n",
" \"measure_noise_learning\": {\n",
" \"num_randomizations\": num_randomizations_learning,\n",
" \"shots_per_randomization\": shots_per_randomization,\n",
" },\n",
" # PEA noise model configuration\n",
" \"layer_noise_learning\": {\n",
" \"max_layers_to_learn\": 10,\n",
" \"layer_pair_depths\": learning_pair_depths,\n",
" \"shots_per_randomization\": shots_per_randomization,\n",
" \"num_randomizations\": num_randomizations_learning,\n",
" },\n",
" },\n",
" # Randomization configuration\n",
" \"twirling\": {\n",
" \"num_randomizations\": num_randomizations,\n",
" \"shots_per_randomization\": shots_per_randomization,\n",
" \"strategy\": \"all\",\n",
" },\n",
" # Experimental settings for PEA method\n",
" \"experimental\": {\n",
" # # Just in case, disable any further qiskit transpilation not related to twirling / DD\n",
" # \"skip_transpilation\": True,\n",
" # Execution configuration\n",
" \"execution\": {\n",
" \"max_pubs_per_batch_job\": max_batch_circuits,\n",
" \"fast_parametric_update\": True,\n",
" },\n",
" # Error Mitigation configuration\n",
" \"resilience\": {\n",
" # ZNE Configuration\n",
" \"zne\": {\n",
" \"amplifier\": \"pea\",\n",
" \"return_all_extrapolated\": True,\n",
" \"return_unextrapolated\": True,\n",
" \"extrapolated_noise_factors\": [0] + noise_factors,\n",
" }\n",
" },\n",
" },\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "5932dc29-7905-4f54-9c97-2bb0fc9c1e39",
"metadata": {},
"source": [
"Finally, we execute the circuits for $\\tilde{S}$ and $\\tilde{H}$ with Estimator."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "01a41068-453f-4ff6-9664-c73be1965dc7",
"metadata": {},
"outputs": [],
"source": [
"# This job required 17 minutes of QPU time to run on ibm_torino. This is only an estimate. Your execution time may vary.\n",
"\n",
"with Batch(backend=backend) as batch:\n",
" # Estimator\n",
" estimator = Estimator(mode=batch, options=options)\n",
"\n",
" job = estimator.run([pub], precision=1)"
]
},
{
"cell_type": "markdown",
"id": "936865a2-828d-4a45-987e-b62cce0535da",
"metadata": {},
"source": [
"### 4.4 Step 4. Post-process and analyze results\n",
"\n",
"What we have obtained from the quantum computer are the individual matrix elements of $\\tilde{S}$ and the commuting Pauli groups that make up the matrix elements of $\\tilde{H}$. These terms must be combined to recover our matrices, so that we can solve the generalized eigenvalue problem."
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "28ed9319-dd36-4104-aba0-f8798cbbd2b6",
"metadata": {},
"outputs": [],
"source": [
"# Store the outputs as 'results'.\n",
"results = job.result()[0]"
]
},
{
"cell_type": "markdown",
"id": "1d3c0af8-fb06-415c-b591-f8e8faedc070",
"metadata": {},
"source": [
"#### Calculate Effective Hamiltonian and Overlap matrices\n",
"\n",
"First calculate the phase accumulated by the $\\vert 0 \\rangle$ state during the uncontrolled time evolution"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "d36e8b32-621d-44da-808d-297c0b60754a",
"metadata": {},
"outputs": [],
"source": [
"prefactors = [\n",
" np.exp(-1j * sum([c for p, c in H_op.to_list() if \"Z\" in p]) * i * dt)\n",
" for i in range(1, krylov_dim)\n",
"]"
]
},
{
"cell_type": "markdown",
"id": "f20402f5-1264-4eeb-9d41-8dd6e53e82d3",
"metadata": {},
"source": [
"Once we have the results of the circuit executions we can post-process the data to calculate the matrix elements of $S$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "d1ff8971-2f80-4130-b7f2-e95c3a1b476d",
"metadata": {},
"outputs": [],
"source": [
"# Assemble S, the overlap matrix of dimension D:\n",
"S_first_row = np.zeros(krylov_dim, dtype=complex)\n",
"S_first_row[0] = 1 + 0j\n",
"\n",
"# Add in ancilla-only measurements:\n",
"for i in range(krylov_dim - 1):\n",
" # Get expectation values from experiment\n",
" expval_real = results.data.evs[0][0][i] # automatic extrapolated evs if ZNE is used\n",
" expval_imag = results.data.evs[1][0][i] # automatic extrapolated evs if ZNE is used\n",
"\n",
" # Get expectation values\n",
" expval = expval_real + 1j * expval_imag\n",
" S_first_row[i + 1] += prefactors[i] * expval\n",
"\n",
"S_first_row_list = S_first_row.tolist() # for saving purposes\n",
"\n",
"\n",
"S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n",
"\n",
"# Distribute entries from first row across matrix:\n",
"for i, j in it.product(range(krylov_dim), repeat=2):\n",
" if i >= j:\n",
" S_circ[j, i] = S_first_row[i - j]\n",
" else:\n",
" S_circ[j, i] = np.conj(S_first_row[j - i])"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "bfe2a427-7ec6-48de-9321-fa6bf2dc7c94",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1.0 & -0.595839719988029 + 0.346522816833994 i & 0.148962122640831 - 0.37923835568426 i & -0.0217605005400922 + 0.0993369468259215 i & 0.167837484202232 + 0.0467433025355653 i\\\\-0.595839719988029 - 0.346522816833994 i & 1.0 & -0.595839719988029 + 0.346522816833994 i & 0.148962122640831 - 0.37923835568426 i & -0.0217605005400922 + 0.0993369468259215 i\\\\0.148962122640831 + 0.37923835568426 i & -0.595839719988029 - 0.346522816833994 i & 1.0 & -0.595839719988029 + 0.346522816833994 i & 0.148962122640831 - 0.37923835568426 i\\\\-0.0217605005400922 - 0.0993369468259215 i & 0.148962122640831 + 0.37923835568426 i & -0.595839719988029 - 0.346522816833994 i & 1.0 & -0.595839719988029 + 0.346522816833994 i\\\\0.167837484202232 - 0.0467433025355653 i & -0.0217605005400922 - 0.0993369468259215 i & 0.148962122640831 + 0.37923835568426 i & -0.595839719988029 - 0.346522816833994 i & 1.0\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 1.0, -0.595839719988029 + 0.346522816833994*I, 0.148962122640831 - 0.37923835568426*I, -0.0217605005400922 + 0.0993369468259215*I, 0.167837484202232 + 0.0467433025355653*I],\n",
"[ -0.595839719988029 - 0.346522816833994*I, 1.0, -0.595839719988029 + 0.346522816833994*I, 0.148962122640831 - 0.37923835568426*I, -0.0217605005400922 + 0.0993369468259215*I],\n",
"[ 0.148962122640831 + 0.37923835568426*I, -0.595839719988029 - 0.346522816833994*I, 1.0, -0.595839719988029 + 0.346522816833994*I, 0.148962122640831 - 0.37923835568426*I],\n",
"[-0.0217605005400922 - 0.0993369468259215*I, 0.148962122640831 + 0.37923835568426*I, -0.595839719988029 - 0.346522816833994*I, 1.0, -0.595839719988029 + 0.346522816833994*I],\n",
"[ 0.167837484202232 - 0.0467433025355653*I, -0.0217605005400922 - 0.0993369468259215*I, 0.148962122640831 + 0.37923835568426*I, -0.595839719988029 - 0.346522816833994*I, 1.0]])"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sympy import Matrix\n",
"\n",
"Matrix(S_circ)"
]
},
{
"cell_type": "markdown",
"id": "2a863b6f-183f-4d7a-bac8-399da3e3578e",
"metadata": {},
"source": [
"And the matrix elements of $\\tilde{H}$"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "c551fdd0-91ef-4531-83d9-ed399b423bb9",
"metadata": {},
"outputs": [],
"source": [
"import itertools\n",
"\n",
"# Assemble S, the overlap matrix of dimension D:\n",
"H_first_row = np.zeros(krylov_dim, dtype=complex)\n",
"H_first_row[0] = H_expval\n",
"\n",
"for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n",
" # Add in ancilla-only measurements:\n",
" for i in range(krylov_dim - 1):\n",
" # Get expectation values from experiment\n",
" expval_real = results.data.evs[2 + 2 * obs_idx][0][\n",
" i\n",
" ] # automatic extrapolated evs if ZNE is used\n",
" expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][\n",
" i\n",
" ] # automatic extrapolated evs if ZNE is used\n",
"\n",
" # Get expectation values\n",
" expval = expval_real + 1j * expval_imag\n",
" H_first_row[i + 1] += prefactors[i] * coeff * expval\n",
"\n",
"H_first_row_list = H_first_row.tolist()\n",
"\n",
"H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n",
"\n",
"# Distribute entries from first row across matrix:\n",
"for i, j in itertools.product(range(krylov_dim), repeat=2):\n",
" if i >= j:\n",
" H_eff_circ[j, i] = H_first_row[i - j]\n",
" else:\n",
" H_eff_circ[j, i] = np.conj(H_first_row[j - i])"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "d1950927-74d6-4012-81f2-3d0b7f476de8",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}10.0 & -3.67474083662792 + 5.79424802432656 i & -2.87080660000195 - 4.50388156185672 i & 3.53539432569443 - 1.04288063424328 i & -0.780365964179053 + 2.94128940749982 i\\\\-3.67474083662792 - 5.79424802432656 i & 10.0 & -3.67474083662792 + 5.79424802432656 i & -2.87080660000195 - 4.50388156185672 i & 3.53539432569443 - 1.04288063424328 i\\\\-2.87080660000195 + 4.50388156185672 i & -3.67474083662792 - 5.79424802432656 i & 10.0 & -3.67474083662792 + 5.79424802432656 i & -2.87080660000195 - 4.50388156185672 i\\\\3.53539432569443 + 1.04288063424328 i & -2.87080660000195 + 4.50388156185672 i & -3.67474083662792 - 5.79424802432656 i & 10.0 & -3.67474083662792 + 5.79424802432656 i\\\\-0.780365964179053 - 2.94128940749982 i & 3.53539432569443 + 1.04288063424328 i & -2.87080660000195 + 4.50388156185672 i & -3.67474083662792 - 5.79424802432656 i & 10.0\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 10.0, -3.67474083662792 + 5.79424802432656*I, -2.87080660000195 - 4.50388156185672*I, 3.53539432569443 - 1.04288063424328*I, -0.780365964179053 + 2.94128940749982*I],\n",
"[ -3.67474083662792 - 5.79424802432656*I, 10.0, -3.67474083662792 + 5.79424802432656*I, -2.87080660000195 - 4.50388156185672*I, 3.53539432569443 - 1.04288063424328*I],\n",
"[ -2.87080660000195 + 4.50388156185672*I, -3.67474083662792 - 5.79424802432656*I, 10.0, -3.67474083662792 + 5.79424802432656*I, -2.87080660000195 - 4.50388156185672*I],\n",
"[ 3.53539432569443 + 1.04288063424328*I, -2.87080660000195 + 4.50388156185672*I, -3.67474083662792 - 5.79424802432656*I, 10.0, -3.67474083662792 + 5.79424802432656*I],\n",
"[-0.780365964179053 - 2.94128940749982*I, 3.53539432569443 + 1.04288063424328*I, -2.87080660000195 + 4.50388156185672*I, -3.67474083662792 - 5.79424802432656*I, 10.0]])"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sympy import Matrix\n",
"\n",
"Matrix(H_eff_circ)"
]
},
{
"cell_type": "markdown",
"id": "7ec60832-f7be-447f-841e-e801241fdeae",
"metadata": {},
"source": [
"Finally, we can solve the generalized eigenvalue problem for $\\tilde{H}$:\n",
"\n",
"$$\\tilde{H} \\vec{c} = c S \\vec{c}$$\n",
"\n",
"and get an estimate of the ground state energy $c_{min}$"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "d6635506-399f-47a3-a837-c5f2558f22b4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The estimated ground state energy is: 10.0\n",
"The estimated ground state energy is: 6.430677870042079\n",
"The estimated ground state energy is: 5.050534767517682\n",
"The estimated ground state energy is: 4.450747729866024\n",
"The estimated ground state energy is: 3.882255130308517\n"
]
}
],
"source": [
"gnd_en_circ_est_list = []\n",
"for d in range(1, krylov_dim + 1):\n",
" # Solve generalized eigenvalue problem\n",
" gnd_en_circ_est = solve_regularized_gen_eig(\n",
" H_eff_circ[:d, :d], S_circ[:d, :d], threshold=1e-1\n",
" )\n",
" gnd_en_circ_est_list.append(gnd_en_circ_est)\n",
" print(\"The estimated ground state energy is: \", gnd_en_circ_est)"
]
},
{
"cell_type": "markdown",
"id": "15865587-85f4-41bc-9bf5-47464cf17298",
"metadata": {},
"source": [
"For a single-particle sector, we can efficiently calculate the ground state of this sector of the Hamiltonian classically"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "3303ae29-c288-417a-9cf2-c82d9770de69",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n_sys_qubits 10\n",
"n_exc 1 , subspace dimension 11\n",
"single particle ground state energy: 2.391547869638771\n"
]
}
],
"source": [
"gs_en = single_particle_gs(H_op, n_qubits)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "2f904ea3-38bc-4841-81ce-cdb69f09a0b7",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"27"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(H_op)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "a5d4a983-1a30-4cea-b695-3e6a67338633",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-2-krylov/extracted-outputs/a5d4a983-1a30-4cea-b695-3e6a67338633-0.avif\" alt=\"Output of the previous code cell\" />"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(\n",
" range(1, krylov_dim + 1),\n",
" gnd_en_circ_est_list,\n",
" color=\"blue\",\n",
" linestyle=\"-.\",\n",
" label=\"KQD estimate\",\n",
")\n",
"plt.plot(\n",
" range(1, krylov_dim + 1),\n",
" [gs_en] * krylov_dim,\n",
" color=\"red\",\n",
" linestyle=\"-\",\n",
" label=\"exact\",\n",
")\n",
"plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))\n",
"plt.legend()\n",
"plt.xlabel(\"Krylov space dimension\")\n",
"plt.ylabel(\"Energy\")\n",
"plt.title(\"Estimating Ground state energy with Krylov Quantum Diagonalization\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "543a2d98-1f72-4df1-ae22-f65d60c0d9c5",
"metadata": {},
"source": [
"## 5. Discussion and extension\n",
"\n",
"To recap, we start with a reference state, then evolve it for different periods of time to generate the unitary Krylov subspace. We project our Hamiltonian onto that subspace. We also estimate the overlaps of the subspace vectors. Finally we solve the lower-dimensional, generalized eigenvalue problem classically.\n",
"\n",
"![A flow-chart overview of QKD: start with a reference state, evolve the state to approximate Krylov vectors, project into the Krylov subspace, diagonalize the projected subspace classically, and determine ground state properties.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig7.avif)\n",
"\n",
"Let’s compare what determines computational costs of using the Krylov technique classically and quantum mechanically. There are not perfect analogs between classical and quantum approaches for all steps. This table collects some scaling of different steps for consideration.\n",
"\n",
"![A table describing scaling of different processes classically and in the quantum approach to Krylov methods. Some quantum steps have no analog. The scalings are the same as those stated in text.](/learning/images/courses/quantum-diagonalization-algorithms/krylov/kqd-fig8.avif)\n",
"\n",
"Recall that Hamiltonians generally have terms that cannot be simultaneously measured (because they do not commute with on another). We sort terms in the Hamiltonian into groups of commuting Pauli operators that can all be measured simultaneously, and we may require many such groups to account for all the terms that do not commute with one another. To build up $\\tilde{H}$ on a quantum computer requires separate measurements for each group of commuting Pauli strings in the Hamiltonian, and each of those requires many shots. We must do this for $r^2$ different matrix elements, corresponding to $r^2$ combinations of different time evolution factors. There are sometimes ways to reduce this, but in this rough treatment, the time required for this scales like $N_\\text{shots}\\times N_\\text{GCP} \\times r^2.$ The elements of $S$ must be estimated, which scales like $O(N_\\text{shots}\\times r^2)$. Finally, solving the generalized eigenvalue problem in the projected space, classically, takes $O(r^3).$\n",
"\n",
"We see that quantum Krylov diagonalization may be useful in cases where the number of commuting Pauli groups in the Hamiltonian is relatively small. These scaling dependencies suggest some applications where the Krylov method can be useful, and others where it likely will not be.\n",
"Some Hamiltonians have high complexity when mapped to qubits, involving many non-commuting Pauli strings that cannot easily be partitioned into a few commuting groups. This is often true of quantum chemistry problems, for example. This complexity presents two primary challenges for near-term quantum computers:\n",
"\n",
"* The estimation of each element of $\\tilde{H}$ becomes computationally expensive due to the large number of terms.\n",
"* The required Trotter circuits become prohibitively deep.\n",
"\n",
"Both of the above points will be less problematic when quantum computers reach fault-tolerance, but they must be considered in the near term. Even systems with “simpler” mappings than those in quantum chemistry may experience the same impediments, if the Hamiltonians have too many non-commuting terms.\n",
"The Krylov method is most useful where the Hamiltonian can be partitioned into relatively few commuting Pauli groups, and where $H$ is easy to implement in trotter circuits. Both of these conditions are satisfied, for example, for many lattice models of interest in physics. KQD is especially useful if very little is known about the ground state. This stems from its inherent convergence guarantees and its applicability in scenarios where alternative methods are untenable due to insufficient ground state knowledge.\n",
"\n",
"While KQD is a powerful tool, the protocol's time-consuming aspects, particularly the estimation of each element of the projected Hamiltonian and the overlap of Krylov states, represent opportunities for improvement. An alternative approach involves leveraging Krylov methods in conjunction with sampling-based methods, which are the subject of the next lesson."
]
},
{
"cell_type": "markdown",
"id": "7fa7609f-9e96-4b0c-b716-ae9242bb40ca",
"metadata": {},
"source": [
"## 6. Appendices\n",
"\n",
"### Appendix I: Krylov subspace from real time-evolutions\n",
"\n",
"The unitary Krylov space is defined as\n",
"\n",
"$$\n",
"\\mathcal{K}_U(H, |\\psi\\rangle) = \\text{span}\\left\\{ |\\psi\\rangle, e^{-iH\\,dt} |\\psi\\rangle, \\dots, e^{-irH\\,dt} |\\psi\\rangle \\right\\}\n",
"$$\n",
"\n",
"for some timestep $dt$ that we will determine later. Temporarily assume $r$ is even: then define $d=r/2$. Notice that when we project the Hamiltonian into the Krylov space above, it is indistinguishable from the Krylov space\n",
"\n",
"$$\n",
"\\mathcal{K}_U(H, |\\psi\\rangle) = \\text{span}\\left\\{ e^{i\\,d\\,H\\,dt}|\\psi\\rangle, e^{i(d-1)H\\,dt} |\\psi\\rangle, \\dots, e^{-i(d-1)H\\,dt} |\\psi\\rangle, e^{-i\\,d\\,H\\,dt} |\\psi\\rangle \\right\\},\n",
"$$\n",
"\n",
"i.e., where all the time-evolutions are shifted backward by $d$ timesteps.\n",
"The reason it is indistinguishable is because the matrix elements\n",
"\n",
"$$\n",
"\\tilde{H}_{j,k} = \\langle\\psi|e^{i\\,j\\,H\\,dt}He^{-i\\,k\\,H\\,dt}|\\psi\\rangle=\\langle\\psi|He^{i(j-k)H\\,dt}|\\psi\\rangle\n",
"$$\n",
"\n",
"are invariant under overall shifts of the evolution time, since the time-evolutions commute with the Hamiltonian. For odd $r$, we can use the analysis for $r-1$.\n",
"\n",
"We want to show that somewhere in this Krylov space, there is guaranteed to be a low-energy state. We do so by way of the following result, which is derived from Theorem 3.1 in [\\[3\\]](#references):\n",
"\n",
"**Claim 1:** there exists a function $f$ such that for energies $E$ in the spectral range of the Hamiltonian (i.e., between the ground state energy and the maximum energy)...\n",
"\n",
"1. $f(E_0)=1$\n",
"2. $|f(E)|\\le2\\left(1 + \\delta\\right)^{-d}$ for all values of $E$ that lie $\\ge\\delta$ away from $E_0$, i.e., it is exponentially suppressed\n",
"3. $f(E)$ is a linear combination of $e^{ijE\\,dt}$ for $j=-d,-d+1,...,d-1,d$\n",
"\n",
"We give a proof below, but that can be safely skipped unless one wants to understand the full, rigorous argument. For now we focus on the implications of the above claim. By property 3 above, we can see that the shifted Krylov space above contains the state $f(H)|\\psi\\rangle$. This is our low-energy state. To see why, write $|\\psi\\rangle$ in the energy eigenbasis:\n",
"\n",
"$$\n",
"|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_k|E_k\\rangle,\n",
"$$\n",
"\n",
"where $|E_k\\rangle$ is the kth energy eigenstate and $\\gamma_k$ is its amplitude in the initial state $|\\psi\\rangle$. Expressed in terms of this, $f(H)|\\psi\\rangle$ is given by\n",
"\n",
"$$\n",
"f(H)|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_kf(E_k)|E_k\\rangle,\n",
"$$\n",
"\n",
"using the fact that we can replace $H$ by $E_k$ when it acts on the eigenstate $|E_k\\rangle$. The energy error of this state is therefore\n",
"\n",
"$$\n",
"\\text{energy error} = \\frac{\\langle\\psi|f(H)(H-E_0)f(H)|\\psi\\rangle}{\\langle\\psi|f(H)^2|\\psi\\rangle}\n",
"$$\n",
"\n",
"$$\n",
"= \\frac{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2}.\n",
"$$\n",
"\n",
"To turn this into an upper bound that is easier to understand, we first separate the sum in the numerator into terms with $E_k-E_0\\le\\delta$ and terms with $E_k-E_0>\\delta$:\n",
"\n",
"$$\n",
"\\text{energy error} = \\frac{\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} + \\frac{\\sum_{E_k> E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2}.\n",
"$$\n",
"\n",
"We can upper bound the first term by $\\delta$,\n",
"\n",
"$$\n",
"\\frac{\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} < \\frac{\\delta\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} \\le \\delta,\n",
"$$\n",
"\n",
"where the first step follows because $E_k-E_0\\le\\delta$ for every $E_k$ in the sum, and the second step follows because the sum in the numerator is a subset of the sum in the denominator. For the second term, first we lower bound the denominator by $|\\gamma_0|^2$, since $f(E_0)^2=1$: adding everything back together, this gives\n",
"\n",
"$$\n",
"\\text{energy error} \\le \\delta + \\frac{1}{|\\gamma_0|^2}\\sum_{E_k>E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0).\n",
"$$\n",
"\n",
"To simplify what is left, notice that for all these $E_k$, by the definition of $f$ we know that $f(E_k)^2 \\le 4\\left(1 + \\delta\\right)^{-2d}$. Additionally upper bounding $E_k-E_0<2\\|H\\|$ and upper bounding $\\sum_{E_k>E_0+\\delta}|\\gamma_k|^2<1$ gives\n",
"\n",
"$$\n",
"\\text{energy error} \\le \\delta + \\frac{8}{|\\gamma_0|^2}\\|H\\|\\left(1 + \\delta\\right)^{-2d}.\n",
"$$\n",
"\n",
"This holds for any $\\delta>0$, so if we set $\\delta$ equal to our goal error, then the error bound above converges towards that exponentially with the Krylov dimension $2d=r$. Also note that if $\\delta<E_1-E_0$ then the $\\delta$ term actually goes away entirely in the above bound.\n",
"\n",
"To complete the argument, we first note that the above is just the energy error of the particular state $f(H)|\\psi\\rangle$, rather than the energy error of the lowest energy state in the Krylov space. However, by the (Rayleigh-Ritz) variational principle, the energy error of the lowest energy state in the Krylov space is upper bounded by the energy error of any state in the Krylov space, so the above is also an upper bound on the energy error of the lowest energy state, i.e., the output of the Krylov quantum diagonalization algorithm.\n",
"\n",
"A similar analysis as the above can be carried out that additionally accounts for noise and the thresholding procedure discussed in the notebook. See [\\[2\\]](#references) and [\\[4\\]](#references) for this analysis.\n",
"\n",
"### Appendix II: proof of Claim 1\n",
"\n",
"The following is mostly derived from [\\[3\\]](#references), Theorem 3.1: let $0 < a < b$ and let $\\Pi^*_d$ be the space of residual polynomials (polynomials whose value at 0 is 1) of degree at most $d$. The solution to\n",
"\n",
"$$\n",
"\\beta(a, b, d) = \\min_{p \\in \\Pi^*_d} \\max_{x \\in [a, b]} |p(x)| \\quad\n",
"$$\n",
"\n",
"is\n",
"\n",
"$$\n",
"p^*(x) = \\frac{T_d\\left(\\frac{b + a - 2x}{b - a}\\right)}{T_d\\left(\\frac{b + a}{b - a}\\right)}, \\quad\n",
"$$\n",
"\n",
"and the corresponding minimal value is\n",
"\n",
"$$\n",
"\\beta(a, b, d) = T_d^{-1}\\left(\\frac{b + a}{b - a}\\right).\n",
"$$\n",
"\n",
"We want to convert this into a function that can be expressed naturally in terms of complex exponentials, because those are the real time-evolutions that generate the quantum Krylov space.\n",
"To do so, it is convenient to introduce the following transformation of energies within the spectral range of the Hamiltonian to numbers in the range $[0,1]$: define\n",
"\n",
"$$\n",
"g(E) = \\frac{1-\\cos\\big((E-E_0)dt\\big)}{2},\n",
"$$\n",
"\n",
"where $dt$ is a timestep such that $-\\pi < E_0dt < E_\\text{max}dt < \\pi$.\n",
"Notice that $g(E_0)=0$ and $g(E)$ grows as $E$ moves away from $E_0$.\n",
"\n",
"Now using the polynomial $p^*(x)$ with the parameters a, b, d set to $a = g(E_0 + \\delta)$, $b = 1$, and d = int(r/2), we define the function:\n",
"\n",
"$$\n",
"f(E) = p^* \\left( g(E) \\right) = \\frac{T_d\\left(1 + 2\\frac{\\cos\\big((E-E_0)dt\\big) - \\cos\\big(\\delta\\,dt\\big)}{1 +\\cos\\big(\\delta\\,dt\\big)}\\right)}{T_d\\left(1 + 2\\frac{1-\\cos\\big(\\delta\\,dt\\big)}{1 + \\cos\\big(\\delta\\,dt\\big)}\\right)}\n",
"$$\n",
"\n",
"where $E_0$ is the ground state energy. We can see by inserting $\\cos(x)=\\frac{e^{ix}+e^{-ix}}{2}$ that $f(E)$ is a trigonometric polynomial of degree $d$, i.e., a linear combination of $e^{ijE\\,dt}$ for $j=-d,-d+1,...,d-1,d$. Furthermore, from the definition of $p^*(x)$ above we have that $f(E_0)=p(0)=1$ and for any $E$ in the spectral range such that $\\vert E-E_0 \\vert > \\delta$ we have\n",
"\n",
"$$\n",
"|f(E)| \\le \\beta(a, b, d) = T_d^{-1}\\left(1 + 2\\frac{1-\\cos\\big(\\delta\\,dt\\big)}{1 + \\cos\\big(\\delta\\,dt\\big)}\\right)\n",
"$$\n",
"\n",
"$$\n",
"\\leq 2\\left(1 + \\delta\\right)^{-d} = 2\\left(1 + \\delta\\right)^{-\\lfloor k/2\\rfloor}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "97d292db-83bb-431c-97e3-91212e168ab8",
"metadata": {},
"source": [
"## References:\n",
"\n",
"[1] https://arxiv.org/abs/2407.14431\n",
"\n",
"[2] https://arxiv.org/abs/1811.09025\n",
"\n",
"[3] https://people.math.ethz.ch/~mhg/pub/biksm.pdf\n",
"\n",
"[4] https://academic.oup.com/book/36426\n",
"\n",
"[5] https://en.wikipedia.org/wiki/Krylov_subspace\n",
"\n",
"[6] Krylov Subspace Methods: Principles and Analysis, Jörg Liesen, Zdenek Strakos https://academic.oup.com/book/36426?login=false\n",
"\n",
"[7] Iterative Methods for Sparse Linear Systems\" by Yousef Saad\n",
"\n",
"[8] \"MINRES-QLP: A Krylov Subspace Method for Indefinite or Singular Symmetric Systems\" by Sou-Cheng Choi, Christopher Paige, and\n",
"Michael Saunders (https://epubs.siam.org/doi/10.1137/100787921)\n",
"\n",
"[9] Ethan N. Epperly, Lin Lin, and Yuji Nakatsukasa. \"A theory of quantum subspace diagonalization\". SIAM Journal on Matrix Analysis and Applications 43, 1263–1290 (2022).\n",
"\n",
"[10] https://link.aps.org/doi/10.1103/PRXQuantum.4.030319"
]
}
],
"metadata": {
"description": "Quantum Krylov Diagonalization (QKD) is described, starting with classical Krylov methods. Convergence and resource intensiveness are discussed.",
"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"
},
"title": "Quantum Krylov Diagonalization"
},
"nbformat": 4,
"nbformat_minor": 2
}
Reference in New Issue View Git Blame Copy Permalink