PennyLane 機能色々(TIPS)

自分用のメモです。TIPSとかチートシートとかいう。
2021/10/25 更新

What is pennylane?

いつものimport

import pennylane as qml
from pennylane import numpy as np
from matplotlib import pyplot as plt

回路関係

基本構文

dev = qml.device("default.qubit", wires=n_qubits)
@qml.qnode(dev)
def circuit(theta):
    qml.RX(theta, wires=range(n_qubits))
    return qml.expval(qml.PauliZ(i))

devとして、シミュレータのタイプや実機のタイプを指定する。
wires は量子ビット数のことと思えば良い。

@qml.qnode(dev)　は　量子回路circuit にdevをくっつける記法。(Pythonのデコレーター機能)
この下に　def された関数（ここではcircuit)は、devと紐づく。

デコレートされたcircuitは qnode というクラスになっている。

type(circuit)
>pennylane.qnode.QNode

circuitを実行すると、指定したデバイスで回路が実行された後の結果が返ってくる。

デコレータを使わずに、以下のように記述しても同じ。

def circuit_alone(theta):
    qml.RX(theta, wires=range(n_qubits))
    return qml.expval(qml.PauliZ(i))

My_Qnode = qml.QNode(circuit_alone,dev)

type(My_Qnode)
>pennylane.qnode.QNode

circuitで使えるゲートはqiskit等とほぼ同じ。wiresとして量子ビット位置を指定する。

回路の戻り値はやや特殊で、基本的には期待値 expval や　サンプル sample や　確率 prob を返す必要がある。
（明示的に戻り値を指定しないとバグる）

状態ベクトル state も使えるよう。

return qml.state()

dev._stateでも出来るらしい？
https://discuss.pennylane.ai/t/pennylane-statevector-results-different-from-qiskit-statevector-issue/1072/28

なおオブザーバブルとしてパウリ行列積の線形和（ハミルトニアン）を使いたい場合は、Hamiltonianを構築してから ExpValCost を使う。
https://pennylane.readthedocs.io/en/latest/code/api/pennylane.ExpvalCost.html
ExpValCostの引数の ansatz は（デコレータの無い)circuitの意味。

##使えるdev

default.qubit, lightning.qubit, qulacs.simulator等。
シミュレータの処理速度としてはlightning.qubitとqulacs.simulatorが現状速い。

##使える勾配計算手法
parameter-shift, backprop, adjoint等
adjointが非常に速いが、シミュレータでしか使えない。
実機で使えるのはparameter-shiftのみ。
devでqulacsを指定するとparameter-shiftになってしまうので、lightning.qubitを指定してadjointを使うのが良い。

どうしてもbackpropを使いたい場合は、interfaceをちゃんと設定する必要がある。
devは.tf を使わないといけないし、デコレーターの定義で interface="tf"が必須。

##さらに早く
変分回路等で回路構造が変わらない場合は　mutable = False に設定すると、都度回路構築をしなくなるのでちょっと早い。

テンプレート回路

Amplitude/Angle embedding, StronglyEntanglingLayers が頻出

期待値

テンソル積の期待値を取りたいとき(量子ビット数が多いのでfor文で）
https://discuss.pennylane.ai/t/how-to-use-a-tensor-measurements-with-variable-number-of-qubits/1287/2

return qml.expval(Tensor(*[qml.PauliZ(i) for i in range(n_qubits)]))

or

obs = reduce(operator.matmul, [qml.PauliZ(i) for i in range(n_qubits)])
return qml.expval(obs)

明示的に書き並べたい場合は、テンソル積記号として @ が使える。

テンソル積ではなく、各1ビットへの作用(例えばIIZII)が取りたい時は以下のようにする。

@qml.qnode(dev)
def qnode(inputs, weights):
    qml.templates.AngleEmbedding(inputs, wires=range(n_qubits))
    qml.templates.StronglyEntanglingLayers(weights, wires=range(n_qubits))
    return [qml.expval(qml.PauliZ(i)) for i in range(n_qubits)]

回路描画

pennylane-qiskit使用

#!pip install pennylane-qiskit pylatexenc
import pennylane as qml
from pennylane import numpy as np

dev = qml.device("qiskit.aer", wires=1)

@qml.qnode(dev)
def f(x):
    qml.RX(x, wires=0)
    return qml.expval(qml.PauliZ(0))

f(0.3)

dev._circuit.draw('mpl')

sampleを集計

def sample_to_counts(sample):
    n_qubits, n_shots = np.shape(sample)
    sample_bin = (1+sample)/2
    weights = [2**i for i in range(n_qubits)]
    sample_dec = np.average(sample_bin,axis=0,weights=weights)*np.sum(weights) # sum b_n*2^n
    sample_dec = sample_dec.astype('int8')
    u, freq = np.unique(sample_dec, return_counts=True) # count frequency
    counts=dict()
    for i,j in zip(u,freq):
        counts.setdefault(format(i, '0'+str(n_qubits)+'b'),j)
    for i in range (2**n_qubits):
        counts.setdefault(format(i, '0'+str(n_qubits)+'b'),0) # zero-padding
    return counts

n_qubits = 3
dev = qml.device("default.qubit", wires=n_qubits, shots=100)

@qml.qnode(dev)
def circuit():
    qml.Hadamard(wires=[0])
    for i in range(0,n_qubits-1):
        qml.CNOT(wires=[i,i+1])
    return [qml.sample(qml.PauliZ(i)) for i in range(n_qubits)]

counts = sample_to_counts(circuit())

import matplotlib.pyplot as plt
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts');
plt.xticks(rotation=90);

アプリケーション別

QCL

w/o Tensorflow.keras

Gradient Descent

import pennylane as qml
from pennylane import numpy as np
from matplotlib import pyplot as plt
num_of_data = 256
X =  np.random.uniform(high=2 * np.pi, size=(num_of_data,1))
Y = np.sin(X[:,0])

########  parameters #############
n_qubits = 2 ## num. of qubit
n_layers = 2 # num of q_layers

dev = qml.device("default.qubit", wires=n_qubits, shots=None) # define a device
#dev = qml.device("lightning.qubit", wires=n_qubits, shots=None) # define a device
# Note: lightning.qubits is faster but "pip install PennyLane-Lightning" is required.

# Initial circuit parameters
var_init = np.random.uniform(high=2 * np.pi, size=(n_layers, n_qubits, 3))

# Definition of a device
@qml.qnode(dev, diff_method='adjoint')
#@qml.qnode(dev, diff_method='adjoint', mutable=False)
#Note: If set the mutable option as False, we can get speeding up but the circuit stracture should be fixed.

# Data encoding and variational ansatz
def quantum_neural_net(var, x):
    qml.templates.AngleEmbedding(x, wires=range(n_qubits))
    qml.templates.StronglyEntanglingLayers(var, wires=range(n_qubits))
    return qml.expval(qml.PauliZ(0))

def square_loss(desired, predictions):
    loss = 0
    for l, p in zip(desired, predictions):
        loss = loss + (l - p) ** 2
    loss = loss / len(desired)
    return loss

def cost(var, features, desired):
    preds = [quantum_neural_net(var, x) for x in features]
    return square_loss(desired, preds)

opt = qml.AdamOptimizer(0.1)
import time

hist_cost = []
var = var_init
for it in range(10):
    t1 = time.time() 
    var, _cost = opt.step_and_cost(lambda v: cost(v, X, Y), var)
    t2 = time.time() 
    elapsed_time = t2-t1
    print("Iter:"+str(it)+", cost="+str(_cost.numpy()))
    print(f"Time：{elapsed_time}")
    hist_cost.append(_cost)

Iter:0, cost=0.614667912125593 Time：12.225000143051147 Iter:1, cost=0.614667912125593 Time：11.945000171661377

scipy.minimize.optimize

# Initial circuit parameters
var_init = np.random.uniform(high=2 * np.pi, size=(n_layers*n_qubits*3)) # one-dimensional array

@qml.qnode(dev, diff_method='adjoint')
def quantum_neural_net(var, x):
    var_3d_array = np.reshape(var,(n_layers,n_qubits,3))
    qml.templates.AngleEmbedding(x, wires=range(n_qubits))
    qml.templates.StronglyEntanglingLayers(var_3d_array, wires=range(n_qubits))
    return qml.expval(qml.PauliZ(0))
        
from scipy.optimize import minimize
hist_cost = []
var = var_init

count = 0
def cbf(Xi):
    global count
    global hist_cost
    count += 1
    cost_now = cost(Xi,X,Y)
    hist_cost.append(cost_now)
    print('iter = '+str(count)+' | cost = '+str(cost_now))
    
result = minimize(fun=cost, x0=var_init, args=(X,Y) ,method='Nelder-Mead', callback=cbf, options={"maxiter":200})
t2 = time.time() 
elapsed_time = t2-t1
print(f"Time：{elapsed_time}")
hist_cost.append(_cost)
print(f"Time per iteration : {elapsed_time/50}")

iter = 1 | cost = 1.004932862586646 iter = 2 | cost = 1.004932862586646 iter = 3 | cost = 1.004932862586646 iter = 4 | cost = 1.004932862586646

w/Tensorflow.keras

Hybrid NN (QNN+NN)

import tensorflow as tf
import keras_metrics

n_qubits = 2
layers = 2
data_dimension = 1 # output
param = {'num_epochs': 128}

dev = qml.device("lightning.qubit", wires=n_qubits)
@qml.qnode(dev, diff_method='adjoint')

def qnode(inputs, weights):
    qml.templates.AngleEmbedding(inputs, wires=range(n_qubits))
    qml.templates.StronglyEntanglingLayers(weights, wires=range(n_qubits))
    return [qml.expval(qml.PauliZ(i)) for i in range(n_qubits)]


weight_shapes = {"weights": (layers, n_qubits,3)}

qlayer = qml.qnn.KerasLayer(qnode, weight_shapes, output_dim=n_qubits)
clayer1 = tf.keras.layers.Dense(n_qubits, activation='linear')
clayer2 = tf.keras.layers.Dense(data_dimension, activation="linear")
model = tf.keras.models.Sequential([clayer1,qlayer,clayer2])

opt = tf.keras.optimizers.Adam(learning_rate=0.01)
model.compile(opt, loss='mse')

hist = model.fit(X, Y, epochs=param['num_epochs'], validation_split=0.1, verbose=1, shuffle='True', batch_size=1024)

loss = hist.history['loss']
val_loss = hist.history['val_loss']

# lossのグラフ
plt.plot(range(param['num_epochs']), 10*np.log10(loss), marker='.', label='loss')
plt.plot(range(param['num_epochs']), 10*np.log10(val_loss), marker='.', label='val_loss')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('epoch')
plt.ylabel('loss')
plt.show()

pred = model.predict(X_test)

data-reuploading

dev = qml.device("lightning.qubit", wires=n_qubits)

@qml.qnode(dev, diff_method='adjoint', immutable=False)
def qnode(inputs, weights):
    weights_each_layer = tf.split(weights, num_or_size_splits=layers, axis=0)
    for i in range(layers):
        qml.templates.AngleEmbedding(inputs, wires=range(n_qubits))
        qml.templates.StronglyEntanglingLayers(weights_each_layer[i], wires=range(n_qubits))
    return qml.expval(qml.PauliZ(0))

weight_shapes = {"weights": (layers, n_qubits,3)}
    
qlayer = qml.qnn.KerasLayer(qnode, weight_shapes, output_dim=n_qubits)
clayer1 = tf.keras.layers.Dense(n_qubits, activation='linear')
clayer2 = tf.keras.layers.Dense(data_dimension, activation="linear")
#model = tf.keras.models.Sequential([clayer1,qlayer,clayer2])
model = tf.keras.models.Sequential([qlayer])

opt = tf.keras.optimizers.Adam(learning_rate=0.05)
model.compile(opt, loss='mse')

hist = model.fit(X, Y, epochs=3, validation_split=0.1, verbose=1, shuffle='True', batch_size=32)

Classifier

Hybrid NN (QNN+NN)

import tensorflow as tf
import keras_metrics

n_qubits = 2
layers = 4
data_dimension = 3 # output
param = {'num_epochs': 128}


dev = qml.device("default.qubit", wires=n_qubits)


@qml.qnode(dev, diff_method='adjoint')
def qnode(inputs, weights):
    qml.templates.AngleEmbedding(inputs, wires=range(n_qubits))
    qml.templates.StronglyEntanglingLayers(weights, wires=range(n_qubits))
    return [qml.expval(qml.PauliZ(i)) for i in range(n_qubits)]


weight_shapes = {"weights": (layers, n_qubits,3)}

qlayer = qml.qnn.KerasLayer(qnode, weight_shapes, output_dim=n_qubits)
clayer1 = tf.keras.layers.Dense(n_qubits, activation='relu')
clayer2 = tf.keras.layers.Dense(3, activation="softmax")
model = tf.keras.models.Sequential([clayer1,qlayer,clayer2])

opt = tf.keras.optimizers.Adam(learning_rate=0.05)
model.compile(opt, loss='categorical_crossentropy')

hist = model.fit(x_train, y_train, validation_split=0.1, epochs=param['num_epochs'], verbose=1, shuffle='True', batch_size=30, callbacks=[early_stopping])

loss = hist.history['loss']
val_loss = hist.history['val_loss']

# lossのグラフ
plt.plot(range(len(loss)), 10*np.log10(loss), marker='.', label='loss')
plt.plot(range(len(val_loss)), 10*np.log10(val_loss), marker='.', label='val_loss')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('epoch')
plt.ylabel('loss')
plt.show()

QSVM

dev_kernel = qml.device("lightning.qubit", wires=n_qubits)

#|00...0><0.....00|
projector = np.zeros((2**n_qubits, 2**n_qubits))
projector[0, 0] = 1

@qml.qnode(dev_kernel)
def kernel(x1, x2):
    """The quantum kernel."""
    AngleEmbedding(x1, wires=range(n_qubits)) #S(x)
    qml.adjoint(AngleEmbedding)(x2, wires=range(n_qubits)) #S^{\dagger}(x')
    return qml.expval(qml.Hermitian(projector, wires=range(n_qubits)))

def kernel_matrix(A, B):
    """Compute the matrix whose entries are the kernel
       evaluated on pairwise data from sets A and B."""
    return np.array([[kernel(a, b) for b in B] for a in A])

svm = SVC(kernel=kernel_matrix).fit(X_train, y_train)

predictions = svm.predict(X_test)
accuracy_score(predictions, y_test)

QAA

from pennylane.templates import ApproxTimeEvolution

T = 2
N = 8
dt = T/N
t = np.arange(0,T,dt)

n_wires = 4
wires = range(n_wires)

dev = qml.device('default.qubit', wires=n_wires)

coeffs = [1, 1, 1, 1]
op_ref = [qml.PauliX(0), qml.PauliX(1), qml.PauliX(2), qml.PauliX(3)]
op_cost = [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliZ(1) @ qml.PauliZ(2), qml.PauliZ(2) @ qml.PauliZ(3), qml.PauliZ(3) @ qml.PauliZ(0)]

hamiltonian_ref = qml.Hamiltonian(coeffs, op_ref)
hamiltonian_cost = qml.Hamiltonian(coeffs, op_cost)

@qml.qnode(dev)
def circuit():
    for i in range(n_wires):
        qml.Hadamard(wires=i)
    for tt in t:
        ApproxTimeEvolution(hamiltonian_ref, (1-tt/T)*dt, 1)
        ApproxTimeEvolution(-1*hamiltonian_cost, tt/T*dt*3, 1)
    return [qml.probs(wires = range(n_wires))]

probs = circuit()
probs = probs.numpy()
probs = probs.flatten()

# プロットする
import matplotlib.pyplot as plt
%matplotlib inline

## z方向に射影測定した時に得られる可能性があるビット列
z_basis = [format(i,"b").zfill(n_wires) for i in range(probs.size)]

plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()

PennyLane-braket周り

braketでaws上タスク回収

def Retrieve_task(task_id):
    from braket.aws import AwsQuantumTask
    # recover task
    task_load = AwsQuantumTask(arn=task_id)

    # print status
    status = task_load.state()
    print('Status of (reconstructed) task:', status)
    print('\n')
    # wait for job to complete
    # terminal_states = ['COMPLETED', 'FAILED', 'CANCELLED']
    if status == 'COMPLETED':
        # get results
        results = task_load.result()
        # print(rigetti_results)

        # get all metadata of submitted task
        metadata = task_load.metadata()
        # example for metadata
        shots = metadata['shots']
        machine = metadata['deviceArn']
        # print example metadata
        print("{} shots taken on machine {}.\n".format(shots, machine))

        # get measurement counts
        counts = results.measurement_counts
        print('Measurement counts:', counts)

        # plot results: see effects of noise
        plt.bar(counts.keys(), counts.values());
        plt.xlabel('bitstrings');
        plt.ylabel('counts');
        plt.tight_layout();
        plt.xticks(rotation=90);

    elif status in ['FAILED', 'CANCELLED']:
        # print terminal message 
        print('Your task is in terminal status, but has not completed.')

    else:
        # print current status
        print('Sorry, your task is still being processed and has not been finalized yet.')

Pennylane側設定例

dev_remote_ionq = qml.device(
    "braket.aws.qubit",
    device_arn=device_arn_ionq,
    wires=n_qubits,
    s3_destination_folder=s3_folder,
    shots = 1024,
    poll_timeout_seconds = 3
)

Pennylaneでaws上マシンを叩いたときのタスクID確認
ionq_task_id = dev_remote_ionq.task.id

PennyLane-qiskit周り

dev = qml.device('qiskit.aer', wires=2, shots=1024)

実機は、

from qiskit import IBMQ
IBMQ.load_account()
dev = qml.device('qiskit.ibmq', wires=2, backend='ibmqx2')

または
dev = qml.device('qiskit.ibmq', wires=2, backend='ibmqx2', ibmqx_token="Your token")