ベスパリブ

プログラミングを主とした日記・備忘録です。ベスパ持ってないです。

AHC026 参加記録

トヨタ自動車プログラミングコンテスト2023#6(AtCoder Heuristic Contest 026)に参加した記録です。

方針1:貪欲法

4時間の短期間コンテストなので、まずは正の得点を得られる簡単な貪欲ができないかを優先した。

とりあえずビジュアライザを動かしてみて、どういう貪欲法がいいかを考える。

ターゲットの箱より上のすべての箱を、最も使わなそうな山を選んで、そこに積んでいく貪欲法をした。

「最も使わなそうな山」は、ターゲットの箱となるものが時間的に最も遠い山を選んだ。

seed0のビジュアライザ

提出コード:https://atcoder.jp/contests/ahc026/submissions/47302082

得点:1,356,657

この得点の人が多かったので、似たような解法をやってる人多そう。

方針2:ソートする?

箱を動かす時、なるべく上から昇順になるようにしたほうが、同じ体力で効率よく動かせるのでは?

イメージ図

結果的に言えば、重さ+1の体力を消費するのでそんなことなかった。

提出コード:https://atcoder.jp/contests/ahc026/submissions/47305951

得点:1,324,254

得点下がっちゃった。

方針3:ビームサーチ

多様性が確保できそうなのでビームサーチいけるかなと思って実装。永遠にバグらせて終了。

終結

方針1のコードが最終提出になった。

得点:1,356,657

順位:426位

強い人の方針

「片方に集める」「寄せる」みたいなのはAHCの典型なので最初そうしようかなと思ったのだけど、後半すっかり忘れていた。

山が10個あるからうまくソート状態を作ればいいのか?(ハノイの塔のイメージ?)

感想

  • 最初の貪欲の方針は悪くなかったと思う。
  • ビームサーチバグらせたのは残念。実装方法のテンプレートをNotionにまとめておこう。
  • でもビームサーチあんまりよいスコアでないらしい
  • AHC楽しいので毎週開催してほしい。
  • レートが単調増加は嬉しい。

AHC025 参加記録

概要

AtCoder Heuristic Contest 025に参加した記録です。

珍しいインタラクション問題。

方針1:ナイーブ実装

とりあえずシンプルな解法を考えて実装するところから着手した。以下のような解法になる。

  1. 最初、各グループは1つずつアイテムを持つ。残ったアイテムはすべてグループ0が持つ。
  2. s=0とする。
  3. グループsが他のグループに順番に1個ずつアイテムをランダムに渡していく。
  4. グループsがグループiにアイテムを渡したとき、グループiがグループsより重さが大きくなったら、s = iになって、再度3を繰り返す。

実装上の注意点として、

  • アイテムが0個のグループを作ってはいけない
    • 他のどのアイテムを足し合わせても敵わないような激重アイテムが1個ある場合、そのアイテム1個を唯一持つグループが出来上がる。このグループからその激重アイテムを取るとアイテムが0個になってしまうので取ってはいけない。

シード1の結果

提出コード:https://atcoder.jp/contests/ahc025/submissions/46658987

得点:1,163,992,256

方針2:クイックソート実装

方針1の問題点は、「アイテムを1個ランダムに選んで渡す」とき、重いアイテムを渡す可能性があることにある。できれば軽いアイテムをやりとりして各グループの重さがならされるようにしたい。

ここで、クエリに2つのアイテムの大小関係を尋ねまくると、(アイテムの重さは不明のままだが)いずれアイテムをソートできることに気づく。アイテムがソートできた状態になると、各グループの一番軽いアイテムがわかるので、軽いアイテムを選択してやりとりすることで微調整できるようになる。

ソートアルゴリズムクイックソートを実装した。

実装上の注意点として、

  • クイックソートは最悪クエリ回数が(N-1)N/2 回かかる
    • クイックソートでクエリ回数を使い果たさないようにしなければならない。今回は実装が簡単のため N*10 <= Q のときクイックソートするのに十分なクエリ回数があると判断するプログラムとした。十分なクエリ回数がない場合、ナイーブ実装に切り替えた。

シード1の結果

提出コード:https://atcoder.jp/contests/ahc025/submissions/46693668

得点:811,781,427

今回の問題は得点が小さいほど良いスコアとなるので、改善されたと言える。

方針3:クイックソート打ち切り

方針2の改良。クイックソートに使えるクエリ回数上限を設定できるようにして、上限に達したら途中でソートを打ち切る。打ち切ると順番が中途半端にわからないアイテムがあるが気にしない。

アイテムの交換に使える最低クエリ回数を D*3 くらい(各グループ3回くらいアイテムのやりとりをするくらい)は確実に残すようにして、あとはクイックソートに使えるように設定した。もちろんクイックソートが早く終わればアイテムの交換に使えるクエリ回数は増える。

これにより方針1のナイーブ実装の実行を完全に排除できた。(運次第だと思うが、ナイーブ実装のほうが良いスコアのサンプルがちらほらあったが……)

提出コード:https://atcoder.jp/contests/ahc025/submissions/46697190

得点:586,051,759

このあと、アイテムの交換に使える最低クエリ回数の値をいろいろ試したところ、D*5 が最もスコアが良かったので最終的にそうした。

方針4:グループの重さのクイックソート実装

Seed 29を考察した。

方針3のシード29の結果

シード29は完全にソートされているにも関わらず、あまり良くないスコアになっている。

たとえば、左から5番目のグループは軽いアイテムをたくさん持っているにも関わらず、それらを貯め込んで一番重いグループになっている。こいつのグループが他のグループに軽いアイテムを配ってくれればもっと良いスコアになるはず。

分析の結果、全1423ターンのうち、509ターン目でソートが完了し、残り約900 ターンでアイテム交換をしていた。アイテム交換に十分のターン数があるにも関わらず、スコアは良くない結果になっている。

もっと深く考察するために、途中結果を出力して確認してみる。

ソート完了直後の状態

ソート完了直後は上の状態となった。

1回目のアイテム交換

一番左と、左から3番目のアイテムの交換が行われた。

2回目のアイテム交換

一番左と、左から4番目のアイテムの交換が行われた。

これを見ると、重いグループ同士で無駄な交換が起きているように見える。なので、クエリ回数が十分に残っている場合(最悪交換回数の(D-1)*D/2以上の場合)、各グループの重さでソートし、「重いグループと軽いグループで交換を行う」ということをすれば良さそうである。

すなわち今回の方針では、グループの重さのクイックソートを実装した。

提出コード:https://atcoder.jp/contests/ahc025/submissions/46748553

得点:602,827,302

スコアが悪くなってしまった。

方針5:方針4の改良(残りクエリ回数を十分取るときだけグループのソートをする)

方針4では、残りクエリ回数が十分にあるときだけグループのクイックソートをして、そうでないときは方針3に切り替えるように調整した。それでもスコアが悪かったので、何がいけないのか出力ファイルを比較して原因を探った。

悪くなったサンプルを比較してみる。

以下の分数は、(アイテムのソートで使ったクエリ回数)/(全クエリ回数)。

シード7の比較。

278/294

左が方針3、右が方針4

シード8の比較。

175/197

左が方針3、右が方針4

シード16の比較。

206/221

左が方針3、右が方針4

良くなったサンプルを比較してみる。

シード1の比較。

306/1120

左が方針3、右が方針4

シード18の比較。

368/1028

左が方針3、右が方針4

他にも比較したが割愛。

比較してわかったのが、 (アイテムのソートで使ったクエリ回数)/(全クエリ回数)の値が大きい (=アイテムのソート後に残ったクエリ回数が少ない)と、グループのソートに使えるクエリ回数が少なすぎて方針3より悪いスコアになってしまうらしいということ。

つまり残りクエリ回数が十分にあるときだけグループのクイックソートをしていたつもりが、まだまだ全然足りなかったということらしい。

残りクエリ回数がQ/2以上のときだけグループのクイックソートをするように変更した。

提出コード:https://atcoder.jp/contests/ahc025/submissions/46748977

得点:543,035,466

方針3よりスコアが良くなった。

上記の方針5では、方針3よりサンプル100個中、

  • 5個スコアが悪くなった
  • 31個スコアが改善された
  • 64個はスコアが変わらなかった

方針6:残りクエリ回数に応じて確率的にグループの重さが切り替わるようなアイテム交換を許す

方針5では、64個はスコアが変わらないのであった。これは、アイテムのソート後に残りクエリ回数が少ないので方針3に切り替わった個数が64個ということであった。

方針3の問題点はなんだろうかと、いくつか悪いスコアのものをピックアップしてみた。

seed62: 157/185

シード62の結果

seed66: 438/482

シード66の結果

ビジュアライザを動かしてみると、残りクエリ回数が少ないときにもグループ間の重さが入れ替わるような重いアイテムの交換を許してしまっていた。

残りクエリ回数がまだ残ってるときはグループの重さが入れ替わってもいいが、残りクエリ回数が少ないときはグループの重さが入れ替わるようなことはせず、軽いアイテムのやりとりだけで済ましてほしい。イメージ的には『安定』したやりとりをしてほしいという欲求がある。

なので、残りクエリ回数に応じて確率的にグループの重さが切り替わるようなアイテム交換を許す実装をした。気持ち的には焼き鈍し法である。

方針7のシード62の結果

方針7のシード66の結果

かなり良くなった感じがしたので、提出した。

提出コード:https://atcoder.jp/contests/ahc025/submissions/46841307

得点:492,973,943

終結

方針6を最終結果として提出した。

A:58,654,377,877

A-Final:300,672,027,666

順位:458

ここまでの感想

  • 方針5を実装した時点での相対スコアが62Gで、順位は380位くらいだった。この得点帯は激戦区で、寝て起きたら順位が結構下がっている。
  • 方針5を実装した時点で頭打ちしてると思ったのだが、上位陣が90Gで途方もなかった
  • 「山登り、焼きなまし、ビームサーチは絶対にできないよな」「いやしかし…」とぐるぐる考えてた
  • 方針6は発想としてはかなり焼き鈍し法なのだが、これに至るまでかなり時間を要した。

公式解説

公式解説動画のリンク

X(Twitter)で調べると強い人はソートして一番重いものと一番軽いものを交換している人が多いイメージ。似たようなことは自分もしているのだが、自分の解法では一番重いやつと軽いやつを交換、2番目に重いやつと軽いやつを交換、…としていたのが良くなかったのか?

公式解説の動画を見てみる。

  • 解法1:クエリを使ってがんばって重さを推測する。その後クエリを使わずに最適化問題を解く
  • 解法2:焼きなまし法(近傍解を求めて、その近傍解が良くなったら採用。悪くなったら悪くなった度合いに応じて確率的に採用)は、今回の天秤を使ったクエリでは「どれだけ悪くなった」かがわからないので難しい(「悪くなった」はわかるが「どれだけ悪くなったか」はわからない)。が、「良くなった」のは判定できるので、山登り法が可能
    • 初期解を適当に作って、山登り法をする
  • 解法3:初期解を貪欲法で作る

公式解説では解法2の山登りを解説してくれた。

山登り法

2つのグループ間でアイテム交換をするの図(公式解説動画から引用)

2つのグループの大小関係がわかっている状態で、アイテム交換をする。アイテム交換後、2つのグループの大小関係が変わっていなかったら、「スコアが良くなった」といえる(2つのグループ間の分散が小さくなっているので)。

基本的にはこれをひたすら山登っていけばいい。

この解法のシンプルな手法ではソートすら必要なくて、アイテムをランダムに1個選び、それをどこのグループに入れるかもランダムに選ぶ。アイテム交換してスコアが良くなるなら交換し、スコアが良くならないなら交換しない。

初期解は各グループ同じアイテムの個数になるように適当に割り振る。

これをすると321位がとれるらしい。

うーーーーーーーーーん。シンプルな解法に見えるが、これでも十分自分より良い点数が取れるところを見ると、考察力が足りないと感じる。やってることは「言われてみれば当然といえば当然」のことなのだが、思いついてないんだよなぁ。

アイテム交換方法2

もうひとつのアイテム交換方法として、グループAとグループBからそれぞれアイテムを1個ずつ取り出して、取り出したアイテムをアイテムa、アイテムbとする。

このとき、「グループA < グループB」「アイテムa < アイテムb」ならば、アイテムaはグループBに入れ、アイテムbはグループAに入れるとスコアが良くなる。

アイテムaはグループBに、アイテムbはグループAに入れるのが良い

このようなアイテム交換はクエリ2回の発行で行える。

「この交換法も50%の確率でする」みたいに組み合わせると218位に入れるらしい。

キャッシュする + 推移律のグラフを作る

一度した交換方法はキャッシュして覚えて、2回もしないようにすると地味に良くなる。

アイテム交換方法2をすると、アイテムの大小関係のデータが蓄積していく。

a < b のような大小関係のグラフを作っておくと、a < b と b < c => a < c のように、クエリを発行しなくても推移律で新たに大小関係が判明する。

推移律からa < cを導く

もっとやろうと思えば以下のように導くこともできる。

a<b と b+c<d から、a+c < d も判明する

こういう不等式系がたくさんあったときにすべて満たす解を求める方法にLP(線形計画法)があるらしい。

それをしなくても「アイテムa+cとアイテムdの交換」といった複数個の交換が可能になる。

これを実装すると200位になる。

グループのソートを実装する

これまでの方法の問題点は、交換するグループをランダムに選んでいたのが問題であった。

そこでグループのソートを実装すると、どのグループがアイテムを交換すべきなのかがわかる。

一番重いグループと一番軽いグループがアイテム交換するのが良い

アイテム交換後のグループのソートは単純にやればO(DlogD)。けどこれは以下のようにして改善できる。

一番重いグループをグループA、一番軽いグループをグループBとする。アイテム交換後、グループAとグループB以外の他のグループはソート済みなので、「グループA(B)は、ソート済みのグループのどこに入るか」を2分探索すればO(logD)で済む。確かに……。

これをすると49位になる。

公式解説の感想

  • 公式解説を見ると、「なぜこれが思いつかなかったのか……」といった感想になった。ソートなどそれっぽいことは思いついて実装できているのだが、初期の解法の土台が良くないのでスコアが伸びなかった印象。論理的な道筋が見えていないのというか。練度が足りない。

AHC024 参加記録

概要

丸紅プログラミングコンテスト2023(AtCoder Heuristic Contest 024)に参加しました。4時間の短期間コンテストでした。

ほぼ何もできなかったといってもいい。

最初の方針

連結成分をノード、「隣接する」というのを辺と考えるとグラフができるので、そのグラフからグリッドを作ろうと思いましたが、実装が大変で2時間経過。スジが悪いと思い直してビジュアライザをいじっていると、与えられたグリッドを白に塗っていくのほうが実装が簡単なことに気づきました。

かなりバグ散らかしつつ終了10分前くらいにかなりナイーブな貪欲ができたのでそれを提出。終了です。

atcoder.jp

seed0のビジュアライザ

白に塗れるところは貪欲に白に塗っていく。

真ん中あたりのエリアが大きいですね。これらを隣接条件等を崩さないように圧縮していければなあと思いつつ終了しました。あと4時間ほしい。

上手い人の解法

中心に寄せる貪欲らしいです。貪欲でもここまでいけるらしい。

https://twitter.com/Jiro_tech15/status/1705889677910348053

左上に寄せていき、最後にぎゅっと圧縮する手法。おお、と思いました。AHC系は「はしっこに寄せる」という発想はよく見る気がします。

消せる行や列を消すという発想。「隣接条件を崩さずにどう圧縮するか」というのを考えたとき、たしかに「1行ごっそり消せるなら消す」というのは実装ラクそう。天才や……

反省

AHCやることメモに「まずビジュアライザをいじる」って書いてあるだろ。なぜそれをしないんだ。

短時間コンテストは実装のしやすさ・速さがかなり重要だと感じました。

あと最近、C++でクラスを書くときはメンバ変数にm_プレフィックスをつけているのだがいい感じ。メンバにアクセスするときは通常this->を使うと思うのだが、長いのがやだ。

Chromeでデベロッパーツールを開くとwasmが遅くなる現象

概要

wasm生成の練習のため、実践Rustプログラミング入門を参考にしてマンデルブロ集合の計算のwasmとjsで速度比較するWebページを作成しました。

taketakeyyy.github.io

しかし開発中、wasmの速度が遅くなったり速くなったりする挙動に悩まされました。色々調べた結果、Chromeデベロッパーツールを開いているか否かが原因だったようです。

比較

デベロッパーツールを開いてない状態

デベロッパーツールを開いている状態

なんでこうなるのか全くわかっていません🤷後述追記した。

この現象はlocalhostや本番環境関係なく発生しています。

デベロッパツールを開いてconsole.logで時間計測する際には注意が必要になりそうです。

「wasm 遅い」で検索すると色々出てきますが、こういうのも原因のひとつになっているんじゃないかと思いました。

2023/09/21 0:56 追記

Twitter(X)でご指摘いただきました。ありがとうございます。

nmi.jp

Chrome において、DevTools をただ開くと wasm の実行が極めて遅くなる現象が確認されております。これはバグではなくて仕様なのですが、結構複雑な仕様で、

DevTools を開くと TurboFan の最適化がキャンセルされて遅くなる しかし Profile タブで Profile を取るときには再度 TurboFan の最適化が適用される という挙動になっております。Chrome では、ただ DevTools を開いているだけで wasm の実行がかなり遅くなります。大変気づきにくいトリガーなので、wasm 系の開発をされる時にはお気をつけください。

そういう仕様らしいです。

その他

2人の子供問題のパラドックスを利用したギャンブル

前提条件

A「ある家庭に2人の子供がいて、1人が男の子の場合、もう1人の子が女の子である確率は?」

B「ある家庭で男の子が生まれる確率と女の子が生まれる確率はそれぞれ1/2ずつで考えていい?」

A「いいよ」

B「うーん。もう1人の子が男の子か女の子かは独立だから、1/2じゃない?」

A「残念。ある家庭に2人の子供がいる場合の(年上, 年下)のペアは(男, 男)、(男, 女)、(女, 男)、(女, 女)の4通りで、1人が男の子で確定したときに(女, 女)は除外される。だから残りの3通りのうち、もう1人が女の子のパターンは(男, 女)、(女, 男)の2通りだから答えは2/3だよ」

B「なるほどなぁ」

ギャンブルに持ち込む

A「よし、今の2人の子供の問題と同じ状況を再現して賭けをしようよ」

B「どんな賭け?」

A「僕は2枚のコインを同時に投げる。1枚のコインはみんなに見えるように。もう1枚のコインは表か裏どちらが出たのか僕も君もわからないようにする。」

B「それで?」

A「見えてるコインが表になったときに賭けは始まるんだ。『見えてないコインがであるか賭けてください』って。裏だったら君の勝ち。そうじゃなかったら僕の勝ちだ。さっきの2人の子供と同じ状況だから、見えてないコインが裏の確率は2/3だよね?」

B「そうかも

A「見えてるコインが裏になったときはコインを投げ直すよ」

B「まあいいでしょう」

A「当たったときの賞金は僕が1000円出す。当たる確率が2/3だから賞金の期待値は666円だね。参加費は1回600円でどう?これは君にバチクソ有利な賭けだよ」

B「やろうやろう」

# coding: utf-8
import random

entry_fee = 600 # 参加費
prize = 1000 # 賞金
bet_num = 0 # 賭けを行った回数
iter_max = 1000 # コインを投げた回数
win_b = 0 # B君が勝った回数

for i in range(iter_max):
    seen = random.randrange(2) # 見えているコイン(0:表, 1:裏)
    unseen = random.randrange(2) # 見えてないコイン(0:表, 1:裏)
    if seen==0:
        # 見えてるコインが表になったら賭けをする
        bet_num += 1
        if unseen==1:
            # 見えてないコインが裏だった場合、B君の勝ち
            win_b += 1

print("賭けを行った回数: {}".format(bet_num))
print("B君の勝率: {}".format(win_b/bet_num))
print("B君が払った参加費の合計: {}".format(entry_fee*bet_num))
print("B君が貰った賞金の合計: {}".format(prize*win_b))
print("B君の収支: {}".format(prize*win_b - entry_fee*bet_num))
賭けを行った回数: 524
B君の勝率: 0.5305343511450382
B君が払った参加費の合計: 314400
B君が貰った賞金の合計: 278000
B君の収支: -36400

1時間後

B「3万円も負けちゃった」

A「そういうときもあるよ」

B「わかった。見えてるコインが裏になったときにコインを投げ直しているのがイカサマだ。試行をごまかそうとしたね?」

A「じゃあ見えているコインが裏になったときも賭けをしよう。そのときは『見えてないコインがであるか賭けてください』だね。表だったら君の勝ち。そうじゃなかったら僕の勝ちだ。これは見えているコインが表のときと同じ賭けだから問題ないよね?」

B「それならいいかも」

A「じゃあ続けるよ」

# coding: utf-8
import random

entry_fee = 600 # 参加費
prize = 1000 # 賞金
bet_num = 0 # 賭けを行った回数
iter_max = 1000 # コインを投げた回数
win_b = 0 # B君が勝った回数

for i in range(iter_max):
    seen = random.randrange(2) # 見えているコイン(0:表, 1:裏)
    unseen = random.randrange(2) # 見えてないコイン(0:表, 1:裏)
    if seen==0:
        # 見えてるコインが表になったら賭けをする
        bet_num += 1
        if unseen==1:
            # 見えてないコインが裏だった場合、B君の勝ち
            win_b += 1
    # 【追加】
    elif seen==1:
        # 見えてるコインが裏になったときも賭けをする
        bet_num += 1
        if unseen==0:
            # 見えてないコインが表だった場合、B君の勝ち
            win_b += 1

print("賭けを行った回数: {}".format(bet_num))
print("B君の勝率: {}".format(win_b/bet_num))
print("B君が払った参加費の合計: {}".format(entry_fee*bet_num))
print("B君が貰った賞金の合計: {}".format(prize*win_b))
print("B君の収支: {}".format(prize*win_b - entry_fee*bet_num))
賭けを行った回数: 1000
B君の勝率: 0.512
B君が払った参加費の合計: 600000
B君が貰った賞金の合計: 512000
B君の収支: -88000

さらに1時間後

B「なるほど、ひとつだけわかったことがある」

A「それはなに?」

B「君との友情はもう壊れたってこと」

A「😅トホホ・・・」

解説

  • 2人の子供問題のパラドックスの本質は、「確率は1/2ではなく、実は確率は2/3でした」というより、「問題文が曖昧なので1/2と解釈することも2/3と解釈することもできる」点にあるだろう
  • その曖昧さを利用して、相手に「実は確率は2/3」であることを納得させる
  • あとは似て非なる「本当の確率は1/2」のギャンブルに持ち込み、ボロ儲けする

B君がやるべきだったギャンブル

# coding: utf-8
import random

entry_fee = 600 # 参加費
prize = 1000 # 賞金
bet_num = 0 # 賭けを行った回数
iter_max = 1000 # 家庭を訪問する回数
win = 0 # B君が勝った回数

for i in range(iter_max):
    pair = random.randrange(4) # ある家庭の2人の子供の(年上,年下)のペア: 0:(男,男),1:(男,女),2:(女,男),3:(女,女)
    if pair!=3:
        # 少なくとも1人が男と判明した場合、賭けが始まる
        bet_num += 1
        if pair!=0:
            # (男,女),(女,男)の場合は、B君の勝ち
            win += 1

print("賭けを行った回数: {}".format(bet_num))
print("B君の勝率: {}".format(win/bet_num))
print("B君が払った参加費の合計: {}".format(entry_fee*bet_num))
print("B君が貰った賞金の合計: {}".format(prize*win))
print("B君の収支: {}".format(prize*win - entry_fee*bet_num))
賭けを行った回数: 748
B君の勝率: 0.6644385026737968
B君が払った参加費の合計: 448800
B君が貰った賞金の合計: 497000
B君の収支: 48200

Qiskit Global Summer School 2023 Lab5 参加記録

各Labの記録

まえがき

Qiskit Global Summer School 2023 Lab5の記録です。

Lab 5: Error mitigation with Qiskit Runtime

  • エラー緩和(Error mitigation)をするラボ。
  • 具体的には、単純な観測値と初期状態を定義し、Estimatorプリミティブを使って期待値を測定する。
  • ノイズの多いシミュレーションを用いて、さまざまなエラー緩和戦略の効果を探る。

Setup

例として簡単なハイゼンベルグハミルトニアン・モデルを定義する。また、簡単な状態準備回路を構築する。

from qiskit import QuantumCircuit, QuantumRegister
from qiskit.quantum_info import SparsePauliOp


def heisenberg_hamiltonian(
    length: int, jx: float = 1.0, jy: float = 0.0, jz: float = 0.0
) -> SparsePauliOp:
    terms = []
    for i in range(length - 1):
        if jx:
            terms.append(("XX", [i, i + 1], jx))
        if jy:
            terms.append(("YY", [i, i + 1], jy))
        if jz:
            terms.append(("ZZ", [i, i + 1], jz))
    return SparsePauliOp.from_sparse_list(terms, num_qubits=length)


def state_prep_circuit(num_qubits: int, layers: int = 1) -> QuantumCircuit:
    qubits = QuantumRegister(num_qubits, name="q")
    circuit = QuantumCircuit(qubits)
    circuit.h(qubits)
    for _ in range(layers):
        for i in range(0, num_qubits - 1, 2):
            circuit.cx(qubits[i], qubits[i + 1])
        circuit.ry(0.1, qubits)
        for i in range(1, num_qubits - 1, 2):
            circuit.cx(qubits[i], qubits[i + 1])
        circuit.ry(0.1, qubits)
    return circuit
length = 5

hamiltonian = heisenberg_hamiltonian(length, 1.0, 1.0)
circuit = state_prep_circuit(length, layers=2)

print(hamiltonian)
circuit.draw("mpl")

ハミルトニアンと準備回路

Calculate exact expectation value (energy)(翻訳)

まず、Estimatorプリミティブのローカルシミュレータ実装を使って正確な期待値を計算する。ハミルトニアンの期待値は "エネルギー "とも呼ばれる。

from qiskit_aer.primitives import Estimator

estimator = Estimator(approximation=True)
job = estimator.run(circuit, hamiltonian, shots=None)
result = job.result()
exact_value = result.values[0]

print(f"Exact energy: {exact_value}")
# Exact energy: 4.290938711029713

Run noisy simulation through Qiskit Runtime

次に、Qiskit Runtimeサービスを初期化し、ノイズを扱えるシミュレータに支えられたEstimatorプリミティブの使用に切り替える。この回路は5量子ビットで動作するが、量子ビット選択の潜在的な効果を後で示すために、6量子ビットでシミュレータを初期化する。

from qiskit_ibm_runtime import QiskitRuntimeService

hub = "ibm-q-internal" # この辺の値は適宜自分用に変える
group = "deployed" #
project = "default" #

service = QiskitRuntimeService(instance=f"{hub}/{group}/{project}")
from qiskit_ibm_runtime import Estimator, Options, Session
from qiskit.transpiler import CouplingMap

backend = service.get_backend("simulator_statevector")
# set simulation options
simulator = {
    "basis_gates": ["id", "rz", "sx", "cx", "reset"],
    "coupling_map": list(CouplingMap.from_line(length + 1)),
}
shots = 10000

No noise

まずはノイズなしのシミュレータで実行してみる。

from qiskit_ibm_runtime import Estimator, Options, Session
from qiskit.transpiler import CouplingMap
import math

options = Options(
    simulator=simulator,
    resilience_level=0,
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 4.275399999999999
Energy error: 0.015538711029713603
Variance: 4.44106724
Standard error: 0.021073839801991474

誤差(Energy error)が 0.015... という値になった。ノイズなしでこのくらいになると覚えておいて、以下を進めていこう。

Readout error

次に、Readout error(読み出しエラー)ありのシミュレータで実行してみる。

ex1

この課題では、本当に悪い Readout error を持つ最初の量子ビットを除いたすべての量子ビットのReadout errorが緩やかなノイズモデルを構築する。

具体的には、以下の特性を持つノイズモデルを構築すること:

  • q0
    • 50%の確率で1を0と読み間違える
    • 20%の確率で0を1と読み間違える
  • それ以外の量子ビット
    • 5%の確率で1を0と読み間違える
    • 2%の確率で0を1と読み間違える

ノイズ作るライブラリは以下たちを参照すると良い。

from qiskit_aer.noise import NoiseModel, ReadoutError

noise_model = NoiseModel()

##### your code here #####
# その他の量子ビットに与えるReadoutError
rest_p0given1 = 0.05 # P(0|1)
rest_p1given0 = 0.02 # P(1|0)
# ReadoutErrorには [[P(0|0), P(1|0)], [P(0|1), P(1|1)]] のフォーマットで与える
rest_readout_error = ReadoutError([[1 - rest_p1given0, rest_p1given0], [rest_p0given1, 1 - rest_p0given1]])
noise_model.add_all_qubit_readout_error(rest_readout_error)

# q0に与えるReadoutError
q0_p0given1 = 0.5
q0_p1given0 = 0.2
q0_readout_error = ReadoutError([[1 - q0_p1given0, q0_p1given0], [q0_p0given1, 1 - q0_p0given1]])
noise_model.add_readout_error(q0_readout_error, (0,))

print(noise_model)
WARNING: Specific readout error on qubits (0,) overrides previously defined all-qubit readout error for these qubits.
NoiseModel:
  Basis gates: ['cx', 'id', 'rz', 'sx']
  Instructions with noise: ['measure']
  Qubits with noise: [0]
  All-qubits errors: ['measure']
  Specific qubit errors: [('measure', (0,))]

まず、読み出しエラーを緩和するために何もせずにシミュレーションを実行してみる。明示的にresilience_level = 0を設定し、Runtimeサービスによるエラー緩和が適用されないようにする。不適切な量子ビットの選択の影響を説明するために、量子ビット0を含む初期レイアウトを明示的に指定する。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0,
    transpilation=dict(initial_layout=list(range(length))),
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 3.5386
Energy error: 0.7523387110297128
Variance: 5.337599239999999
Standard error: 0.023103244880319302

この誤差(Energy error)はかなり大きい。これを改善するために、量子ビット0を避けるような量子ビットレイアウトを選んでみよう。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0,
    transpilation=dict(initial_layout=list(range(1, length + 1))), # initial_layoutからq0を取り除く
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 3.9506
Energy error: 0.3403387110297129
Variance: 4.962461319999999
Standard error: 0.02227658259248936

誤差(Energy error)は小さくなったが、それでもまだ大きい。resilience_level = 1に設定して、読み出しエラー緩和を有効にしてみよう。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=1, # 1に設定
    transpilation=dict(initial_layout=list(range(1, length + 1))),
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 4.303468384959293
Energy error: 0.012529673929580376
Variance: 7.169951196036428
Standard error: 0.02677676454696577

誤差(Energy error)はかなり小さくなり、読み出しエラーの影響はほぼ完全に緩和されたことがわかった。この緩和は「タダ」ではない。特に、

  • 読み出しエラーを緩和するために、ランタイムサービスは追加のキャリブレーション回路を実行する必要があるため、全体的な実行時間が長くなる可能性がある。
  • 推定値の分散が大きくなり、平均値の標準誤差が大きくなる。その結果、所定の標準誤差を達成するためには、より多くのショットを指定する必要がある。

通常このようなコストは比較的小さいため、読み出しエラーの緩和を可能にすることは、ほとんどの場合、価値がある。

ex2

読み出しエラーの緩和をオンにすると、推定値の分散(Variance)が2倍増加すると仮定する。 元々10,000ショットで実験を行った場合、平均の標準誤差(Standard error)を同じにするためには今何ショットを使用すべきか?

bellcurve.jp

元々 n=10,000ショットで実験を行った場合の標準偏差 \sigmaとする。

読み出しエラーの緩和をオンにしたときのショット数を n'標準偏差 \sigma'とする。

読み出しエラーの緩和をオンにすると分散が2倍になるので、標準偏差 \sqrt{2}倍になる。よって \sigma' = \sqrt{2} \sigma

標準誤差 SE = \frac{\sigma}{\sqrt{n}} より、

 \frac{\sigma}{\sqrt{10000}} = \frac{\sqrt{2}\sigma}{\sqrt{n'}} を解けばいい。 よって  n' = 20,000

new_shots: int

##### your code here #####
new_shots = 20000

Depolarizing error and zero-noise extrapolation

このセクションでは、ゼロノイズ外挿(zero-noise extrapolation)を使ってどのように脱分極エラー(depolarizing error)を軽減できるかを見ていく。Qiskit Runtimeのゼロノイズ外挿機能はまだベータ版であるため、現在いくつかの制限がある。

特に、本稿執筆時点では、ゼロノイズ外挿機能は読み出しエラーを緩和しない。

したがって以下の例では、ノイズモデルから読み出しエラーを取り除く。

ex3

各CNOTゲートの後に2量子ビットの脱分極エラーを加えるノイズモデルを構築し、エラーチャネルが1%の確率で入力量子状態を完全混合状態にマップするようにする。

qiskit.org

  • depolarizing_error
    • 第1引数:確率
    • 第2引数:n-qubit用の脱分極エラーかを指定する

今回CNOTゲートは2-qubitゲートなので、第2引数は2を指定する。

詳しい使い方はNoise Models (qiskit_aer.noise)の「Example: depolarizing noise model」に書いてある。

from qiskit_aer.noise import depolarizing_error

noise_model = NoiseModel()

##### your code here #####
depolar_p = 0.01
depolar_error = depolarizing_error(depolar_p, 2)

noise_model.add_all_qubit_quantum_error(depolar_error, ['cx']) # CNOTゲートに脱分極エラーを適用

print(noise_model)
NoiseModel:
  Basis gates: ['cx', 'id', 'rz', 'sx']
  Instructions with noise: ['cx']
  All-qubits errors: ['cx']

resilience_level = 0でestimatorを実行してみる。

# resilience_level = 0で実行
options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 4.0624
Energy error: 0.2285387110297128
Variance: 4.79633856
Standard error: 0.021900544650761543

次はresilience_level = 1で読み出しエラーの緩和をオンにしてestimatorを実行してみよう。今このノイズモデルには読み出しエラーの緩和は含まれていないため、これによる効果は期待できない。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=1, # 1に設定して読み出しエラーの緩和をオンにする
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 4.0586
Energy error: 0.23233871102971282
Variance: 4.8106246
Standard error: 0.02193313611866757

予想通り、resilience_level = 0のときと誤差(Energy error)は変わらない結果となった。

では、resilience_level = 2に設定して、ゼロノイズ外挿をオンにして実行してみよう。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # 2に設定してゼロノイズ外挿をオンにする
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 4.206533333333335
Energy error: 0.08440537769637846
Variances: [4.84564196, 5.46333092, 5.91612628]
Standard error: [0.022012818901721787, 0.02337376931519604, 0.024323088372984216]

誤差(Energy error)がかなり小さくなったことがわかる。脱分極ノイズの影響はほぼ完全に緩和されたおかげである。

定量に対する単一の分散値を得る代わりに、外挿のために測定された各データポイントに対する分散のリストを返すようになったことに注意。Qiskit Runtimeの将来のバージョンでは、これらの分散も外挿され、最終的な推定量に対して単一の分散が返されるようになるだろう。

ex4(ungraded)

脱分極エラー(depolarizing error)以外に、どのような種類のノイズがゼロノイズ外挿によって緩和できるだろうか?他のノイズモデルを構築し、ゼロノイズ外挿を使った場合と使わない場合でシミュレーションを行い、あなたの提案をテストしてみよ。

Noise Models (qiskit_aer.noise) の「Quantum Error Functions」の項目に色々エラーの種類が書いてある。

以下の記事を大いに参考にした。

Coherent Noiseをrx,ry,rzゲートに使ってみる。

# Coherent Noise
from qiskit_aer.noise import coherent_unitary_error
from qiskit.circuit.library import RYGate, RXGate, RZGate

noise_model = NoiseModel()

epsilon = math.pi/30 # over rotation amount

# RYゲートのコヒーレントノイズ
epsilon_rotation = RYGate(epsilon).to_matrix() # get matrix representation
over_rotation = coherent_unitary_error(epsilon_rotation)
# for i in range(5):
#     noise_model.add_quantum_error(over_rotation, ['ry'], qubits=[i])
noise_model.add_all_qubit_quantum_error(over_rotation, ['ry'])

# RXゲートのコヒーレントノイズ
epsilon_rotation = RXGate(epsilon).to_matrix()
over_rotation = coherent_unitary_error(epsilon_rotation)
noise_model.add_all_qubit_quantum_error(over_rotation, ['rx'])

# RZゲートのコヒーレントノイズ
epsilon_rotation = RZGate(epsilon).to_matrix()
over_rotation = coherent_unitary_error(epsilon_rotation)
noise_model.add_all_qubit_quantum_error(over_rotation, ['rz'])

print(noise_model)
NoiseModel:
  Basis gates: ['cx', 'id', 'rx', 'ry', 'rz', 'sx']
  Instructions with noise: ['ry', 'rx', 'rz']
  All-qubits errors: ['ry', 'rx', 'rz']

まずはresilience_level = 0で実行。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # resilience_level = 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 3.9776
Energy error: 0.3133387110297132
Variance: 4.8817328
Standard error: 0.02209464369479626

次にresilience_level = 2で実行。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # resilience_level = 2で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 3.6169500000000014
Energy error: 0.6739887110297116
Variances: [5.236187, 5.256725159999999, 5.26025432]
Standard error: [0.022882716184928747, 0.02292754928028724, 0.02293524431960558]

ゼロノイズ外挿によって緩和しなかった。というか誤差(Energy error)がひどくなってる。

pauli_errorをid,rx,ry,rzゲートに使ってみる

ある確率でbit-flipエラーやphase-flipエラーを起こすノイズを追加する。

from qiskit_aer.noise import pauli_error

# Construct a 1-qubit bit-flip and phase-flip errors
p_error = 0.01

# 確率p_errorでXゲートが適用される。確率1-p_errorでIゲートが適用される(何もしない)
bit_flip_error = pauli_error([('X', p_error), ('I', 1 - p_error)])
# 確率p_errorでZゲートが適用される。確率1-p_errorでIゲートが適用される(何もしない)
phase_flip_error = pauli_error([('Z', p_error), ('I', 1 - p_error)])

noise_model = NoiseModel()

noise_model.add_all_qubit_quantum_error(bit_flip_error, ['id', "rx", "ry", "rz"])
noise_model.add_all_qubit_quantum_error(phase_flip_error, ['id', "rx", "ry", "rz"])

print(noise_model)
WARNING: all-qubit error already exists for instruction "id", composing with additional error.
WARNING: all-qubit error already exists for instruction "rx", composing with additional error.
WARNING: all-qubit error already exists for instruction "ry", composing with additional error.
WARNING: all-qubit error already exists for instruction "rz", composing with additional error.
NoiseModel:
  Basis gates: ['cx', 'id', 'rx', 'ry', 'rz', 'sx']
  Instructions with noise: ['ry', 'id', 'rx', 'rz']
  All-qubits errors: ['id', 'rx', 'ry', 'rz']

まずはresilience_level = 0で実行。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # resilience_level = 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 2.5622
Energy error: 1.7287387110297132
Variance: 6.7096326
Standard error: 0.025902958518285127

次にresilience_level = 2で実行。

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # resilience_level = 2で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 2.467666666666668
Energy error: 1.823272044363045
Variances: [6.82206672, 6.840397360000001, 6.86765392]
Standard error: [0.02611908635461815, 0.026154153322178107, 0.026206209035264907]

ゼロノイズ外挿によって緩和しなかった。

thermal_relaxation_error(熱緩和エラー)をid,rx,ry,rzゲートに使ってみる

from qiskit_aer.noise import thermal_relaxation_error

t_error = thermal_relaxation_error(161.83e3, 61.64e3, 368)

noise_model = NoiseModel()

noise_model.add_all_qubit_quantum_error(t_error, ['id', "rx", "ry", "rz"])

print(noise_model)

resilience_level = 0で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # resilience_level = 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 3.908
Energy error: 0.3829387110297131
Variance: 5.06420016
Standard error: 0.022503777816180112

resilience_level = 2で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # resilience_level = 2で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 3.8287833333333356
Energy error: 0.4621553776963774
Variances: [5.20887776, 5.223925960000001, 5.2486150799999995]
Standard error: [0.022822965977278238, 0.022855909432792212, 0.022909856132241425]

エラー緩和しなかった。

thermal_relaxation_error(熱緩和エラー)をCNOTゲートに使ってみる

今までrxゲートなどにノイズをのせていたが、よくよく考えたらLabではCNOTゲートにノイズを乗せていたのであった。なのでそうしてみる。

import numpy as np

# 5量子ビット系
qn = 5

# T1 and T2 values for qubits 0-1
T1s = np.random.normal(50e3, 10e3, qn) # Sampled from normal distribution mean 50 microsec
T2s = np.random.normal(70e3, 10e3, qn)  # Sampled from normal distribution mean 50 microsec

# Truncate random T2s <= T1s
T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(qn)])

# Instruction times (in nanoseconds)
time_u1 = 0   # virtual gate
time_u2 = 50  # (single X90 pulse)
time_u3 = 100 # (two X90 pulses)
time_cx = 300
time_reset = 1000  # 1 microsecond
time_measure = 1000 # 1 microsecond

# QuantumError objects
errors_reset = [thermal_relaxation_error(t1, t2, time_reset)
                for t1, t2 in zip(T1s, T2s)]
errors_measure = [thermal_relaxation_error(t1, t2, time_measure)
                  for t1, t2 in zip(T1s, T2s)]
errors_u1  = [thermal_relaxation_error(t1, t2, time_u1)
              for t1, t2 in zip(T1s, T2s)]
errors_u2  = [thermal_relaxation_error(t1, t2, time_u2)
              for t1, t2 in zip(T1s, T2s)]
errors_u3  = [thermal_relaxation_error(t1, t2, time_u3)
              for t1, t2 in zip(T1s, T2s)]
errors_cx = [[thermal_relaxation_error(t1a, t2a, time_cx).expand(
             thermal_relaxation_error(t1b, t2b, time_cx))
              for t1a, t2a in zip(T1s, T2s)]
               for t1b, t2b in zip(T1s, T2s)]
print(errors_cx)
[[QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c674160>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c91ad10>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c91bf70>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c675360>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c91b130>, 1.0)])], [QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7827d0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c780ac0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c783b50>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c781ab0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7819c0>, 1.0)])], [QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c91ace0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c675e70>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c783d00>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7830d0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7808b0>, 1.0)])], [QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c783dc0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c783550>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c781e40>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7832b0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c783310>, 1.0)])], [QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c782ce0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c780f10>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c780fa0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff68c7801c0>, 1.0)]), QuantumError([(<qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7ff6a8cae230>, 1.0)])]]
noise_model = NoiseModel()

for j in range(4):
    # noise_thermal.add_quantum_error(errors_reset[j], "reset", [j])
    # noise_thermal.add_quantum_error(errors_measure[j], "measure", [j])
    # noise_thermal.add_quantum_error(errors_u1[j], "u1", [j])
    # noise_thermal.add_quantum_error(errors_u2[j], "u2", [j])
    # noise_thermal.add_quantum_error(errors_u3[j], "u3", [j])
    for k in range(4):
        noise_model.add_quantum_error(errors_cx[j][k], "cx", [j, k])

resilience_level = 0で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # resilience_level = 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")

resilience_level = 2で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # resilience_level = 2で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 4.286916666666666
Energy error: 0.004022044363047428
Variances: [4.57962836, 4.97150324, 5.158112879999999]
Standard error: [0.02140006626157966, 0.022296868031183213, 0.022711479212063665]

エラー緩和した。ゼロノイズ外挿はCNOTゲートにしか効かない……?知識不足でわからない。

Coherent NoiseをCNOTゲートに使ってみる

# Coherent Noise
from qiskit_aer.noise import coherent_unitary_error
from qiskit.circuit.library import RYGate, RXGate, RZGate, CUGate

noise_model = NoiseModel()

epsilon = math.pi/30 # over rotation amount

# CUゲートのコヒーレンスノイズ
epsilon_rotation = CUGate(epsilon, epsilon, epsilon, epsilon).to_matrix() # get matrix representation
over_rotation = coherent_unitary_error(epsilon_rotation)
noise_model.add_all_qubit_quantum_error(over_rotation, ['cx'])
        
print(noise_model)
NoiseModel:
  Basis gates: ['cx', 'id', 'rz', 'sx']
  Instructions with noise: ['cx']
  All-qubits errors: ['cx']

resilience_level = 0で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=0, # resilience_level = 0で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variance = result.metadata[0]["variance"]
std = math.sqrt(variance / shots)

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variance: {variance}")
print(f"Standard error: {std}")
Estimated energy: 3.9168000000000003
Energy error: 0.37413871102971274
Variance: 5.076300159999999
Standard error: 0.02253064615140897

resilience_level = 2で実行

options = Options(
    simulator=dict(noise_model=noise_model, **simulator),
    resilience_level=2, # resilience_level = 2で実行
)

with Session(service=service, backend=backend):
    estimator = Estimator(options=options)
    job = estimator.run(circuit, hamiltonian, shots=shots)

result = job.result()
experiment_value = result.values[0]
error = abs(experiment_value - exact_value)
variances = result.metadata[0]["zne"]["noise_amplification"]["variance"]
std = [math.sqrt(var/shots) for var in variances]

print(f"Estimated energy: {experiment_value}")
print(f"Energy error: {error}")
print(f"Variances: {variances}")
print(f"Standard error: {std}")
Estimated energy: 4.7794000000000025
Energy error: 0.4884612889702895
Variances: [5.0520554, 6.75144548, 7.730010119999999]
Standard error: [0.022476777793981058, 0.025983543792177387, 0.02780289574846476]

エラー緩和しなかった。これはCUGate使ってるのが間違ってそう。でもうーん、わからん。

Qiskit Global Summer School 2023 Lab4 参加記録

各Labの記録

まえがき

Qiskit Global Summer School 2023 Lab4の記録です。

Lab 4 : Iterative phase estimation

背景(翻訳)

量子位相推定(QPE)アルゴリズムは、最も重要で有名な量子アルゴリズムの1つである。Shorの因数分解アルゴリズムの重要なサブルーチンであり、量子シミュレーションのアルゴリズムでもある。このアルゴリズムの教科書的なバージョンは、必要な精度に応じてスケーリングされる補助的な量子ビットの数を使用するため、量子ビットの数や接続性が制限された今日のノイズの多いデバイスで実行するのは困難な回路になる。

反復位相推定(IPE)は、補助量子ビットを1つだけ必要とするQPEの変種である。IPEでは、補助量子ビットが繰り返し測定され、その測定結果が将来の量子演算の指針となる。このような古典的なフィードフォワードは、以前はIBMの量子プロセッサーで実行することは不可能だったが、最近導入された動的回路機能を使えば可能になる。

他の位相推定アルゴリズムと同様に、IPEは以下の問題を解決するように設計されている:

問題文: ユニタリー行列Uと、未知の固有値 e^{i2πφ}を持つUの固有状態|Ψ⟩が与えられたとき、φの値を推定せよ。

この問題では、いくつかの重要な点を明確にする必要がある。つまりUと|Ψ⟩がどのように指定されるかである。

UはUを実装する量子回路として与えられ、実際、我々は正の整数tに対してcontrolled- U^{2^{t}} の演算を効率的に実装する能力を持っていると仮定する。固有状態も量子回路として与えられる。

簡単のために、まずφが厳密な二項展開ができると仮定してみよう。

φの二項展開

ここで、最後の等式では基数2の「小数点」表記を用いている。簡単のために、Uが1量子ビットに作用するユニタリー演算子であるとする(ここで述べることは、Uが複数の量子ビットに作用する場合にも当てはまる)。IPEは補助量子ビットを必要とするので、q0とq1の2つの量子ビットの系が必要になる。q0は補助量子ビットで、q1はUが作用する物理系を表す。

ここで、q0を状態|+⟩に、q1を状態|Ψ⟩に初期化したとする。q0をコントロール、q1をターゲットとしてcontrolled- U^{2^{t}}ゲートを適用するとどうなるか?|Ψ⟩は固有値 e^{i2πφ}を持つUの固有状態なので、次のようになる。

こうなる。

つまり、システム量子ビットの状態は変化せず、 e^{i2π2^{𝑡}φ}の位相が補助量子ビットの状態に「キックバック」されたことになる。

さて、次のことに注意してほしい。

小数点表記に書ける

最後の等式では、任意の整数nに対して e^{i2πn}=1であるため、位相の「10進数」表現の整数部分は消えている。たとえば、

  • t=0の場合、位相は e^{i2\pi2^{0}\varphi} = e^{i2\pi\varphi} = e^{i2\pi0.\varphi_1 \varphi_2 ... \varphi_m}

  • t=1の場合、位相は e^{i2\pi2^{1}\varphi} = e^{i2\pi\varphi_1} e^{i2\pi0.\varphi_2 \varphi_3 ... \varphi_m} = e^{i2\pi0.\varphi_2 \varphi_3 ... \varphi_m}

  • t=2の場合、位相は e^{i2\pi2^{2}\varphi} = e^{i2\pi0.\varphi_3 \varphi_4 ... \varphi_m}

  • t=m-1の場合、位相は  e^{i2\pi2^{m-1}\varphi} = e^{i2\pi0.\varphi_m}

最後のt=m-1の場合、位相は e^{i2\pi0.\varphi_m}となり、 \varphi_m = 0の場合は1、 \varphi_m = 1の場合は-1となる。最初のケースでは、補助量子ビットq0は|+⟩状態にあり、2番目のケースでは|-⟩状態にある。したがって、Pauli X基底で量子ビットを測定すれば、最初のケースでは0を測定し、2番目のケースでは1を測定するようになり、100%の成功率でこれらのケースを区別できる。測定されたビットは \varphi_mに等しくなる。

Pauli X基底での測定は、量子ビットを測定する前にアダマールゲートを実行することで行われる。

実装

固有値 e^{i\pi/2}=e^{i2\pi・1/4}を持つ固有状態|Ψ⟩=|1⟩を使う。つまり𝜑=1/4=0.01=0.𝜑1𝜑2となる。

𝜑は2ビットで正確に表現できるので、量子回路の実装では2ビットの古典レジスタを使う。

q0が補助量子ビット、q1が|Ψ⟩。

ex1

UゲートとしてSゲートを使って反復位相推定を行う。

Sゲート

controlled-SゲートはCPhaseGateを使えば良い。

CPhaseGateの行列

from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
import numpy as np


def step_1_circuit(qr: QuantumRegister, cr: ClassicalRegister) -> QuantumCircuit:
    # qr is a quantum register with 2 qubits
    # cr is a classical register with 2 bits

    qc = QuantumCircuit(qr, cr)

    ####### your code goes here #######
    qc.h(0)
    qc.x(1)
    qc.cp(np.pi, 0, 1)
    qc.h(0)
    qc.measure(0,0)

    return qc


qr = QuantumRegister(2, "q")
cr = ClassicalRegister(2, "c")
qc = QuantumCircuit(qr, cr)
qc = step_1_circuit(qr, cr)
qc.draw("mpl")

qc回路

ex2

def step_2_circuit(qr: QuantumRegister, cr: ClassicalRegister) -> QuantumCircuit:
    # qr is a quantum register with 2 qubits
    # cr is a classical register with 2 bits

    # begin with the circuit from Step 1
    qc = step_1_circuit(qr, cr)

    ####### your code goes here #######
    qc.reset(qr[0])
    qc.h(qr[0])
    with qc.if_test((cr[0], 0)) as else_:
        # 0を測定したら、正しいので何もしない
        pass
    with else_:
        # 1を測定したら、-pi/2で位相補正する
        qc.rz(-np.pi/2, qr[0])

    return qc


qr = QuantumRegister(2, "q")
cr = ClassicalRegister(2, "c")
qc = QuantumCircuit(qr, cr)
qc = step_2_circuit(qr, cr)
qc.draw("mpl")

ex2の回路

さて、この回路をAerSimulatorで実行してみる。

from qiskit_aer import AerSimulator

sim = AerSimulator()
job = sim.run(qc, shots=1000)
result = job.result()
counts = result.get_counts()
counts
# {'01': 1000}

結果は 01 となり、𝜑=1/4=0.01=0.𝜑1𝜑2 を正しく推定できる回路を作ることができた。

ex3

Tゲートの位相を推定するIPEを作れという問題。

これ以降は以前の問題と同じと気づいたので詳細は省略。この回路は何ビット必要か?固有状態は?など、詳しくはIBM Quantum Challenge: Spring 2023 参加記録のLab3のExercise 3を参照。

TゲートのIPE回路は以下になる。

from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
import numpy as np


def t_gate_ipe_circuit(qr: QuantumRegister, cr: ClassicalRegister) -> QuantumCircuit:
    # qr is a quantum register with 2 qubits
    # cr is a classical register with 3 bits

    qc = QuantumCircuit(qr, cr)

    ####### your code goes here #######
    """𝜑1のステップ"""
    # 初期化
    qc.h(qr[0])
    qc.x(qr[1])
    
    # 𝜑1の推定
    qc.cp(np.pi, qr[0], qr[1])
    qc.h(qr[0])
    qc.measure(qr[0], cr[0])
    
    """𝜑2のステップ"""
    # 初期化
    qc.reset(qr[0])
    qc.h(qr[0])
    
    # Z軸周りの角度wの回転で位相補正
    with qc.if_test((cr[0], 1)):
        qc.p(-2*np.pi/(2**2), qr[0])
    
    # 𝜑2の推定
    qc.cp(np.pi/(2**1), qr[0], qr[1])
    qc.h(qr[0])
    qc.measure(qr[0], cr[1])
    
    """𝜑3のステップ"""
    # 初期化
    qc.reset(qr[0])
    qc.h(qr[0])     

    # Z軸周りの角度wの回転で位相補正
    with qc.if_test((cr[0], 1)):
        qc.p(-2*np.pi/(2**3), qr[0])
    with qc.if_test((cr[1], 1)):
        qc.p(-2*np.pi/(2**2), qr[0])

    # 𝜑3の推定
    qc.cp(np.pi/(2**2), qr[0], qr[1])
    qc.h(qr[0])
    qc.measure(qr[0], cr[2])   
    
    return qc


qr = QuantumRegister(2, "q")
cr = ClassicalRegister(3, "c")
qc = QuantumCircuit(qr, cr)
qc = t_gate_ipe_circuit(qr, cr)
qc.draw("mpl")

TゲートのIPE回路

from qiskit_aer import AerSimulator

sim = AerSimulator()
job = sim.run(qc, shots=1000)
result = job.result()
counts = result.get_counts()
counts
# {'001': 1000}

ex4

詳しくはIBM Quantum Challenge: Spring 2023 参加記録のLab3のExercise 4を参照。

ex5

詳しくはIBM Quantum Challenge: Spring 2023 参加記録のLab3のExercise 5を参照。

ex6(ungraded)

ここまでに構築したIPE回路は、特定のゲートと特定の精度ビット数用に設計されたものであった。ここで一般化して、異なるゲートと精度レベルを扱える一般的なIPEルーチンを実装してみよう。

次の関数を完成させて、一般化された IPE ルーチンを実装しよう。以下の入力を受け取る:

from qiskit.circuit import Gate


def iterative_phase_estimation(
    qr: QuantumRegister,
    cr: ClassicalRegister,
    controlled_unitaries: list[Gate],
    state_prep: Gate,
) -> QuantumCircuit:
    qc = QuantumCircuit(qr, cr)

    ####### your code goes here #######
    for i, u_gate in enumerate(controlled_unitaries[::-1]):
        # リセット処理
        if i != 0:
            qc.reset(qr[0])
            
        qc.h(qr[0])
        if i == 0:
            qc.x(qr[1])
        else:
            # Z軸周りの角度wの回転で位相補正
            for j in range(0, i):
                with qc.if_test((cr[j], 1)):
                    qc.p(-np.pi /(2**(i-j)), qr[0])

        qc.append(u_gate, [qr[0], qr[1]])

        # X基底でq0を測定
        qc.h(qr[0])
        qc.measure(qr[0], cr[i])

    return qc

この関数を使用してS ゲートのIPE回路を生成する。

from qiskit.circuit.library import CPhaseGate, XGate

qr = QuantumRegister(2, "q")
cr = ClassicalRegister(2, "c")

s_angle = np.pi / 2
controlled_unitaries = [CPhaseGate(s_angle * 2**k) for k in range(2)]
qc = iterative_phase_estimation(qr, cr, controlled_unitaries, XGate())
qc.draw()

iterative_phase_estimationにSゲートを与えたときの回路

回路を実行する。結果は01になるはずである。

sim = AerSimulator()
job = sim.run(qc, shots=1000)
result = job.result()
counts = result.get_counts()
counts
# {'01': 1000}

合ってる。

Tゲートでも検証してみる。

qr = QuantumRegister(2, "q")
cr = ClassicalRegister(3, "c")

t_angle = np.pi / 4
controlled_unitaries = [CPhaseGate(t_angle * 2**k) for k in range(3)]
qc = iterative_phase_estimation(qr, cr, controlled_unitaries, XGate())
qc.draw()

iterative_phase_estimationにTゲートを与えたときの回路

回路を実行する。結果は001になるはずである。

sim = AerSimulator()
job = sim.run(qc, shots=1000)
result = job.result()
counts = result.get_counts()
counts
# {'001': 1000}

合ってる。