크릴로프 양자 대각화

이번 크릴로프 양자 대각선화(KQD) 강의에서는 다음과 같은 질문에 답해드리겠습니다:

• 일반적으로 크릴로프 방식이란 무엇인가요?
• 크릴로프 방식은 왜 작동하며 어떤 조건에서 작동하나요?
• 양자 컴퓨팅은 어떤 역할을 하나요?

계산의 양자 부분은 주로 참고 문헌 [1] 의 연구를 기반으로 합니다.

아래 동영상은 클래식 컴퓨팅의 크릴로프 메서드에 대한 개요와 사용 동기, 그리고 해당 워크스트림에서 양자 컴퓨팅이 어떤 역할을 할 수 있는지 설명합니다. 이어지는 텍스트는 더 자세한 내용을 제공하며 고전적인 방법과 양자 컴퓨터를 사용하여 크릴로프 방법을 구현합니다.

1. 크릴로프 방법 소개

크릴로프 부분 공간 방법은 크릴로프 부분 공간 이라고 불리는 것을 중심으로 구축된 여러 가지 방법을 지칭할 수 있습니다. 이에 대한 전체 검토는 이 강의의 범위를 벗어나지만 참고자료 [2-4 ]를 참고하면 훨씬 더 많은 배경 지식을 얻을 수 있습니다. 여기서는 크릴로프 부분공간이 무엇인지, 고유값 문제를 푸는 데 유용한 방법과 이유는 무엇인지, 마지막으로 양자 컴퓨터에서 어떻게 구현할 수 있는지에 대해 집중적으로 살펴보겠습니다.

정의: 대칭적이고 양의 반정적 $N\times N$ 행렬 $A$ 이 주어지면, $r$ 차수의 크릴로프 공간 $\mathcal{K}^r$ 은 행렬 $A$ 의 상위 거듭제곱을 $r-1\leq N$ 까지 곱하여 얻은 벡터와 기준 벡터 $\vert v \rangle$ 가 걸쳐 있는 공간입니다.

위의 벡터는 크릴로프 아공간이라고 부르는 공간에 걸쳐 있지만, 직교할 것이라고 생각할 이유는 없습니다. 종종 그램-슈미트 직교 정규 화와 유사한 반복적인 직교 정규화 프로세스를 사용합니다. 여기서는 각각의 새 벡터가 생성될 때 다른 벡터와 직교하도록 만들어지므로 프로세스가 약간 다릅니다. 이러한 맥락에서 이를 아놀디 반복이라고 합니다. 초기 벡터 $|v\rangle$ 로 시작하여 다음 벡터 $A|v\rangle$ 를 생성한 다음 $|v\rangle$ 에서 투영을 빼서 이 두 번째 벡터가 첫 번째 벡터에 직교하도록 합니다. 즉

이제 $|v_0\rangle \perp |v_1\rangle,$ 이후

다음 벡터에 대해서도 동일한 작업을 수행하여 앞의 두 벡터와 직교하도록 합니다:

모든 $r$ 벡터에 대해 이 과정을 반복하면 크릴로프 공간에 대한 완전한 정형 기반이 생깁니다. $m$ 직교 벡터는 반드시 전체 공간에 걸쳐야 하므로 여기서 직교화 프로세스를 수행하면 $r>m$ 이 0을 반환한다는 점에 유의하세요. 이후의 모든 벡터는 해당 벡터의 배수가 되므로 어떤 벡터가 $A$ 의 고유 벡터인 경우에도 이 프로세스는 0을 산출합니다.

1.1 간단한 예시: 수작업으로 크릴로프

그 과정을 볼 수 있도록 아주 작은 행렬에서 크릴로프 부분공간 생성의 한 세대를 살펴봅시다. 관심 있는 초기 매트릭스 $A$ 부터 시작합니다:

이 작은 예제에서는 손으로도 쉽게 고유 벡터와 고유값을 구할 수 있습니다. 여기에 수치화된 솔루션을 보여드립니다.

# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use
# linalg.eig in this lesson.

import numpy as np

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

Output:

The eigenvalues are  [2.58578644 4.         5.41421356]
The eigenvectors are  [[ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]
 [ 7.07106781e-01  1.37464400e-16 -7.07106781e-01]
 [ 5.00000000e-01  7.07106781e-01  5.00000000e-01]]

나중에 비교할 수 있도록 여기에 기록합니다:

Krylov 하위 공간인 $r$ 의 차원을 늘릴 때 이 프로세스가 어떻게 작동하는지(또는 실패하는지) 연구해보고자 합니다. 이를 위해 이 프로세스를 적용해 보겠습니다:

• 무작위로 선택한 벡터 $|v\rangle$ (위와 같이 이미 정규화된 경우 $|v_0\rangle$ )로 시작하여 전체 벡터 공간의 부분공간을 생성합니다.
• 전체 행렬 $A$ 을 해당 부분 공간에 투영하고 투영된 행렬의 고유값을 구합니다 $\tilde{A}$.
• 그램-슈미트 직교화와 유사한 프로세스를 사용하여 더 많은 벡터를 생성하고 직교가 되도록 하여 하위 공간의 크기를 늘립니다.
• $A$ 을 더 큰 부분 공간에 투영하고 결과 행렬의 고유값( $\tilde{A}$ )을 구합니다.
• 고유값이 수렴할 때까지(또는 이 장난감의 경우 원래 행렬의 전체 벡터 공간에 걸친 벡터를 생성할 때까지 $A$ ) 이 작업을 반복합니다.

일반적인 크릴로프 방법의 구현에서는 모든 크릴로프 하위 공간에 투영된 행렬의 고유값 문제를 풀지 않아도 됩니다. 원하는 차원의 하위 공간을 구성하고 해당 하위 공간에 행렬을 투영한 다음 투영된 행렬을 대각선화할 수 있습니다. 각 하위 공간 차원에서의 투영 및 대각선화는 수렴을 확인하기 위해서만 수행됩니다.

$r=1$ 치수:

임의의 벡터를 선택합니다

아직 정규화되지 않은 경우 정규화합니다.

이제 행렬 $A$ 을 이 하나의 벡터의 부분 공간에 투영합니다:

이것은 행렬에 벡터가 하나만 포함된 경우 행렬을 크릴로프 부분공간에 투영한 것입니다( $|v_0\rangle$ ). 이 행렬의 고유값은 사소하게도 4입니다. 이를 $A$ 의 고유값(이 경우 단 하나)에 대한 0차 추정치라고 생각할 수 있습니다. 비록 잘못된 추정치이기는 하지만 정확한 크기 순서입니다.

$r=2$ 치수:

이제 이전 벡터에 $A$ 연산을 통해 하위 공간에서 다음 벡터를 생성합니다:

이제 직교성을 보장하기 위해 이 벡터를 이전 벡터에 투영한 값을 뺍니다.

아직 정규화되지 않은 경우 정규화합니다. 이 경우 벡터는 이미 정규화되었으므로

이제 행렬 A를 이 두 벡터의 부분 공간에 투영합니다:

이 행렬의 고유값을 결정하는 문제는 여전히 남아 있습니다. 하지만 이 매트릭스는 전체 매트릭스보다 약간 작습니다. 매우 큰 행렬을 포함하는 문제에서는 이 작은 하위 공간으로 작업하는 것이 매우 유리할 수 있습니다.

여전히 좋은 추정치는 아니지만 0차 추정치보다는 낫습니다. 프로세스를 명확히 하기 위해 이 작업을 한 번 더 반복할 예정입니다. 그러나 이렇게 하면 다음 반복에서 3x3 행렬을 대각선으로 만들게 되므로 시간이나 계산 능력을 절약하지 못하므로 이 방법의 요점이 약해집니다.

$r=3$ 치수:

이제 이전 벡터에서 A를 연산하여 하위 공간에 다음 벡터를 생성합니다:

이제 직교성을 보장하기 위해 이 벡터를 이전 두 벡터에 투영한 값을 뺍니다.

아직 정규화되지 않은 경우 정규화합니다. 이 경우 벡터는 이미 정규화되었으므로

이제 행렬 $A$ 을 이 벡터의 하위 공간에 투영합니다:

이제 고유값을 결정합니다:

이 고유값은 원래 행렬 $A$ 의 고유값과 정확히 일치합니다. 원래 행렬 $A$ 의 전체 벡터 공간에 걸쳐 있도록 크릴로프 부분공간을 확장했기 때문에 반드시 그래야 합니다.

이 예시에서는 크릴로프 방식이 직접 대각선화보다 특별히 쉬워 보이지 않을 수 있습니다. 실제로 뒷부분에서 살펴보겠지만 크릴로프 방법은 특정 행렬 차원 이상에서만 유리하며, 이는 매우 큰 행렬의 고유값/고유벡터 문제를 해결하는 데 도움을 주기 위한 것입니다.

여기서는 '수작업'으로 작업한 예제만 보여드리지만, 아래 섹션 2에서는 계산 예제를 보여드리겠습니다.

용어 설명

흔히 오해하는 것은 주어진 문제에 대해 하나의 크릴로프 하위 공간만 존재한다는 것입니다. 물론 행렬을 적용할 수 있는 초기 벡터가 많기 때문에 가능한 크릴로프 부분공간도 많습니다. 여기서는 특정 예제에 대해 이미 정의된 특정 크릴로프 하위 공간을 지칭할 때만 " 크릴로프 하위 공간"이라는 표현을 사용합니다. 일반적인 문제 해결 접근 방식에 대해서는 " 크릴로프 하위 공간"을 참조하겠습니다. 마지막으로 '크릴로프 스페이스 '를 참조하는 것이 유효하다는 점을 명확히 합니다. 초기 공간에서 행렬을 부분 공간으로 투영하는 맥락에서 사용되기 때문에 흔히 '크릴로프 부분 공간 '이라고 부릅니다. 이러한 맥락에 따라 여기서는 주로 하위 공간으로 지칭하겠습니다.

이해도 점검

크릴로프 아공간 $r$ 의 차원을 관심 행렬의 차원 $N$ 이상으로 확장하는 것이 (가) 유용하지 않은 이유와 (나) 가능한 이유를 설명하십시오.

• 답변

  (a) 벡터를 생성하는 과정에서 정교범준화를 수행하고 있으므로, 이러한 벡터들로 구성된 집합 $N$ 은 완전 기저를 이룬다. 즉, 이 벡터들의 선형 조합을 사용하여 해당 공간 내의 모든 벡터를 생성할 수 있다.

  (b) 직교화 과정은 새로운 벡터의 모든 기존 벡터에 대한 투영을 빼는 것으로 이루어진다. 만약 이전의 모든 벡터들이 전체 벡터 공간을 생성한다면, 전체 부분공간으로의 투영을 뺀 결과는 항상 영벡터가 된다.

동료 연구원이 작은 예제 행렬에 크릴로프 방법을 적용해 보여주는 상황을 가정해 보자. 매트릭스 $A$ 와 초기 벡터 $|\psi\rangle$ 의 선택에 문제가 있는 건가요?

및

• 답변

  동료가 실수로 초기 벡터로 고유 벡터를 선택했습니다. 초기 벡터에 행렬을 적용하면 단순히 고유값으로 스케일링된 동일한 벡터를 반환합니다. 이렇게 하면 차원이 증가하는 하위 공간이 생성되지 않습니다. 동료에게 다른 초기 벡터를 선택하도록 조언하고 고유 벡터가 아닌지 확인합니다.

주어진 행렬에 크릴로프(Krylov) 방법을 적용하되, 적절한 새로운 초기 벡터를 선택하십시오. 크릴로프 부분공간에 대한 0차 및 1차 근사값을 적어 보십시오.

• 답변

  초기 벡터의 선택에 따라 다양한 답이 나올 수 있습니다. 저희가 선택하겠습니다:

  $$|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}.$$

  $|v_1\rangle$ 를 얻으려면 $A$ 을 $|v_0\rangle$ 에 한 번 적용한 다음 $|v_1\rangle$ 을 다음과 같이 직교시킵니다 $|v_0\rangle.$

  $$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}$$$$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}$$

  0차에서 크릴로프 부분공간으로의 투영은 다음과 같습니다

  $$\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}$$

  1 차에서 이 크릴로프 부분공간에 투영하는 것은 다음과 같습니다

  $$\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}$$

  이 작업은 손으로 할 수도 있지만 넘피를 사용하면 가장 쉽게 할 수 있습니다:

  import numpy as np
  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)]]
  )
  A = np.array([[1, 1, 0],
                [1, 1, 1],
                [0, 1, 1]])
  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)]])
  proj = vstar@A@v
  print(proj)
  eigenvalues, eigenvectors = np.linalg.eig(proj)
  print("The eigenvalues are ", eigenvalues)
  print("The eigenvectors are ", eigenvectors)

  outputs:

  [[ 2.33333333  0.47140452]
   [ 0.47140452 -0.33333333]]
  The eigenvalues are  [ 2.41421356 -0.41421356]
  The eigenvectors are  [[ 0.98559856 -0.16910198]
   [ 0.16910198  0.98559856]]

  최소 고유값 추정치는 -0.414 입니다.

1.2 크릴로프 방법의 유형

"크릴로프 아공간 방법"은 대규모 선형 시스템과 고유값 문제를 푸는 데 사용되는 여러 반복 기법 중 하나를 지칭할 수 있습니다. 이들 모두의 공통점은 크릴로프 부분 공간에서 대략적인 솔루션을 구축한다는 것입니다

$\mathcal{K}^r(A,|v\rangle ) = \text{span}\{|v\rangle, A|v\rangle, A^2|v\rangle, ..., A^{r-1}|v\rangle\},$

여기서 $|v\rangle$ 은 초기 추측입니다(참고 [5] 참조). 이 하위 공간에서 최적의 근사치를 선택하는 방식이 다르며 수렴 속도, 메모리 사용량, 전체 계산 비용과 같은 요소의 균형을 맞추는 방식이 다릅니다. 이 강의의 초점은 크릴로프 아공간 방법의 맥락에서 양자 컴퓨팅을 활용하는 것입니다. 이러한 방법에 대한 자세한 논의는 이 강의의 범위를 벗어납니다. 아래의 간략한 정의는 맥락을 위한 것이며, 이러한 방법을 더 자세히 조사하기 위한 몇 가지 참고 자료가 포함되어 있습니다.

공액 그라디언트(CG) 방법 : 이 방법은 대칭, 양의 정적 선형 시스템을 푸는 데 사용됩니다 [6]. 각 반복에서 오차의 A-노름을 최소화하므로 이산 타원형 PDE로 인해 발생하는 시스템에 특히 효과적입니다 [7]. 다음 섹션에서는 이 접근법을 사용해 왜 크릴로프 아공간이 선형 시스템의 개선된 솔루션을 탐색하는 데 효과적인 아공간이 될 수 있는지에 대해 설명하겠습니다.

일반화된 최소 잔차(GMRES) 방법 : 일반적인 비대칭 선형 시스템을 풀기 위해 설계되었습니다. 각 반복에서 크릴로프 공간에 대한 잔여 노멀을 최소화하여 강력하지만 대규모 시스템에서는 잠재적으로 메모리 집약적일 수 있습니다 [7].

최소 잔차(MINRES) 방법 : 이 방법은 대칭 무기한 선형 시스템을 푸는 데 사용됩니다. GMRES와 유사하지만 계산 비용을 줄이기 위해 행렬 대칭성을 활용합니다 [8].

주목할 만한 다른 접근법으로는 고유값 문제에 대한 아놀디의 방법과 밀접한 관련이 있는 완전 직교화 방법(FOM), 이중 공액 그라데이션( BiCG ) 방법, 유도 차원 축소(IDR) 방법 등이 있습니다.

1.3 크릴로프 부분공간 방법이 작동하는 이유

여기에서는 크릴로프 아공간 방법이 가파른 하강의 렌즈를 통해 행렬 고유 벡터 근사치를 반복적으로 개선하여 행렬 고유값을 근사화하는 효율적인 방법이라는 것을 설명합니다. 접지 상태에 대한 초기 추측이 주어졌을 때 가장 빠르게 수렴하는 초기 추측을 연속적으로 수정하는 공간은 크릴로프 부분공간이라고 주장할 것입니다. 컨버전스 동작에 대한 엄격한 증명이 부족합니다.

관심 행렬( $A$ )이 대칭이고 양정확하다고 가정합니다. 따라서 위의 CG 메서드와 가장 관련이 있습니다. 여기서는 희소성에 대해 가정하지 않으며, $A$ 이 반드시 해밀턴주의자이어야 한다고 주장하지도 않습니다(해밀턴주의자라면 그래야 합니다).

일반적으로 다음과 같은 형태의 문제를 해결하고자 합니다

$|b\rangle=c|x\rangle$ 여기서 $c$ 은 고유값 문제에서와 같이 어떤 상수라고 생각할 수 있습니다. 하지만 현재로서는 문제 진술이 더 일반적입니다.

대략적인 솔루션인 $|x_0\rangle$ 벡터로 시작합니다. 섹션 1.1 의 $|x_0\rangle$ 과 $|v_0\rangle$ 사이에는 유사점이 있지만 여기서는 이를 활용하지 않습니다. $|x_0\rangle$ 에 오류가 있습니다 $|e_0\rangle:$

또한 잔여 $R_0:$

여기서는 자본 $R$ 을 사용하여 잔여를 크릴로프 하위 공간 $r$ 의 차원과 구분합니다.

이제 다음과 같은 형식의 수정 단계를 만들고 싶습니다

를 사용하면 근사치가 개선될 것으로 기대합니다. 여기 $|p_0\rangle$ 는 아직 결정되지 않은 벡터입니다. $|e_1\rangle$ 을 수정 후 오류로 설정합니다. 그러면

진정한

우리는 매트릭스에 의해 변환될 때 오류가 어떻게 작동하는지 궁금합니다. 따라서 오류의 $A$ -규범을 계산해 보겠습니다. 즉

여기서 $A$ 의 대칭을 사용했고 $A |e_0\rangle = |R_0\rangle.$ 여기서 $d$ 는 $|p_0\rangle$ 와 독립적인 상수입니다. 섹션 1.2 에서 언급했듯이 오류의 $A$ -규범은 우리가 최소화하기 위해 선택할 수 있는 유일한 양은 아니지만, 좋은 양입니다. 선택한 보정 벡터에 따라 이 수량이 어떻게 달라지는지 확인하고 싶습니다 $|p_0\rangle.$ 따라서 $f$ 함수를 정의하여 다음과 같이 설정합니다

$f$ 는 $A$ -norm에서 측정된 수정값 $|p_0\rangle$ 의 함수인 $|e_1\rangle$ 에 불과합니다. 따라서 $f(|p_0\rangle)$ 이 가능한 한 작도록 $|p_0\rangle$ 을 선택하려고 합니다. 이를 위해 $f$ 의 그라데이션을 계산합니다. $A$ 의 대칭을 사용하여 다음과 같이 계산합니다

기울기는 가장 가파르게 상승하는 방향을 가리키며, 그 반대는 함수가 가장 많이 감소하는 방향, 즉 가장 가파르게 하강하는 방향을 나타냅니다. 초기 추측으로는 $|x_0\rangle$, 여기서 $|p_0\rangle=0$ 이 있습니다 $\nabla f(0) = -2|R_0\rangle.$ 따라서 $f$ 함수는 잔류 방향에서 가장 많이 감소합니다 $|R_0\rangle.$ 따라서 초기 선택은 일부 스칼라 $|p_0\rangle=\alpha_0 |R_0\rangle$ 에 대해 벡터 $\alpha_0$ 를 추가하는 것이 가장 유리할 것입니다.

다음 단계에서는 다시 벡터 $|p_1\rangle$ 를 선택하고 현재 근사값에 해당 값을 더합니다. 이전과 동일한 인수를 사용하여 일부 스칼라 $\alpha_1$ 에 대해 $|p_1\rangle = \alpha_1 |R_1\rangle$ 를 선택합니다. 이런 식으로 계속 진행하면 벡터의 $k^\text{th}$ 반복은 다음과 같습니다

마찬가지로 $|R_0\rangle$, $|R_1\rangle$ 등을 순서대로 추가하여 개선된 견적을 선택할 수 있는 공간을 구축하고자 합니다. $k^\text{th}$ 예상 벡터는 다음과 같습니다

이제 다음과 같은 관계를 사용하여

우리는

즉, 정답 $|x\rangle$ 에 가장 효율적으로 근사하는 공간은 $|R_0\rangle.$ 에서 행렬 $A$ 을 연속적으로 연산하여 구축한 공간입니다. 크릴로프 부분공간은 가장 가파른 하강 방향의 벡터가 연속적으로 걸쳐 있는 공간입니다.

마지막으로, 저희는 이 접근 방식의 확장성에 대해 수치적으로 주장하지 않았으며, 희소 행렬의 비교 이점에 대해서도 논의하지 않았음을 다시 한 번 강조합니다. 이는 크릴로프 서브스페이스 메서드 사용에 대한 동기를 부여하고 직관적인 의미를 더하기 위한 것입니다. 이제 이러한 메서드의 동작을 수치로 살펴보겠습니다.

이해도 점검

위의 워크플로에서는 오류의 $A$ -규범을 최소화하는 방법을 제안했습니다. 접지 상태와 고유값을 구할 때 최소화할 수 있는 다른 수량은 무엇일까요?

• 답변

  오차의 $A$ -규범 대신 잔차 벡터를 사용하는 것을 상상할 수 있습니다. 오차 벡터 자체를 고려하는 것이 유용한 경우가 있을 수 있습니다.

2. 고전적 계산에서의 크릴로프 방법

이 섹션에서는 아놀디 반복을 계산적으로 구현하여 고유값 문제를 푸는 데 크릴로프 부분공간을 활용할 수 있도록 합니다. 먼저 이를 소규모 예제에 적용한 다음 관심 행렬의 크기가 커짐에 따라 계산 시간이 어떻게 확장되는지 살펴보겠습니다. 여기서 핵심 아이디어는 크릴로프 공간에 걸친 벡터의 생성이 필요한 총 컴퓨팅 시간에 큰 영향을 미친다는 것입니다. 필요한 메모리는 특정 크릴로프 메서드에 따라 다릅니다. 하지만 메모리 제약으로 인해 기존 크릴로프 메서드의 사용이 제한될 수 있습니다.

2.1 간단한 소규모 예시

크릴로프 부분공간을 만드는 과정에서 부분공간에 있는 벡터를 직교 정규화해야 합니다. 하위 공간 vknown (정규화되지 않은 것으로 가정)에서 기존 벡터와 하위 공간 vnext 에 추가할 후보 벡터를 가져와 vnext 를 vknown 에 직교하고 정규화하는 함수를 정의해 보겠습니다. 크릴로프 부분 공간에서 확립된 모든 벡터에 대해 이 과정을 거쳐 완전한 정규 집합을 보장하는 함수를 추가로 정의해 보겠습니다.

# vknown is some established vector in our subspace. vnext is one we wish to add,
# which must be orthogonal to vknown.


def orthog_pair(vknown, vnext):
    vknown = vknown / np.sqrt(vknown.T @ vknown)
    diffvec = vknown.T @ vnext * vknown
    vnext = vnext - diffvec
    return vnext


# v is the candidate vector to be added to our subspace. s is the existing subspace.


def orthoset(v, s):
    v = v / np.sqrt(v.T @ v)
    temp = v
    for i in range(len(s)):
        temp = orthog_pair(s[i], temp)
    v = temp / np.sqrt(temp.T @ temp)
    return v

이제 크릴로프 벡터의 공간이 원래 행렬의 전체 공간에 걸쳐 있을 때까지 반복적으로 점점 더 큰 크릴로프 부분공간을 구축하는 함수를 정의해 보겠습니다. 이렇게 하면 크릴로프 부분공간 방법을 사용하여 얻은 고유값이 크릴로프 부분공간 차원에 따라 정확한 값과 얼마나 잘 일치하는지 확인할 수 있습니다. 중요한 것은 krylov_full_build 함수가 크릴로프 벡터, 투영된 해밀토니안, 고유값, 필요한 시간을 반환한다는 점입니다.

# Necessary imports and definitions to track time in microseconds
import time


def time_mus():
    return int(time.time() * 1000000)


# This function constructs a Krylov subspace that spans the whole space of the original matrix.
#     Input:
#       v0          : initial vector
#       matrix      : original matrix to be diagonalized
#     Output:
#       ks          : Krylov vectors
#       Hs          : projected Hamiltonians
#       eigs        : eigenvalues
#       k_tot_times : time required for the operation


def krylov_full_build(v0, matrix):
    t0 = time_mus()
    b = v0 / np.sqrt(v0 @ v0.T)
    A = matrix
    ks = []
    ks.append(b)
    Hs = []
    eigs = []
    Hs.append(b.T @ A @ b)
    eigs.append(np.array([b.T @ A @ b]))
    k_tot_times = []

    for j in range(len(A) - 1):
        vec = A @ ks[j].T
        ortho = orthoset(vec, ks)
        ks.append(ortho)
        ksarray = np.array(ks)
        Hs.append(ksarray @ A @ ksarray.T)
        eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)
        k_tot_times.append(time_mus() - t0)

    # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues,
    # and the total time required.
    return (ks, Hs, eigs, k_tot_times)

아직은 매우 작지만 수작업으로 할 수 있는 것보다 큰 매트릭스에서 이를 테스트해 보겠습니다.

# Define our small test matrix
test_matrix = np.array(
    [
        [4, -1, 0, 1, 0],
        [-1, 4, -1, 2, 1],
        [0, -1, 4, 3, 3],
        [1, 2, 3, 4, 0],
        [0, 1, 3, 0, 4],
    ]
)

# Give the test matrix and an initial guess as arguments in the function defined above.
# Calculate outputs.
test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(
    np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix
)

마지막 단계(크릴로프 공간이 원래 행렬의 전체 벡터 공간인 경우)에서 크릴로프 방법의 고유값이 정확한 수치 대각선화의 고유값과 정확히 일치하는지 확인하여 함수를 확인할 수 있습니다:

print(np.linalg.eig(test_matrix).eigenvalues)
print(test_eigs[len(test_matrix) - 1])

Output:

[-1.36956923  8.43756009  2.9040308   5.34436028  4.68361806]
[-1.36956923  8.43756009  2.9040308   4.68361806  5.34436028]

성공적이었습니다. 물론 실제로 중요한 것은 크릴로프 아공간 차원에 따라 근사치가 얼마나 좋은지입니다. 우리는 종종 기저 상태 및 기타 최소 고유값을 찾는 데 관심이 있기 때문에(그리고 아래에 설명된 다른 대수적 이유 때문에) 크릴로프 아공간 차원의 함수로서 최저 고유값의 추정치를 살펴 보겠습니다. 즉

def errors(matrix, krylov_eigs):
    targ_min = min(np.linalg.eig(matrix).eigenvalues)
    err = []
    for i in range(len(matrix)):
        err.append(min(krylov_eigs[i]) - targ_min)
    return err

import matplotlib.pyplot as plt

krylov_error = errors(test_matrix, test_eigs)

plt.plot(krylov_error)
plt.axhline(y=0, color="red", linestyle="--")  # Add dashed red line at y=0
plt.xlabel("Order of Krylov subspace")  # Add x-axis label
plt.ylabel("Error in minimum eigenvalue")  # Add y-axis label
plt.show()

Output:

크릴로프 부분공간이 $\mathcal{K}^2,$ 까지 커지면 최소 고유값이 상당히 정확하게 도달하고 다음과 같이 완벽해짐을 알 수 있습니다 $\mathcal{K}^3.$

2.2 시간 스케일링과 행렬 차원

다음과 같은 방식으로 크릴로프 방식이 정확한 수치 아이겐솔버보다 유리할 수 있음을 확신해 보겠습니다:

• 무작위 행렬 구성(희소 행렬이 아니며, KQD에 이상적인 애플리케이션은 아님)
• NumPy 을 직접 사용하는 방법과 크릴로프 부분공간을 사용하는 두 가지 방법을 사용하여 고유값을 결정합니다.
• 크릴로프 추정치를 수락하기 전에 고유값이 얼마나 정확해야 하는지에 대한 컷오프를 선택합니다.
• 이 두 가지 방법으로 해결하는 데 필요한 벽 시간을 비교하세요.

주의 사항: 아래에서 자세히 설명하겠지만, 크릴로프 양자 대각선화는 행렬 표현이 희소하거나 소수의 통근 파울리 연산자 그룹을 사용하여 작성할 수 있는 연산자에 가장 잘 적용됩니다. 여기서 사용하는 무작위 행렬은 이러한 설명에 맞지 않습니다. 이는 기존 크릴로프 메서드가 유용할 수 있는 규모를 조사하는 데에만 유용합니다. 둘째, 크릴로프 방법을 사용할 때 다양한 크기의 크릴로프 부분공간을 사용하여 고유값을 계산합니다. 접지 상태 고유값에 대해 필요한 정확도를 달성하는 최소 차원 크릴로프 부분공간에 필요한 시간을 보고합니다. 다시 말하지만, 이것은 정확한 해를 사용하여 필요한 차원을 평가하기 때문에 정확한 아이겐솔버가 풀기 어려운 문제를 푸는 것과는 약간 다릅니다.

먼저 무작위 행렬 세트를 생성합니다.

import numpy as np

# Set the random seed
np.random.seed(42)

# how many random matrices will we make
num_matrix = 200

matrices = []
for m in range(1, num_matrix):
    matrices.append(np.random.rand(m, m))

이제 numpy를 사용하여 각 행렬을 직접 대각선화합니다. 나중에 비교하기 위해 대각선화에 필요한 시간을 계산합니다.

matrix_numpy_times = []
matrix_numpy_eigs = []
for mm in range(num_matrix - 1):
    t0 = time_mus()
    matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))
    matrix_numpy_times.append(time_mus() - t0)

plt.plot(matrix_numpy_times)
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

위 이미지에서 125 차원 주변의 비정상적으로 높은 시간은 행렬의 무작위 특성 또는 사용된 클래식 프로세서에서의 구현으로 인한 것일 수 있지만 재현할 수는 없습니다. 코드를 다시 실행하면 다른 비정상 피크가 있는 다른 프로필이 생성됩니다.

이제 각 행렬에 대해 크릴로프 부분공간을 구축하고 단계적으로 고유값을 계산합니다. 각 단계에서 지정된 절대 오차 범위 내에서 가장 낮은 고유값이 얻어졌는지 확인합니다. 지정된 오차 내에서 고유값을 먼저 제공하는 하위 공간은 계산 시간을 기록할 하위 공간입니다. 이 셀을 실행하는 데는 프로세서 속도에 따라 몇 분 정도 걸릴 수 있습니다. 평가를 건너뛰거나 대각선 행렬의 최대 차원을 자유롭게 줄이세요. 미리 계산된 결과를 보는 것으로 충분합니다.

# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size
# at which that error is achieved.
abserr = 0.05
accept_subspace_size = []

# Lists to store total time spent on the Krylov method, and the subset of that time spent on
# diagonalizing the projected matrix.
matrix_krylov_tot_times = []
matrix_krylov_dim = []

# Step through all our random matrices
for mm in range(0, num_matrix - 1):
    test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(
        np.ones(len(matrices[mm])), matrices[mm]
    )
    # We have not yet found a Krylov subspace that produces our minimum eigenvalue to
    # within the required error.
    found = 0
    for j in range(0, len(matrices[mm]) - 1):
        # If we still haven't found the desired subspace...
        if found == 0:
            # ...but if this one satisfies the requirement, then record everything
            if (
                abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])
                < abserr
            ):
                accept_subspace_size.append(j)
                matrix_krylov_tot_times.append(test_k_tot_times[j])
                matrix_krylov_dim.append(mm)
                found = 1

이 두 가지 방법에서 얻은 시간을 비교하기 위해 그래프로 그려 보겠습니다:

plt.plot(matrix_numpy_times, color="blue")
plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

이는 실제 소요 시간이지만, 논의의 편의를 위해 인접한 몇 개의 점/행렬 차원에 대한 평균을 구하여 이러한 곡선을 평활화해 보겠습니다. 이 작업은 아래에서 수행됩니다:

smooth_numpy_times = []
smooth_krylov_times = []

# Choose the number of adjacent points over which to average forward;
# the same will be used backward.
smooth_steps = 10

# We will do this smoothing for all points/matrix dimensions
for i in range(len(matrix_krylov_tot_times)):
    # Ensure we don't exceed the boundaries of our lists
    start = max(0, i - smooth_steps)
    end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)

    # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov
    # and the NumPy calculations.
    smooth_count = 0
    smooth_numpy_sum = 0
    smooth_krylov_sum = 0

    for j in range(start, end):
        smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]
        smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]
        smooth_count = smooth_count + 1

    # Appending the averaged adjacent values to our new smooth lists
    smooth_numpy_times.append(smooth_numpy_sum / smooth_count)
    smooth_krylov_times.append(smooth_krylov_sum / smooth_count)

plt.plot(smooth_numpy_times, color="blue")
plt.plot(smooth_krylov_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (smoothed, microsec)")  # Add y-axis label
plt.show()

Output:

크릴로프 하위 공간의 구축에 필요한 시간은 처음에는 눔피의 전체 대각선화에 필요한 시간을 초과한다는 점에 유의하세요. 그러나 행렬의 크기가 커질수록 크릴로프 방식이 유리해집니다. 이는 허용 오차를 낮추더라도 마찬가지이지만 행렬 크기가 커질수록 이점이 더 커집니다. 이 부분은 자세히 살펴볼 가치가 있습니다.

수치 대각선화의 시간 복잡도는 $O(n^3)$ (알고리즘에 따라 약간의 차이가 있음)입니다. $n$ 벡터의 직교 기저를 생성하는 시간 복잡도도 $O(n^3)$ 입니다. 따라서 크릴로프 방법의 장점은 $\it{some}$ 정규 기저를 사용하는 것이 아니라 관심 있는 고유값을 효과적으로 골라내는 특정 정규 기저를 사용하는 것과 관련이 있습니다. 우리는 이미 이 강의의 첫 번째 섹션에서 증명의 스케치에서 이를 살펴봤으며, 이는 크릴로프 메서드의 수렴 보장에 매우 중요합니다.

지금까지의 진행 상황을 살펴보겠습니다:

• 매우 큰 행렬의 경우 크릴로프 부분 공간 방법은 기존의 대각선화 알고리즘보다 필요한 오차 범위 내에서 대략적인 고유값을 더 빠르게 구할 수 있습니다.
• 이렇게 매우 큰 행렬의 경우, 크릴로프 부분공간 생성은 크릴로프 부분공간 방법에서 가장 많은 시간이 소요되는 부분입니다.
• 따라서 크릴로프 서브스페이스를 효율적으로 생성하는 방법은 매우 유용할 것입니다. 바로 이 지점에서 양자 컴퓨터가 등장합니다.

이해도 점검

위의 대각화 시간에 대한 행렬 차원의 평활화 그래프를 참조하십시오.

(a) 이 그래프에 따르면, 크릴로프(Krylov) 방법이 어느 정도의 행렬 차원부터 더 빨라지기 시작했는가?

(b) 크릴로프(Krylov) 방법이 더 빨라지는 그 기준치는 계산의 어떤 측면에서 달라질 수 있는가?

• 답변

  (a) 계산을 다시 수행하면 결과가 달라질 수 있지만, 크릴로프(Krylov) 방법은 차원이 약 80~85가 되면 계산 속도가 빨라집니다.

  (b) 가능한 답은 여러 가지가 있다. 중요한 요소로는 우리가 요구하는 정밀도와 대각화 대상 행렬의 희소성이 있습니다.

3. 시간 진화를 통한 크릴로프

지금까지 설명한 모든 작업은 고전적으로 수행할 수 있습니다. 그렇다면 언제 어떻게 양자 컴퓨터를 사용할 수 있을까요? 매우 큰 행렬의 경우, 크릴로프 방식은 긴 계산 시간과 많은 양의 메모리가 필요할 수 있습니다. $|v\rangle$ 에서 $H$ 의 행렬 연산에 필요한 시간은 최악의 경우 $O(N^2)$ 와 같이 확장됩니다. 벡터에 희소 행렬을 곱하는 것(고전적인 크릴로프 유형 솔버의 일반적인 경우)도 $O(N)$ 과 같은 시간 복잡도 스케일링이 있습니다. 이 작업은 하위 공간에서 원하는 모든 벡터에 대해 수행됩니다. 하위 공간 차원 $r$ 은 일반적으로 $N$ 의 상당 부분이 아니며 종종 $\log(N)$ 과 같이 스케일링됩니다. 따라서 모든 벡터를 생성하면 최악의 경우 $O(N^2 \log(N))$ 처럼 스케일링됩니다. 직교화와 같은 다른 단계도 있지만, 이것이 염두에 두어야 할 주요 스케일링입니다.

양자 컴퓨팅을 사용하면 문제의 속성을 변경하여 필요한 시간과 리소스의 규모를 결정할 수 있습니다. 전반적으로 행렬 크기( $N$ )에 의존하는 대신 해밀턴을 구성하는 샷 수와 비교대 폴리 항의 수와 같은 것을 볼 수 있습니다. 어떻게 작동하는지 살펴보겠습니다.

3.1 시간에 따른 변화

양자 상태를 시간 진화시키는 연산자는 $e^{-iHt/\hbar}$ (특히 양자 컴퓨팅에서는 $\hbar$ 을 표기에서 빼는 것이 매우 일반적입니다)라는 점을 기억하세요. 이러한 연산자의 지수 함수를 이해하고 심지어 실현하는 한 가지 방법은 테일러 급수 확장을 살펴보는 것입니다. 이 연산은 어떤 초기 벡터 $|v\rangle$ 에 작용하여 초기 상태에 적용된 $H$ 의 거듭제곱이 증가하는 항의 합을 산출합니다. 초기 추측 상태를 시간 진화시켜 크릴로프 하위 공간을 만들 수 있을 것 같습니다!

주의할 점은 실제 양자 컴퓨터에서 시간 진화를 실현하는 것입니다. 해밀턴의 많은 용어는 서로 통용되지 않습니다. 따라서 $e^{-iZ}$ 같은 간단한 지수 연산자는 간단한 회로에 해당하지만, 일반적인 해밀턴 연산자는 그렇지 않습니다. 그리고 지수에는 통근이 아닌 용어가 포함되어 있기 때문에 숫자처럼 단순히 지수를 단순 곱으로 분해할 수 없습니다.

따라서 이것은 사소한 일이 아니지만 양자 컴퓨팅에서 잘 연구된 과정입니다. 양자 컴퓨터에서 트로터라이제이션이라는 프로세스를 사용하여 시간 진화를 수행하는데, 이는 그 자체로 풍부한 주제입니다 [10]. 그러나 매우 높은 수준에서는 시간 진화를 매우 작은 단계(예: $m$ 단계의 크기 $dt$ )로 나누면 용어의 비공변성의 영향을 제한할 수 있습니다.

여기서 $dt = t/m$.

고전적인 맥락에서 H의 거듭제곱을 직접 사용하여 생성한 차수 r의 크릴로프 부분공간을 "파워 크릴로프 부분공간"이라고 부르겠습니다.

이제 단일 시간 진화 연산자 $U \equiv e^{-iHt}$ 를 사용하여 유사한 공간을 생성합니다. 이를 "단일 크릴로프 공간" $\mathcal{K}_U^r$ 이라고 부릅니다. 우리가 고전적으로 사용하는 파워 크릴로프 아공간 $\mathcal{K}_P^r$ 은 $H$ 이 단일 연산자가 아니므로 양자 컴퓨터에서 직접 생성할 수 없습니다. 단일 크릴로프 부분공간을 사용하면 파워 크릴로프 부분공간과 유사한 수렴 보장, 즉 초기 상태 $|v\rangle$ 가 기하급수적으로 소멸하지 않는 실제 기저 상태와 겹치고 고유값 사이에 충분한 간격이 있는 한 기저 상태 오차가 효율적으로 수렴한다는 것을 보여줄 수 있습니다. 컨버전스에 대한 보다 정확한 논의는 참고 [1] 을 참조하세요.

여기서 $U$ 의 거듭제곱은 서로 다른 시간 단계가 됩니다( $U$ 의 $k^\text{th}$ 거듭제곱은 시간만큼 앞으로 나아갑니다 $k \times dt$ ). 총 시간 $k dt$ 동안 시간 진화하는 하위 공간의 요소에 $|\psi_k\rangle$ 이라는 레이블을 붙일 수 있습니다.

해밀턴 H를 단일 크릴로프 부분공간인 $\mathcal{K}_U^r$ 에 투영할 수 있습니다. 즉, $\mathcal{K}_U^r$ 기준으로 $H$ 의 각 행렬 요소를 계산합니다. 이 예상 매트릭스를 $\tilde{H}$ 이라고 합니다.

3.2 양자 컴퓨터에서 구현하는 방법

$\tilde{H}$ 의 행렬 요소는 양자 컴퓨터를 사용하여 추정할 수 있는 기대값 $\langle \psi_m |H| \psi_n\rangle$ 으로 주어집니다. $H$ 은 서로 다른 큐비트에 대한 폴리 연산자의 합으로 쓸 수 있으며, 모든 폴리 연산자를 동시에 측정할 수 있는 것은 아니라는 점에 유의하세요. 폴리 용어를 통근 용어 그룹으로 분류하고 모든 용어를 한 번에 측정할 수 있습니다. 하지만 모든 용어를 다루기 위해서는 이러한 그룹이 많이 필요할 수 있습니다. 따라서 기간을 분할할 수 있는 별개의 출퇴근 그룹( $N_\text{GCP}$ )의 수가 중요해집니다.

여기서 $P_\alpha$ 는 $P_\alpha \sim IZIXII...YZXIX$ 형태의 파울리 벡터이거나, 서로 교환 관계를 이루는 이러한 파울리 벡터들의 집합이다. $H$ 를 측정 가능한 연산자들의 이러한 합으로 표현할 수 있다는 점을 고려할 때, $\tilde{H}$ 의 행렬 요소에 대한 다음 표현들은 IBM Quantum 기본 추정기를 사용하여 구현할 수 있다.

여기서 $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ 은 단일 크릴로프 공간의 벡터이고 $t_n = n dt$ 은 선택한 시간 단계 $dt$ 의 배수입니다. 양자 컴퓨터에서 각 행렬 요소의 계산은 양자 상태 간의 중첩을 얻을 수 있는 모든 알고리즘으로 수행할 수 있습니다. 이 레슨에서는 하다마드 테스트에 중점을 두겠습니다. $\mathcal{K}_U$ 의 차원이 $r$ 이라고 가정하면, 하위 공간에 투영된 해밀턴은 $r \times r$ 의 차원을 갖습니다. $r$ 이 충분히 작으면(일반적으로 $r<<100$ 이 고유값 추정치의 수렴을 얻기에 충분합니다), 투영된 해밀토니안 $\tilde{H},$ 을 고전적으로 쉽게 대각선화할 수 있습니다. 하지만 크릴로프 공간 벡터의 비직교성 때문에 $\tilde{H}$ 을 직접 대각선으로 만들 수는 없습니다. 겹치는 부분을 측정하고 행렬을 구성해야 합니다 $\tilde{S}$

이를 통해 비직교 공간에서 고유값 문제(일반화된 고유값 문제라고도 함)를 해결할 수 있습니다

그런 다음 이 일반화된 고유값 문제의 해를 살펴봄으로써 $H$ 의 고유값과 고유 상태의 추정치를 구할 수 있습니다. 예를 들어, 접지 상태 에너지의 추정치는 가장 작은 고유값 $E$ 과 해당 고유 벡터 $\vec{c}$ 에서 접지 상태를 구하여 얻을 수 있습니다. $\vec{c}$ 의 계수는 $\mathcal{K}_U$ 에 걸쳐 있는 여러 벡터의 기여도를 결정합니다.

일반화된 고유값 문제

$\tilde{H}$ 을 단순히 대각선화할 수 없는 이유는 무엇인가요? $\tilde{S}$ 에는 크릴로프 기저의 기하학적 정보(아주 특별한 경우를 제외하고는 모두 직교하지 않음)가 포함되어 있으므로 $\tilde{H}$ 자체는 전체 해밀턴의 투영을 설명하지 않으므로 고유값은 전체 해밀턴의 고유값과 특별한 관계가 없으며 임의의 값이 될 수 있습니다. 일반화된 고유값 문제를 풀려면 전체 해밀턴을 크릴로프 공간에 투영한 것에 해당하는 근사 고유값과 고유 벡터를 구해야 합니다.

그림은 서로 다른 양자 상태 간의 중첩을 계산하는 데 사용되는 방법인 수정된 하다마드 테스트의 회로 표현을 보여줍니다. 각 행렬 요소 $\tilde{H}_{i,j}$ 에 대해 상태 $\vert \psi_i \rangle$, $\vert \psi_j \rangle$ 사이의 Hadamard 테스트가 수행됩니다. 이는 그림에서 행렬 요소의 색 구성표와 해당 $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ 연산으로 강조 표시되어 있습니다. 따라서 투영된 해밀토니안 $\tilde{H}$ 의 모든 행렬 요소를 계산하려면 크릴로프 공간 벡터의 가능한 모든 조합에 대한 일련의 하다마드 테스트가 필요합니다. 하다마르 테스트 회로의 상단 와이어는 X 또는 Y 기준으로 측정되는 안실라 큐비트이며, 그 기대값에 따라 상태 간 중첩 값이 결정됩니다. 아래쪽 와이어는 시스템 해밀턴의 모든 큐비트를 나타냅니다. $\text{Prep} \; \psi_i$ 연산은 안실라 큐비트의 상태에 의해 제어되는 $\vert \psi_i \rangle$ 상태의 시스템 큐비트를 준비하고( $\text{Prep} \; \psi_j$ )와 마찬가지로 $P$ 연산은 시스템 해밀턴 $H = \sum_i P_i$ 의 파울리 분해를 나타냅니다. 이를 양자 컴퓨터에서 구현하는 방법은 아래에 자세히 설명되어 있습니다.

4. 양자 컴퓨터에서의 크릴로프 양자 대각화

이제 실제 양자 컴퓨터에서 크릴로프 양자 대각선화를 구현해 보겠습니다. 몇 가지 유용한 패키지를 가져오는 것부터 시작하겠습니다.

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import warnings

from qiskit.quantum_info import SparsePauliOp, Pauli
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

# from qiskit.providers.fake_provider import Fake20QV1
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator, Batch

import itertools as it

warnings.filterwarnings("ignore")

방금 설명한 일반화된 고유값 문제를 해결하기 위해 아래 함수를 정의합니다.

def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in our Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of our Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem


    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array([vec for val, vec in zip(s_vals, s_vecs) if val > threshold])
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]

적어도 초기 벤치마킹에서는 컨버전스 동작을 확인하기 위해 정확한 고전적인 솔루션을 아는 것이 유용합니다. 아래 함수는 해밀턴과 큐비트 수를 인수로 사용하여 해밀턴의 기저 상태 에너지를 계산합니다.

def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]  # list all of the possible sets of n_exc indices of 1s in n_exc-particle states

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

4.1 1단계: 문제를 양자 회로 및 연산자로 매핑하기

이제 해밀턴을 정의하겠습니다. 위의 함수는 해밀턴을 인수로 받아 기저 상태만 반환한다는 점에서 위의 함수와 다르며, 고전적인 방식으로 동작합니다. 여기서 정의하는 해밀턴은 모든 에너지 고유 상태의 에너지 레벨을 결정하며, 이 해밀턴은 파울리 연산자를 사용하여 구성하고 양자 컴퓨터에서 구현할 수 있습니다.

공간에서 어떤 방향이든 가질 수 있는 스핀의 사슬에 해당하는 해밀턴을 '하이젠베르크 사슬'이라고 합니다. $i^\text{th}$ 스핀은 가장 가까운 이웃( $(i-1)^\text{th}$ 및 $(i+1)^\text{th}$ 스핀)의 영향을 받을 수 있지만 더 먼 이웃의 영향을 받지는 않는다고 가정합니다. 또한 스핀이 서로 다른 축을 가리킬 때 스핀 간의 상호작용이 달라질 가능성도 허용합니다. 이러한 비대칭성은 예를 들어 스핀이 내장된 결정 격자의 구조로 인해 때때로 발생합니다.

# Define problem Hamiltonian.
n_qubits = 10
# coupling strength for XX, YY, and ZZ interactions
JX = 1
JY = 3
JZ = 2

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [
    (term, JZ)
    if term.count("Z") == 2
    else (term, JY)
    if term.count("Y") == 2
    else (term, JX)
    for term in H_int
]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)