Matrix Product Stateによる量子ゲート計算をJuliaでやってみる

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

IBM-Q の53qubit量子コンピュータのシミュレータを作ってみた
という記事で紹介されていた、Matrix Product State(MPS)という状態の記述方式を利用して量子状態を必要最小限のパラメータ数で表す手法について。

非常に勉強になる記事ですが自分の場合はある程度手を動かさないと理解が追いつかなそうだったため、少し実装してみました。

手法の解説については同記事でかなり丁寧にされていると思いますので、本記事は自身の学習進捗記録を兼ねて「やってみた」のノリでいきます。
不正確な点があったらすみません。

実装は上記事で紹介されていた論文を参考にしました。
http://arxiv.org/abs/1008.3477
言語くらいはまだ書かれてなさそうなものを、ということでどこかで使ってみたかったJuliaを試しています。

モジュールの用意

今回使う関数(+α)をJuliaのmoduleにまとめました。

module MPSforQuantum
    export MPS, restore, OneQubitGate, dstack, ddstack, mps_size
    export CX, inner_product
    using LinearAlgebra

    function SVD(C::Array{ComplexF64,2}, eps::Float64)
        #eps = 1e-2
        A = svd(C)
        filter!((x) -> x > eps, A.S)
        l = length(A.S)
        return A.U[:, 1:l], diagm(A.S) * A.Vt[1:l, :]
    end

    function SVD(C::Array{ComplexF64,2}, D::Int64)
        A = svd(C)
        if length(A.S) > D
            return A.U[:, 1:D], diagm(A.S[1:D]) * A.Vt[1:D, :]
        else
            return A.U[:, :], diagm(A.S) * A.Vt[:, :]
        end
    end

    function MPS(C::Array{ComplexF64,1}, param)
        # C = copy(C0)
        d = 2
        N = Int64(log2(size(C)[1]))
        arrs = []
        #eps = 1e-2
        #D = 550

        r = 1
        for i = 1:N-1
            C = reshape(C, d * r, d^(N - i))
            tmp_A, C = SVD(C, param)
            r = size(tmp_A, 2)
            col = convert(Int64, size(tmp_A, 1) / 2)
            push!(arrs, [tmp_A[1:col, :], tmp_A[(col+1):col*2, :]])
        end

        Ct = transpose(C)
        col2 = convert(Int64, size(C, 2) / 2)
        push!(arrs, [C[:, 1:col2], C[:, (col2+1):col2*2]])
        return arrs
    end

    function mps_size(arrs::Array{Any, 1})
        params = 0
        N = size(arrs)[1]
        for i = 1:N
            println("array $i's size: ", size(arrs[i][1]))
            params += length(arrs[i][1]) * 2
        end
        println("Num of parameters: ", params)
        println("2^N: ", 2^N)
    end

    # Restore prob. ampl. of a state from MPS
    function restore(arrs::Array{Any, 1}, n::Int64) # n is decimal representation of the state '0110...'
        N = Int64(size(arrs)[1])
        s = bitstring(n)[end - (N - 1):end] # 後ろからN bit目まで
        phys_idx = [convert(Int64, s[i]) - 48 for i = length(s):-1:1] # 後ろが最上位ビット
        return prod([arrs[i][phys_idx[i] + 1] for i = 1:length(phys_idx)])[1]
    end

    function OneQubitGate(arrs::Array{Any, 1}, O::Array{Complex{Float64},2}, n::Int64)
        arrs_ = similar(arrs[n + 1])
        arrs__ = copy(arrs)
        arrs_[1] = arrs[n + 1][1] * O[1, 1] + arrs[n + 1][2] * O[2, 1]
        arrs_[2] = arrs[n + 1][1] * O[1, 2] + arrs[n + 1][2] * O[2, 2]
        arrs__[n + 1] = arrs_
        return arrs__
    end

    function dstack(A::Array{ComplexF64,2}, B::Array{ComplexF64,2})
        return cat(A, B, dims = 3)
    end
    
    function ddstack(A::Array{ComplexF64,2}, B::Array{ComplexF64,2}, 
                    C::Array{ComplexF64,2}, D::Array{ComplexF64,2})
        AC = cat(A, C, dims = 3)
        BD = cat(B, D, dims = 3)
        return cat(AC, BD, dims = 4)
    end

    function inner_product(arrs1::Array{Any, 1}, arrs2::Array{Any, 1})
        N = Int64(size(arrs1)[1])
        ip = arrs1[1][1]' * arrs2[1][1] + arrs1[1][2]' * arrs2[1][2]
        for i = 2:N
            phys0 = arrs1[i][1]' * ip * arrs2[i][1]
            phys1 = arrs1[i][2]' * ip * arrs2[i][2]
            ip = phys0 + phys1
        end
        return ip
    end

    function CX(arrs::Array{Any, 1}, n::Int64, t::Number)
        tmp = copy(arrs)
        arrs_ = [tmp[n + 1][1]*tmp[n + 2][1] tmp[n + 1][2]*tmp[n + 2][2]
                tmp[n + 1][1]*tmp[n + 2][2] tmp[n + 1][2]*tmp[n + 2][1]]
        arrs__ = MPSforQuantum.SVD(arrs_, t)
        col = convert(Int64, size(arrs__[1], 1) / 2)
        tmp[n + 1][1] = arrs__[1][1:col, :]
        tmp[n + 1][2] = arrs__[1][(col+1):(2*col), :]
        tmp[n + 2][1] = arrs__[2][:, 1:col]
        tmp[n + 2][2] = arrs__[2][:, (col+1):(2*col)]
        return tmp
    end
end

計算

まずは上のモジュールを"using ~"で使えるようにします。
"MPS.jl"はモジュールを記述しているファイル名です。

using LinearAlgebra
include("./MPS.jl")
using .MPSforQuantum

まずは簡単な状態として$\ket{000}$をMPSに変換してみます。

実行コード

N = 20
C0 = zeros(ComplexF64, 2^N)
C0[1] = 1 # '000'

#D = 10
eps = 1e-3
mps = MPS(C0, eps) # convert to MPS
mps_size(mps) # print the size of MPS

n = 0 # check the state |bin(n)>
println("Restore from MPS: ", restore(mps, n))
println("Original state: ", C0[n + 1])

結果

array 1's size: (1, 1)
array 2's size: (1, 1)
array 3's size: (1, 1)
array 4's size: (1, 1)
array 5's size: (1, 1)
array 6's size: (1, 1)
array 7's size: (1, 1)
array 8's size: (1, 1)
array 9's size: (1, 1)
array 10's size: (1, 1)
array 11's size: (1, 1)
array 12's size: (1, 1)
array 13's size: (1, 1)
array 14's size: (1, 1)
array 15's size: (1, 1)
array 16's size: (1, 1)
array 17's size: (1, 1)
array 18's size: (1, 1)
array 19's size: (1, 1)
array 20's size: (1, 1)
Num of parameters: 40
2^N: 1048576
Restore from MPS: 1.0 + 0.0im
Original state: 1.0 + 0.0im

20 qubitからなる量子状態は20個の行列の積に分解され、各行列のサイズを表示しています(ここでは1x1なのでスカラーと変わらないですが)。
実際は各qubitが|0>, |1>の2状態をとりうるため、$20 \times 2$ 個の行列が保持されています。

変数epsについて、MPSへの変換時に一定よりweightの小さい次元を削除してパラメータ数を圧縮しているのですが、その閾値です。
圧縮した結果、$2^{20} = 1048576$ 次元ヒルベルト空間に張られた状態がたった40個のパラメータで表わせました。

最下部二行が、実行コードで指定(n = 0 ⇔ $\ket{00...0}$)したある1状態の係数です。
ここが全てのn (0 ~ 2^20 - 1)について一致していればMPSは元の量子状態の係数を保持できています。

次は試しに、ランダムな係数をもつ20 qubit量子状態をMPSに変換してみます。

実行コード

N = 20
C0 = normalize!(rand(ComplexF64, 2^N)) # Changed!

eps = 1e-3
mps = MPS(C0, eps) # convert to MPS
mps_size(mps) # print the size of MPS

n = 1000 # check the state |bin(n)>
println("Restore from MPS: ", restore(mps, n))
println("Original state: ", C0[n + 1])

出力

array 1's size: (1, 2)
array 2's size: (2, 4)
array 3's size: (4, 8)
array 4's size: (8, 16)
array 5's size: (16, 32)
array 6's size: (32, 64)
array 7's size: (64, 128)
array 8's size: (128, 256)
array 9's size: (256, 512)
array 10's size: (512, 982)
array 11's size: (982, 512)
array 12's size: (512, 256)
array 13's size: (256, 128)
array 14's size: (128, 64)
array 15's size: (64, 32)
array 16's size: (32, 16)
array 17's size: (16, 8)
array 18's size: (8, 4)
array 19's size: (4, 2)
array 20's size: (2, 1)
Num of parameters: 2710184
2^N: 1048576
Restore from MPS: 0.0003481771015525017 + 2.4616598541730162e-5im
Original state: 0.00035035163995938485 + 2.5586690026457734e-5im

MPSのパラメータ数が2710184個もあります(下から4行目)。
20 qubitなので本来パラメータは$2^{20} = 1048576$個もあれば十分なはず。

これは変換元の量子状態の係数を完全ランダムとしているためです。
このように生成された量子状態は一般に各qubitが複雑にエンタングルしており、MPSに変換時の圧縮がほぼ効ききません。

次に、HadamardゲートとCXゲートを使ってBell状態を作ってみます。

実行コード

Hadamard = convert(Array{ComplexF64,2}, [1 1; 1 -1] ./ sqrt(2))

N = 2
C0 = zeros(ComplexF64, 2^N)
C0[1] = 1 # '000'
eps = 1e-3

mps1 = MPS(C0, eps)
println("Initial MPS")
mps1 = OneQubitGate(mps1, Hadamard, 0) # -> '001' & '000'
mps_size(mps1)

arr = CX(mps1, 0, eps)
println("\nBell state MPS")
mps_size(arr)
println()
for i=0:3
    res = restore(arr, i)
    println("|", bitstring(i)[end - 1:end], ">: ",  res)
end

出力

Initial MPS # 追記：|00>です
array 1's size: (1, 1)
array 2's size: (1, 1)
Num of parameters: 4
2^N: 4

Bell state MPS
array 1's size: (1, 2)
array 2's size: (2, 1)
Num of parameters: 8
2^N: 4

|00>: 0.7071067811865475 + 0.0im
|01>: 0.0 + 0.0im
|10>: 0.0 + 0.0im
|11>: 0.7071067811865475 + 0.0im

MPSをBell状態へと操作することでエンタングルメントが付与された結果、MPSのパラメータ数が増えました(4→8)。
(仮にBell状態を作ってからMPSに変換しても同じ結果になります。)
原則として、量子状態をMPSで(一定の精度で)表すのに必要なパラメータは状態のエンタングルメント量によって決まるようです。

まとめ

以上の例は極端ですが、qubit間のエンタングルメントが小さい状態ほどMPSに圧縮した際にパラメータ数が削減できていることを確認できたと思います。

もう少し言うと、エンタングルメントが小さいqubitsは $2^N$次元ヒルベルト空間のほんの一部分しか占めず、じゃあそれを実際に必要な最小限のパラメータで表そう、という話になります。
そうすることで、古典計算機でより大きな量子系をシミュレートできる可能性が出てきます。

長くなってしまったため今回はここまでです。
最後までありがとうございました。

参考にしたもの

• IBM-Q の53qubit量子コンピュータのシミュレータを作ってみた
• http://arxiv.org/abs/1008.3477