量子アニーリングによる作曲 ( メロディの生成 ) の再現実験を行います。具体的には、OpenJij を用いて複数個の音符を同時生成し、メロディが生成されることを確認します。また、D-Wave Ocean の機能を使用して、マルコフ確率場を用いたメロディの生成を行います。
文献情報
- タイトル : Music Composition Using Quantum Annealing
- 著者 : Ashish Arya, Ludmila Botelho, Fabiola Cañete, Dhruvi Kapadia, Özlem Salehi
- 書誌情報 : https://doi.org/10.48550/arXiv.2201.10557
問題
変数の定義
まずは、変数 $x_{i,j}$ を以下のように定義します。
$$
\begin{equation}x_{i,j} = \left \{\begin{array}{l}1 \ \ {\rm if}\ \ (\ i\ 番目の音符は\ j\ ) \\0 \ \ {\rm otherwize} \end{array}\right.\end{equation}
$$
ここで、$i \in [n]$ かつ $j \in P$ とします。$[n]$ は集合 $\{ 1,2,\dots n \}$ で、$n$ は生成する音符の数を表します。また、$P$ はピッチ ( $C$ や $D$ などの音名 )の集合です。
制約
次に、メロディの生成に必要な 4 つの音楽規則を、最適化問題における制約として定式化します。
制約 1 : $i$ 番目の音符のピッチは 1 つだけ
$n$ 個の音符それぞれに 1 つずつピッチを指定する制約です。この制約は以下のように定式化することができます。
$$ \sum_{j\in P} x_{i,j}=1, \forall i \in [n] $$
この等式制約は、以下の最小化問題に緩和することができます。
$$\min_{\vec{x}} \ \lambda_1 \sum_{i\in [n]} \left( 1-\sum_{j\in P}{x_{i,j}} \right)^2$$
制約 2 : 同じ音を 3 回以上連続して出さない
3 回以上同じ音が連続しないようにする制約です。この制約は以下のように定式化することができます。
$$\exists i \in [n],\ \ \forall j \in P,\ \ x_{i,j}=x_{i+1,j}=x_{i+2,j}=1$$
今、$x \in \{0,1 \}$ であるので、上の式は以下のように表されます。
$$\forall i \in [n-2],\ \ \forall j \in P,\ \ x_{i,j}+x_{i+1,j}+x_{i+2,j} \le 2$$
この不等式は、補助変数 $y_{i,j} \in \{ 0,1,2\}$ を用いて以下の最小化問題に緩和することができます。
$$ \min_{\vec{x}} \ \lambda_2 \sum_{i\in [n-2]} \sum_{j \in P} \left(x_{i,j}+x_{i+1,j}+x_{i+2,j}-y_{i,j} \right)^2 $$
さらに、補助変数 $ y \in \{ 0,1,2\} $ を新たな 2 つの二値変数 $ z_{i,j},z’_{i,j} \in \{0,1\}$ を用いて表すことにより、最終的に以下のような最小化問題になります。
$$ \min_{\vec{x}} \ \lambda_2 \sum_{i\in [n-2]} \sum_{j \in P} \left\{ \left(x_{i,j}+x_{i+1,j}+x_{i+2,j}-2^1z_{i,j}-2^0z’_{i,j}\right)^2 + \lambda_2’z_{i,j}z’_{i,j} \right\}$$
制約 3 : 最初と最後の音を指定
最初と最後の音符を音階の 1 度 ( 今回は $C$ ) に指定する制約です。この制約は以下のように定式化することができます。
$$x_{1, 1}=x_{n,1}=1$$
この等式制約は以下の最小化問題に緩和することができます。
$$\min_{\vec{x}} \ \lambda_3 (-x_{1,C}-x_{n,C})$$
制約 4 : 傾性音の取り扱い
音階の ｢ 2、4、6 度の次には 1 度だけ下の音 ｣ を、｢ 7 度の次には 1 度だけ上の音 ｣ を指定する制約です。具体的には、$D$ の後には $C$ を、$E$ の後には $D$ を、 $G$ の後には $F$ を、 $B$ の後には $C$ が鳴るようにします ( 今回は 1 オクターブの半音階を使用するので、$B$ の後はオクターブ上の $C$ ではなく 1 度の C に戻します )。
$$\begin{equation} \left \{\begin{array}{l} x_{i,D}=1\ ならば\ x_{i+1,C}=1 \\ x_{i,E}=1\ ならば\ x_{i+1,D}=1 \\ x_{i,G}=1\ ならば\ x_{i+1,F}=1 \\ x_{i,B}=1\ ならば\ x_{i+1,C}=1 \end{array}\right.\end{equation}$$
この等式制約は以下の最小化問題に緩和することができます。
$$\min_{\vec{x}} \ \lambda_4 \sum_{i \in [n-1]} \left\{ x_{i,D}\ \left(1-x_{i+1,C} \right) + x_{i,E}\ \left(1-x_{i+1,D} \right) + x_{i,G}\ \left(1-x_{i+1,F}\right) + x_{i,B}\ \left(1-x_{i+1,C}\right) \right\}$$
報酬
特定の音符のペアが現れた時に目的関数を減少させる報酬項を追加します。
$$\min_{\vec{x}} \left( -\lambda’ \sum_{ \substack {i \in [n-1] \\ j, j’ \in P} }{W_{j, j’} \ x_{i,j} \ x_{i+1, j’}}\right)$$
$W_{j, j’}$ は $i$ 番目の音が $j$、$i+1$ 番目の音が $j’$ であったときの報酬の大きさ ( 重み ) を表す係数行列です。この行列は既存の楽曲を分析し、連続するピッチのペアを調べることによって作成します。
なお、これまで導入してきた制約項・報酬項の先頭にある係数 $\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda’ $ は、各項の影響の強さを決める係数です。
マルコフ確率場
マルコフ確率場とは
マルコフ確率場とは、マルコフ特性を満たす、無効グラフで表される確率変数の集合のことを言います。確率変数はグラフの頂点に割り当てられ、辺によって繋がった頂点同士はポテンシャルと呼ばれる値によって依存関係を持つことになります。
以下にマルコフ確率場の一例を示します。
マルコフ確率場の一例 ( https://doi.org/10.48550/arXiv.2201.10557 )
マルコフ確率場では「クリーク」と呼ばれる、全ての頂点 ( ノード ) が繋がっている部分集合を考えます。上のグラフでは、{ A, B } や { C, E }, { A, C, D } などがクリークにあたります。
マルコフ確率場におけるエネルギーは、これらクリークに対して定義されるポテンシャルの合計です。このエネルギーが低いほど、ある確率変数の組み合わせが起こる確率は高くなります。
例えば、上のグラフのクリーク { C, E } を選び、C $\in$ { 0, 1 }, E $\in$ { 0, 1 } とすると、ポテンシャルは以下のように定義できます。
| C | D | p(C,E) |
|---|---|---|
| 0 | 0 | 0.3 |
| 0 | 1 | 0.9 |
| 1 | 0 | 2.6 |
| 1 | 1 | 5.0 |
マルコフ確率場では、最も発生確率の高い ( エネルギーの低い ) 構成に注目します。そして、関心のある状態 ( 今回は変数 $x_{i,j}$ ) をノードとして符号化し、量子アニーリングなどの最適化アルゴリズムを用いて最もエネルギーの低い状態が得られるようにポテンシャルを定義します。
マルコフ確率場と量子アニーリング
イジングモデルは、クリークのサイズが 1 または 2 で、状態空間が { -1, +1 } であるようなマルコフ確率場の特殊ケースと見なすことができます。この関連性から、D-Wave のモジュールを用いることで、マルコフ確率場のポテンシャルから QUBO モデルを作成することができます。
例えば、$\phi_{ij}$ を $ X_1=i $ と $ X_2=j $ に対するポテンシャルとすると、ポテンシャル $ \phi_{00}, \ \phi_{01} \ \phi_{10} \ \phi_{11} $ を持つ確率変数 $X_1$ および $X_2$ に対応する QUBO 形式は以下のように与えられます。
$$
\left(\phi_{10}-\phi_{00}\right)x_1+\left(\phi_{10}-\phi_{00}\right)x_2+\left(\phi_{11}-\phi_{10}-\phi_{01}+\phi_{00}\right)x_1x_2+\phi_{00}
$$
従って、全ての繋がっているノード間のポテンシャルを定義しさえすれば、D-Wave のモジュールでマルコフ確率場を生成することができ、さらに、量子アニーリングによってエネルギーが最も低い確率変数の組み合わせを取得することができます。
実験
実験 0 : 準備
ライブラリのインストール
はじめに、実験を行うために必要なライブラリをインストールします。
!pip -q -q -q install numpy
!pip -q -q -q install matplotlib
!pip -q -q -q install mido python-rtmidi # Midi形式の楽曲を処理するためのライブラリ
!pip -q -q -q install openjij
!pip -q -q -q install dwave-ocean-sdk
!pip -q -q -q install pyqubo
!pip -q -q -q install --upgrade music21
import re
import numpy as np
import IPython
import matplotlib.pyplot as plt
import openjij as oj
import dwave_networkx as dnx
import networkx as nx
from mido import MidiFile
from pyqubo import Array, Constraint, solve_qubo, Placeholder
続いて、楽譜を五線譜として表示するための設定を行います。インストールの途中でエンターの入力を求められるので、エンターキーを押してインストールを続行します。この部分はこちらのページを参考にしました。[1]
!sudo add-apt-repository ppa:mscore-ubuntu/mscore-stable
!sudo apt-get update
!sudo apt-get install musescore -y
# 仮想フレームバッファの設定
!apt-get install xvfb -y
!sh -e /etc/init.d/x11-common start
import os
os.putenv("DISPLAY", ":99.0")
!start-stop-daemon --start --pidfile /var/run/xvfb.pid --make-pidfile --background --exec /usr/bin/Xvfb -- :99 -screen 0 1024x768x24 -ac +extension GLX +render -noreset
# パスの設定
from music21 import *
us = environment.UserSettings()
us["musescoreDirectPNGPath"] = "/usr/bin/mscore"
us["musicxmlPath"] = "/usr/bin/mscore"
us["directoryScratch"] = "/tmp"
楽曲の再生
楽曲を再生するための関数 play_music を定義します。なお、この関数はこちらのページを参考にして作成しました。[2]
freqs = [0] + [880.0 * 2.0 ** ((i - 9) / 12.0) for i in range(12)]
notes = ["R", "C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"]
# 音の生成
def generate_sound(note, bpm, rate=44100):
dic = {s: i for i, s in enumerate(notes)}
freq = freqs[dic[note]]
qn_duration = 60.0 / bpm
t = np.linspace(0.0, qn_duration, int(rate * qn_duration))
signal = 0.5 * np.sin(2 * np.pi * freq * t)
return signal
# 音の再生
def play_music(score, bpm):
melody = np.concatenate(
[generate_sound(note, bpm) for note in score]
) # 全体の音の合成
return IPython.display.Audio(data=melody, rate=44100)
実験 1 : 短いフレーズの使用
実験 1 では 12 個の音符からなる短いフレーズを使用して重み行列 $ W_{j, j’} $ を作成し、それを元にメロディを生成します。
変数と集合の定義
まずは、二値変数 $x_{i,j} \in \{ 0,1 \}$ および補助変数 $ z_{i,j},\ z’_{i,j} \in \{0,1\}$ を定義します。
ここで、$i \in [n], \ j \in P$ であるので、初めに生成する音符の数 $n$ とピッチの集合 $P$ を定義します。今回は、$n=52$、$P=\left\{ C,C\#,D,D\#,E,F,F\#,G,G\#,A,A\#,B \right\}$ とします。従って、50 個の音符からなるメロディを生成します。また、その音符は 1 オクターブの半音階から選択されます。
n = 52 # 生成する音符の数
Pitch = [
"C",
"C#",
"D",
"D#",
"E",
"F",
"F#",
"G",
"G#",
"A",
"A#",
"B",
] # ピッチの集合
len_pitch = len(Pitch) # ピッチの数（12個）
重み行列の作成
次に、重み行列 $ W_{j, j’}$ を作成します。本実験では、以下のフレーズを使用します。
score1 = [
"E",
"D",
"C",
"E",
"D",
"C",
"C",
"C",
"D",
"E",
"E",
"F",
"E",
"C",
"E",
"E",
"D",
"C",
"C",
]
この曲を上で定義した play_music 関数で再生してみましょう（このサイト上では再生ができないため、音声を鳴らしたい場合には Open in Colab ボタンから Google Colab を開き、コードを実行してみて下さい）。今回はメロディのみについて考えているため、それぞれの音符が鳴る長さは一定です。
play_music(
score1, 120
) # 第一引数にリスト， 第二引数にbpm ( 曲を流す速さ ) の数字を入れる
このフレーズの隣り合う音符のペアを数え上げ、重み行列を作成していきます。
W1 = np.zeros([len_pitch, len_pitch])
for i in range(len(score1) - 1):
W1[Pitch.index(score1[i])][Pitch.index(score1[i + 1])] += 1
作成した重み行列 W1 を確認してみます。
print(W1)
今回のフレーズでは、半音階 $P=\{ C,C\#,D,D\#,E,F,F\#,G,G\#,A,A\#,B \}$ のうち $C,D,E,F$ という 4 つのピッチしか登場していないため、W1 の要素のほとんどが 0 になっていることが分かります。
QUBOの定式化
4 つの制約項と報酬項を足し合わせたコスト関数を QUBO 行列として定義する関数 make_pyqubo_model を定義します。
def make_pyqubo_model(n, len_pitch, W):
x = Array.create("x", shape=(n, len_pitch), vartype="BINARY")
z = Array.create("z", shape=(n, len_pitch), vartype="BINARY")
z_ = Array.create("z_", shape=(n, len_pitch), vartype="BINARY")
# 制約項
H1 = np.sum(
[(1 - np.sum([x[i, j] for j in range(len_pitch)])) ** 2 for i in range(n)]
)
H2 = np.sum(
[
np.sum(
[
(x[i, j] + x[i + 1, j] + x[i + 2, j] - 2 * z[i, j] - z_[i, j]) ** 2
+ 3 * z[i, j] * z_[i, j]
for j in range(len_pitch)
]
)
for i in range(n - 2)
]
)
H3 = -x[0, 0] - x[n - 1, 0]
H4 = np.sum(
[
x[i, 2] * (1 - x[i + 1, 0])
+ x[i, 4] * (1 - x[i + 1, 2])
+ x[i, 7] * (1 - x[i + 1, 5])
+ x[i, 11] * (1 - x[i + 1, 0])
for i in range(n - 1)
]
)
# 報酬項
H5 = -np.sum(
[
np.sum(
[
np.sum(
[W[j, j2] * x[i, j] * x[i + 1, j2] for j in range(len_pitch)]
)
for j2 in range(len_pitch)
]
)
for i in range(n - 1)
]
)
H = (
Placeholder("lambda_1") * Constraint(H1, "H1")
+ Placeholder("lambda_2") * Constraint(H2, "H2")
+ Placeholder("lambda_3") * Constraint(H3, "H3")
+ Placeholder("lambda_4") * Constraint(H4, "H4")
+ Placeholder("lambda_5") * H5
)
return H.compile()
この関数を用いて QUBO 行列を作成します。それぞれの項に掛ける係数は適切に調整します。
model = make_pyqubo_model(n, len_pitch, W1)
feed_dict = dict(
lambda_1=20.0, lambda_2=1.0, lambda_3=2.0, lambda_4=2.0, lambda_5=3.0
) # 各項の係数を設定
bqm = model.to_bqm(feed_dict=feed_dict)
print("num_variables:", bqm.num_variables)
結果の出力
OpenJij の SQASampler ( シミュレーテッド量子アニーリング ) を用いて、上で定義した QUBO 行列に対する最適化問題を解きます。今回は、その結果を 10 回読み出し、その中でコスト関数が最も小さくなった解（ decoded_sample[0] ）を採用します。
sampler = oj.SQASampler()
response = sampler.sample(bqm, num_reads=10).aggregate()
decoded_sample = model.decode_sampleset(response, feed_dict)
得られた解において、制約が守られているかどうかを確認します。
print(decoded_sample[0].constraints(only_broken=True))
制約 2 に対応する項が正の値になっているため、制約 2 が破られていることが分かります。
次に、$\vec{x} ∈ \{ 0,1 \}^n$ で得られた結果をアルファベット表記の楽譜に変換します。
result = []
for i in range(n):
for j in range(len_pitch):
if decoded_sample[0].array("x", (i, j)) == 1:
result.append(Pitch[j])
では、生成されたメロディを聴いてみましょう（このサイト上では再生ができないため、音声を鳴らしたい場合には Open in Colab ボタンから Google Colab を開き、コードを実行してみて下さい）。
play_music(result, 120)
生成されたメロディを五線譜として表示する関数 show_music_score を定義します。関数内では、”tinyNotation: 4/4″ から始まる特定の文字列に書き換えます。
def show_music_score(result):
staff = "tinyNotation: 4/4 "
for res in result:
staff += res.lower() + "4 " # 全て四分音符で表示
cp = converter.parse(staff)
cp.show()
生成されたメロディを五線譜で表示します。
show_music_score(result)
音声データや五線譜から、重み行列の元にしたフレーズに多く登場する $E,D,C$ というパターンが繰り返されていることが分かります。また、$C$ が何度も続く箇所が多く見られます。これは、$C,C$ のペアに対応する重みが W1[0][0] = 3 という最大の値になっていることが要因だと考えられます。しかし、制約 2 の効果により、$C$ が 3 つ以上続くのは 1 箇所だけに制限されています。
実験 2 : 多数の曲の使用
実験 1 では、短いフレーズを使用して重み行列を作成しました。その結果、元にしたフレーズに多く登場するパターンのみを繰り返す、単調なメロディとなってしまいました。
そこで、実験 2 では、重み行列に使用するメロディの素材を増やします。具体的には、短いフレーズの代わりに 2 ~ 3 分程度の曲にし、さらに 9 つの曲を使用します。9 つの曲のメロディから、隣り合う音符のペアの出現回数を数え上げて重み行列に足し込んでいきます。従って、実験 1 に比べて、重み行列に反映される音符のペアの種類が増えることが期待されます。
重み行列の作成
Midi ファイルの読み込み
重み行列 $W_{j, j’}$ を作成するために、既存の楽曲からメロディを読み込みます。今回は Midi 形式の楽曲ファイルを使用し、メロディに含まれる音符のペアの出現回数を数え上げていきます。本実験では、ドビュッシー ( Debussy ) による楽曲の Midiファイルを 9 つ利用します。
wget コマンドを使用し、こちらのページから「debussy.zip」というファイルを読み込みます。[3]
!wget http://www.piano-midi.de/zip/debussy.zip
Zipファイルを解凍し、フォルダに含まれる 9 つのファイル ( DEB_CLAI.MID、deb_menu.mid、DEB_PASS.MID、deb_prel.mid、debussy_cc_1.mid、debussy_cc_2.mid、debussy_cc_3.mid、debussy_cc_4.mid、debussy_cc_6.mid ) をアップロードします。
!unzip debussy.zip
メロディの抽出
読み込んだ 9 つの Midi ファイルからメロディを抽出します。まずは、Midi ファイル内の音符を読み取る関数を定義します。本実験で使用する Midi ファイルでは、tracks[1] が ‘Piano right’ ( 右手 ) パートになっています。従って、tracks[1] に含まれる note ( 音符 ) の列をメロディとみなして抜き出していきます。なお、note は 0 から 127 までのいずれかの値を取り、60 が $C4$ の音になっています。従って、72 は 1 オクターブ上の $C5$ になります。[4]
本実験では 1 オクターブの半音階を使用するため、note の値を 12 で割った余りを求めることによって 0 から 11 までの値に直します。
def convert(mid):
staff = []
for msg in mid.tracks[1]: # 'Piano right' から抜き出す
if msg.type == "note_on":
staff.append(
Pitch[msg.note % len_pitch]
) # note の値を len_pitch = 12 で割った余りを求める
return staff
上で定義した convert 関数を用いて、実際に 9 つの Midi ファイルからメロディを読み出します。
midi_list = [
"deb_menu.mid",
"DEB_CLAI.MID",
"DEB_PASS.MID",
"deb_prel.mid",
"debussy_cc_1.mid",
"debussy_cc_2.mid",
"debussy_cc_3.mid",
"debussy_cc_4.mid",
"debussy_cc_6.mid",
]
# アルファベット表記の音符列
score2 = [convert(MidiFile(midi)) for midi in midi_list]
メロディが抽出できていることを確認するために、score2[0] ( 1 番目のファイル ‘deb_menu.mid’ の内容に対応 ) を表示してみます。
print(score2[0])
きちんとアルファベット表記の音符列になっていることが確認できました。
また、score2[4] ( ‘debussy_cc_1.mid’ に対応 ) の 0 ~ 100 番目の音符を再生してみます（このサイト上では再生ができないため、音声を鳴らしたい場合には Open in Colab ボタンから Google Colab を開き、コードを実行してみて下さい）。
play_music(score2[4][0:100], 240)
重み行列の作成
抽出したメロディから、連続する音符のペアを数え上げていき、重み行列 W2 を構成していきます。
W2 = np.zeros([len_pitch, len_pitch])
for j in range(9): # 9 曲分
for i in range(len(score2[j]) - 1):
W2[Pitch.index(score2[j][i])][Pitch.index(score2[j][i + 1])] += 1
作成した重み行列 W2 を確認してみます。
print(W2)
行列のどの成分も正の値になっていることから、半音階 $P=\{ C,C\#,D,D\#,E,F,F\#,G,G\#,A,A\#,B \}$ で取り得る全てのペア ( $C,D$ と $D,C$ のように順序が逆になっているペアや、$D,D$ のような同じピッチのペアも含む ) が登場していることが分かります。
次に、W2 のヒートマップを表示させます。
plt.imshow(W2)
plt.colorbar()
plt.show()
ヒートマップから、行列の対角成分の値が大きいことが分かります。これは、$C,C$ や $D\#,D\#$ といった同じピッチのペアが多く登場していることを意味します。
QUBOの定式化
実験 1 と同様に、make_pyqubo_model 関数を用いて QUBO 行列を定義し、各項の係数を設定します。今回は重み行列 W2 の要素の最大値が 562 であるため、その値で報酬項を割っています。
model2 = make_pyqubo_model(n, len_pitch, W2)
feed_dict2 = dict(
lambda_1=20.0, lambda_2=1.0, lambda_3=2.0, lambda_4=2.0, lambda_5=3 / 562
)
bqm2 = model2.to_bqm(feed_dict=feed_dict2)
print("num_variables:", bqm2.num_variables)
結果の出力
実験 1 と同様に OpenJij を用いて、上で定義した QUBO 行列に対する最適化問題を解きます。
response2 = sampler.sample(bqm2, num_reads=10).aggregate()
decoded_sample2 = model.decode_sampleset(response2, feed_dict2)
得られた解において、制約が守られているかどうかを確認します。
print(decoded_sample2[0].constraints(only_broken=True))
制約 2 と制約 4 に対応する項が正になっているため、制約 2 と制約 4 が破られていることが分かります。
$\vec{x} ∈ \{ 0,1 \}^n$ で得られた結果をアルファベット表記の楽譜に変換します。
result2 = []
for i in range(n):
for j in range(len_pitch):
if decoded_sample2[0].array("x", (i, j)) == 1:
result2.append(Pitch[j])
生成されたメロディを再生します（このサイト上では再生ができないため、音声を鳴らしたい場合には Open in Colab ボタンから Google Colab を開き、コードを実行してみて下さい）。
play_music(result2, 120)
生成されたメロディを五線譜で表示します。
show_music_score(result2)
実験 1 とは異なり、様々なピッチが登場していることが分かります。また、同じパターンを何度も繰り返している箇所が減りました。加えて、半音階の様々なピッチが登場しています。これは、重み行列に様々なメロディを反映させた結果だと考えられます。
また、音声データを聴くと、実験 1 のメロディに比べて耳障りの悪いメロディになっている印象があります。これは、単なるピッチのペアに対して重み付けする定式化をしているためだと考えられます。すなわち、題材にしたメロディに頻繁に登場したピッチのペアが次々に鳴っているだけで、9 曲のメロディそれぞれの全体的な特徴を捉えているとは言い難いためです。
従って、より良いメロディを生成していくためには、どんなメロディをどれだけ重み行列に入れるか、ということよりも、重み行列 $W$ の根本的な定式化を見直すことが必要だと考えられます。
実験 3 : マルコフ確率場の使用
実験 3 では、D-Wave の機能を用いてマルコフ確率場を生成し、その最低エネルギー状態を OpenJij によって求めます。これにより、マルコフ確率場の各ノードに割り当てた確率変数の値を { 0, 1 } で取得します。従って、各ノードに変数 $x_{i,j}$ を割り当てることにより、メロディを生成することができます。
D-wave Ocean には、繋がっているノード間のポテンシャルを定義することにより、それに応じたマルコフ確率場を生成するモジュールがあります。従って、以下では各ノード間の全ての接続に対するポテンシャルを定義していきます。
なお、今回は制約 1 ( $i$ 番目の音符のピッチは 1 つだけ ) のみが守られるようにポテンシャルを設定し、その他 3 つの制約は考慮しないものとします。また、音符の数を 20 個とします。従って、変数 $x_{i,j}$ にあたるノードの数は 20 ( 音符の数 ) × 12 ( ピッチの数 ) = 240 個になります。
重み行列 $W$ を用いたポテンシャルの定義
実験 2 で作成した重み行列をポテンシャルに反映させることを考えます。ポテンシャルは大きい値であるほど、それに対応するノードのペアに割り当てられた確率変数が 1 を取りやすくなります。
従って、ここでは W2 の逆数を取ることにします。W2 の要素は、それぞれのピッチのペアの出現回数なので、W2 の逆数を取ることにより、相性の良いペアのポテンシャルに小さい値を設定することができます。
values = {}
for i in range(len_pitch):
for j in range(len_pitch):
values[(Pitch[i], Pitch[j])] = 562 / W2[i][j] + 20 # W2 の逆数を取る
なお、実際には W2 の逆数を取るだけでなく、562 ( W2 の要素の最大値 ) を掛けて 20 を足しています。これは、制約 1 を守った出力が OpenJij から得られるようにポテンシャルを調整した結果です。従って、これらの数を変えた場合、制約 1 がうまく守られなくなる可能性があります。
上で定義したポテンシャルは、$i$ 番目の音符と $i+1$ 番目の音符に対応するノード間のポテンシャルです。従って、$i$ 番目の音符同士のポテンシャルも定義する必要があります。ただし、制約 1 より、$i$ 番目の音符が 2 つ以上なるのは避けなければなりません。従って、それらのポテンシャルには大きい値 ( 以下では 150 ) を設定します。
n_nodes = 20 # 音符の数
potential = {}
for i in range(n_nodes - 1):
for j in Pitch:
for m in Pitch:
for n in Pitch[Pitch.index(m) + 1 : len_pitch]:
potential[(f"{i}_{m}", f"{i}_{n}")] = {
(0, 0): 50,
(0, 1): 50,
(1, 0): 50,
(1, 1): 150,
} # i 番目の音符のピッチは 1 つだけ
for k in Pitch:
potential[(f"{i}_{j}", f"{i+1}_{k}")] = {
(0, 0): 50,
(0, 1): 50,
(1, 0): 50,
(1, 1): values[(j, k)],
} # i 番目の音符と i+1 番目の音符の間のポテンシャル
for j in Pitch:
for k in Pitch[Pitch.index(j) + 1 : len_pitch]:
potential[(f"{n_nodes-1}_{j}", f"{n_nodes-1}_{k}")] = {
(0, 0): 50,
(0, 1): 50,
(1, 0): 50,
(1, 1): 150,
} # 最後 ( 20 番目 ※ 0 から数えると 19 番目 ) の音符
マルコフ確率場の生成
D-Wave のモジュールを用いて、上で定義したポテンシャルに対応するマルコフ確率場を生成します。
MN = dnx.markov_network(potential)
生成されたマルコフ確率場のグラフを確認してみます。
nx.draw(MN, node_size=1, width=0.01, node_color="magenta")
ノードの数がきちんと 240 個になっていることが確認できます。また、各ノードは 12 個 ( ピッチの個数分 ) のノードからなる層に分かれており、その層が 20 層 ( 音符の個数分 ) あることが分かります。すなわち、右端 ( あるいは左端 ) から数えて $i$ 番目の層は $ i $ 番目の音符のピッチを表すノードの層です。そして、一つ隣にずれると、$i+1$ 番目 ( あるいは $i-1$ 番目 ) の音符のピッチを表す層になります。また、$i$ 番目の層に含まれるノード 12 個と $ i+1 $ 番目の層に含まれるノード 12 個は全て繋がっていることが確認できます。これらの辺それぞれに、上で定義してきたポテンシャルが割り当てられています。
結果の出力
OpenJij を用いて、上で作成したマルコフ確率場のエネルギー ( ポテンシャルの総和 ) を最小にするような確率変数の組み合わせを $\{0,1\} $ で取得します。これにより、1 ~ 20 番目までの音符のピッチが同時生成され、メロディが得られることになります。
今回も 10 回結果を 10 回読み出し、その中でコスト関数が最も小さくなった解を採用します。
sampler = oj.SASampler()
samples = dnx.sample_markov_network(MN, sampler)
print(samples[0])
得られた結果を確認してみます。
for key, value in samples[0].items():
if value == 1:
print(key)
制約 1 を満たす 20 個の音符が得られたことが分かります。ただし、現状では ” ( 音符の位置 )_( ピッチ ) ” という文字列になっているので、( 音符の位置 ) をインデックス、( ピッチ ) を要素とするリストに書き直し、実験 1 や実験 2 と同様の楽譜を作成します。
score3 = [""] * n_nodes
for key, value in samples[0].items():
if value == 0:
continue
temp = re.sub(r"\D", "", key) # 数字部分を取り出す
idx = key.find("_")
score3[int(temp)] = key[
idx + 1 :
] # 数字部分をリストのインデックスにし、"_" の後の文字をそのインデックスの要素にする
では、生成されたメロディを再生してみます（このサイト上では再生ができないため、音声を鳴らしたい場合には Open in Colab ボタンから Google Colab を開き、コードを実行してみて下さい）。
play_music(score3, 120)
続いて、生成されたメロディを五線譜で表示します。
show_music_score(score3)
実験 2 と同じ重み行列 W2 を使用しているのにも関わらず、実験 2 で得られたメロディとは異なる雰囲気のメロディになっていることが分かります。その理由として、今回は W2 をマルコフ確率場のポテンシャルに反映させているため、重み行列の効き方が異なっていることが考えられます。また、実験 3 では制約 1 以外を考慮していないことも一因であると予想されます。
まとめ
実験 1, 2, 3 を通じて、量子アニーリングによるメロディの生成を行ってきました。いずれの実験においても重み行列 $W_{j, j’}$ を使用してきましたが、この $W_{j, j’}$ は隣り合う音符のペアの出現回数をそのまま重みにしているに過ぎないため、実験 1 では、重みの大きい音符のペアを繰り返すだけのメロディが生成されました。また、実験 2 , 3 においても、元にした楽曲のメロディの特徴がうまく掴めているとは言い難い結果となりました。
従って、量子アニーリングを用いてより良いメロディを生成するためには、重み行列 $W_{j, j’}$ の定式化を見直す必要があると考えられます。
あと書き
本実験では、論文に記載の定式化のみを用いましたが、特に重み行列に関しては改善が必要だということが分かりました。今回の実装を通して、このような論文を読むだけでは気づくことができなかった部分が見えてきたという点で、自分で手を動かして再現実験をしてみることの重要性を感じました。