数値計算

ブルースカイ・カタストロフィ

ブルースカイ・カタストロフィ(つよそう)

スクリーンショット 2018-05-05 17.10.41

↑は参考にあった Blue sky orbit in the Gavrilov-A. Shilnikov model でパラメータも同じもの.

以下は, \(\varepsilon\) をそのまま固定しながら, \(\mu\) を \([0.1, 1.1]\) の範囲で動かしたもの(0.02 間隔で 80000 回のルンゲクッタ.最初の 1000 回は捨てている).

ローレンツアトラクタとカオス

こんにちは.TRSasasusu です.

記事を更新していないだけで,KCS は活動しています.新歓の準備も進んでいるようです.

さて,自分は生命情報学科に所属していますが,この前のレポートにてカオスで遊んだので,今回はその話をします.

スクリーンショット 2018-02-01 15.46.11

これはローレンツ方程式(↓の式)のアトラクタです.この軌道に沿って点が移動していきます.


\begin{align}
\frac{dx}{dt} &= -\sigma x + \sigma y \\
\frac{dy}{dt} &= -xz + rx -y \\
\frac{dz}{dt} &= xy – bz
\end{align}

特にこのようなフラクタル構造を持つアトラクタはストレンジアトラクタと呼ばれ,カオスを生み出します.蝶の羽みたいで美しいですね.

そして,パラメータ(今回は \(\sigma,\ r,\ b\))を変化させると,カオスになったりならなかったりします.これを表すのが分岐図です.

スクリーンショット 2018-02-01 16.07.43

なんかすごいですね.そして,カオスを定量的に表すのがリアプノフ指数(多次元でまとめたものがリアプノフスペクトラム)です.力学系を表すヤコビ行列を時系列順に並べて QR 分解して(時刻 \(0\) では時刻 \(0\) のヤコビ行列のみを,それ以降ではヤコビ行列に前時刻の \(Q\) を右から掛けたものをそれぞれ QR 分解する),各時刻の \(R\) のそれぞれの対角要素の \(\log\) をとって全部足し合わせて \(N\) で割ると,以下のようになります.

スクリーンショット 2018-02-02 1.06.55

1 つでも正のリアプノフ指数があればカオスとなるようなので,この結果は分岐図のものと一致していますね.

最後に,ストレンジアトラクタの動画を載せて終わりにします(動画はレポートには載せられない).

量子コンピュータSimulatorを開発しました+量子アルゴリズムの練習をした→悲しい事実に気づく(1)

なぜ量子コンピュータを作ったのか

こんにちは.

今回は,唐突に思いついて作った量子コンピュータシミュレータと,その量子コンピュータを使った量子アルゴリズムの構築法を勉強するための練習についてのお話です.

まず,「量子コンピュータ」と聞くと,一見難しいものと思えるかもしれません.実際,僕も数日前までにはそう思ってました.ところが,数日前に,古典コンピュータ(今皆さんが使っているものです)で数値計算を行うことによって,量子コンピュータの挙動を再現することが可能だということを知りました.このときに,量子コンピュータが実際どういう仕組みで動いているかよくわかっていないことに気付いて,文献を漁って調べてみることにしました.すると,なんだか難しそうだなぁと勝手に思っていたものは,実際にはそんなに難しくなく,再現するための数値計算の手法も難しくないことを知りました.

という経緯があって,量子コンピュータが(数値計算上では)誰にでもつくれるという考えに至り,実際に作ってみようということになりました.

量子コンピュータは何で動くか

古典コンピュータは,AND回路やOR回路などの,複数の古典ビットを何らかの物理的相互作用を用いて変化させることで,様々な情報を処理しています.このことからもわかるように,古典コンピュータの情報の最小単位は,1古典Bitという,0か1かの区別であり,これより小さな情報の単位は存在しません.普段我々が目にしている古典コンピュータは,古典Bitを用いてどのように情報を処理しているのか見ることはありませんが,中身は古典Bitの膨大な操作によって,動いています.これらは,古典Bitをどのように操作するかの技術(≒アルゴリズム)によって支えられています.

一方,量子コンピュータの動作は古典Bitによるものではありません.量子コンピュータの情報の最小単位は「量子ビット(qubit)」と呼ばれるものです.量子ビットは,0か1かの厳密な区別である古典ビットと違い,0か1かの確率的な状態を保持しています.量子ビットは古典ビットではないので,表記法も異なります.例えば,1量子ビットは,a|0>+b|1>のように表記します(ディラックのブラケット記法).ここで,a,bは複素数です.つまり,量子ビットは,0と1がそれぞれa,bだけ含まれる,という状態を保持することができるのです.もちろん,量子ビットは量子(の性質を持つもの)ですから,「観測」を行うことで0か1かが確定します.すなわち,|0>か|1>になるということです.どちらに確定するかは,もともとの量子ビットの状況(a|0>+b|1>)に応じて確率的に決まります.具体的には,a|0>+b|1>である量子ビットが|0>になる確率は,|a|^2で,|1>になる確率は,|b|^2です(このことからもわかるように,|a|^2 + |b|^2 = 1 が満たされる必要があります).ここで重要なポイントは,我々は量子ビットの状態を「観測」せずに知ることができないということです.これは,a,bを直接的に知ることができないということを表しています.

もう一つ,量子ビットについての重要なポイントがあります.それは,量子の「もつれ(entanglement)」です.量子のもつれとは,簡単に言うと,複数の量子Bitがあったとき,複数の量子Bitを観測したときの確率論的独立性がない状態のことです.具体例をあげてみます.今,2つの量子ビットX,Yがあったとします.X,Yがもつれていないときは,XとYを観測したとき,観測されたXの値が0,1どちらであっても,観測されたYが0,1を取る確率は,変わりません.これは,2つのさいころをどうじに投げた状況に似ています.一方のサイコロの目が1だろうが,2だろうが,もう一方のサイコロの目は等しく1/6ですよね.しかしながら,X,Yがもつれているときは,Xの値が0,1どちらになるかによって,Yが0,1どちらになるかの確率が変わってしまうのです.サイコロの例でいうと,2つのサイコロを同じ向きで完全に固定したときに似ています.一方のサイコロの値が1であったら,もう一方は1でしょう.もちろん,一方のサイコロの値が1であったら,もう一方は奇数が出やすくなる,のような例もありですね.このようなもつれ,の状態は,もはや|X>=a|0>+b|1>,|Y>=c|0>+d|1>のような表記はできません.量子がもつれているような状態を表記するためには,|XY>=e|00>+f|01>+g|10>+h|11>のようにします.このような状態の量子ビット群では(X,Y)が(観測したときに)(0,0)を取る確率がe^2,(0,1)を取る確率がf^2,…ということを意味しています.ここまでくればお気づきでしょうが,(任意にもつれさせることのできる)量子ビットがn個あれば,2^n通りの観測結果の確率を表現できます.このことは,非常に重要なことです.すなわち,適切に量子ビットを操作することで,確率的に存在しているあらゆる組み合わせの計算を一度に行える可能性があるということです.

これらの話から,量子コンピュータは,どのようなアルゴリズムで動くのか予想がついた方もいるでしょう.すなわち,量子コンピュータは,量子ビットを操作し,もつれさせ,観測を行うことで,結果を得る,という仕組みです.

量子コンピュータをシミュレートする

と,いうのが量子コンピュータの簡単な説明なのですが,文中に会ったように,結局は,量子の状態は.量子ビット数nに対して,2^n個の複素数で表現できるということです.量子の操作は,この2^n個の複素数を,量子の操作に対応する関数を用いて,写像させることで表現できます.つまり,古典的コンピュータで,量子コンピュータがシミュレート(数値計算で仮想的に実現するということ)することができるということです.よって,作りました(演繹的).

今回は,「量子コンピュータが実装できたよ」ということを示すために,簡単な量子プログラムを作って,結果どうなったかを示したいと思います.そのために,まず(このシミュレータが)量子操作としてどのような操作が行えるのかを,簡単に書いておきます.

.量子ビット数は11qbit.(なぜ,11なのかは,特別な理由はない)

・量子ビット初期化:量子ビットを,任意の確定値にする(|XYZ…>=|01010…>)

・1つの量子ビット(任意)を,パウリ行列(X,Y,Z)とかアダマール行列(H)を用いて変換する.(詳しくは後述)

・1つの制御量子ビット,1つの目標量子ビットを,パウリ行列(X,Y,Z)とかアダマール行列(H)を用いて,制御変換する.(詳しくは後述)

・2つの制御量子ビット,1つの目標量子ビットを,パウリ行列(X)を用いて,制御制御変換する.(詳しくは後述する.)

imagee
量子コンピュータシミュレータ(開発中)の画面

これらの操作機能を持った,量子コンピュータを使って,(僕が勝手に考えた問題)4=x+yを解いてみようと考えた.

…のだが,量子プログラムを書いてからとても残念なことに気づいてしまったのである.

文章が予想以上にながくなったので,これを前半として,前後半に分けて投稿することにする.つづきは,後半で.

 

 

 

 

フレーム独立GMM-based mappingによる声質変換

こんにちは。1年のSannkoです。AI班です。

初投稿になります。

宜しくお願いします。

 

UMUさんの記事にもありますが、AI班は声質変換の活動をしています。

まだどのような手法で声質変換をしていくかは検討中ですが、既存の手法を試してみようと思います。

 

戸田さんの論文の2章に書いてある、提案手法ではない古い手法をやってみます。

この手法ではフレームごとにGMMでソースとターゲットの同時分布を推定して、そこから条件付き確率やら周辺確率が出せるのでそれを使って変換を行います。

フレーム間の関係は全く考えられていないので、最近の手法と比べるとかなり古典的な手法ですね。

 

https://pdfs.semanticscholar.org/d419/ceb2753232373fd4ab9534b371e017cd9dc1.pdf

 

データはこのサイトのものを使わせてもらいました。

ありがたいですね、こういう研究室は。

8.00モーラ/秒の25文のデータを使っています。

 

http://www.it.ice.uec.ac.jp/SRV-DB/

 

女性1から男性への変換をやってみたいと思います。

変換元の女性の声はこれ

 

 

ターゲットはこの男性です。

 

 

この文章のデータは訓練データからは除いてあります。

変換を行った結果がこれです。

基本周波数はそのままなので、声が高くなっているのがわかります。

ただ、声質は男性のものに近くなっている気がしますね(定性的)。

 

 

基本周波数の変換についてはフィルタとの関係などを考えてきちんと検討すべきですが、(正直よく分からないので)今回はとりあえず単純に切片0の回帰曲線でモデル化しました。

単純ですが、画像をみるとわかる通り、それなりに妥当です(定性的2)。

(追記)貼ってから気づいたんですが、この画像は切片を0にする前のやつでした。下のデータに使ったモデルはちゃんと切片0になっています。

figure_1-%e3%81%ae%e3%82%b3%e3%83%94%e3%83%bc

このモデルを使って基本周波数も変換をした結果がこれです。

ピッチを適当に扱ったのでノイズが増えてしまったような気がします(定性的3)。

しかしかなりターゲットの声に近づいたと思います(定性的4)。

 

 

正直、この単純な手法でここまで変換できたので驚いています。

今後は戸田さんの提案手法など、もう少し高度な手法を試してみたいですね。

音声分析合成システムWORLDをPythonに移植した

最近,音声の声質変換を行うために,Pythonで使えるライブラリを探していました.matlab上で使える便利な音声分析合成ソフトとして「WORLD」が存在しますが,これをPython上で使おうとすると少々めんどくさいということがありました.そこで,科学研究用のスタンダードなPython開発環境Anacondaがあれば簡単にWORLDを使えるように,matlabのコードを参考にPythonへ完全移植を行っています.とりあえず,基本的な分析→合成までの部分の移植が完了したので,移植したコードがまともに動くかどうかの検証も含めて記事にしてみました.

 

まずは,元の音声データです.この音声データをもとに,パラメータを抽出(分析)し,さらに再合成することで音声データが再構築されます.再構築する際に,パラメータを統計的に変化させて別人の声を生成する方法が,「統計的声質変換」です.

 

次に,正規版(matlab)の変換結果です.

 

最後に,移植版(Python)の変換結果です.

 

聞いてみると,正規版と移植版の違いが分かると思いますが,これは,matlab版で使用されている一部の関数が,pythonでは存在しないために,なるべく同じような挙動をするようにコードを書いていますが,やや違いがでてきてしまうためです.

なお,この移植版を公開していいかわからないため,公開するかどうかは後日決定します.

今後,移植を完成させつつ,声質変換についての記事も書いていこうと思います.BYE!

 

 

–ここから残念なお話–

まず,matlabからpythonに移植するとき,配列のインデックスの書き方の違いで,とても混乱した.matlabでは,配列のはじめは1から始まるが,Pythonでは一般的な言語と同じく0である.これだけなら,数値を1ずらすだけでいいのだが(実際にはそれだけでもかなりのストレスである.),配列から一部の列を取り出す際,例えば,[0,1,2,3,4,5]という配列から[2,3,4]を取り出すとき,取り出しの記述はPythonでは2:5(2番目から5番目の意)となり,Matlabでは3:5(3番目から5番目の意)とかく(Pythonでは,一番最後のインデックスは無視される.)これが非常に頭を混乱させる要因となっており,一番混乱したときは,matlabのコードを実行したときの図を紙に書いて,紙に書かれた内容をPythonコードで書くという所業をしてしまった.

さらに,面倒くさいのが,matlabで定義されている関数を,pythonで移植する場合,matlabの関数となるべく近い挙動をするpythonコードを書く必要がある.pythonで完全一致する関数があればいいのだが(うれしいことに,汎用的なものは共通化されているようだ),すべてきれいに一致するとは限らないので,なるべく同じ計算結果となるようにコードを書く必要があった.

一部は,面倒くさいので放置しているが,聞こえ上は似ているものができたので,まずは良しとすることにした.

 

LOGSPECTROGRAMさむねいる

反応拡散系シミュレーション

Simulation 1

\[\frac{\partial{u}}{\partial{t}}=20 Hill(u,0.2,2)/v – 80U +1 \]

\[\frac{\partial{v}}{\partial{t}}=20 (u-v) \]

Simulation 2

\[\frac{\partial{u}}{\partial{t}}=20 Hill(u,0.2,2)/v – 80U +1 \]

\[\frac{\partial{v}}{\partial{t}}=30 (u-v) \]

Simulation 3

\[\frac{\partial{u}}{ \partial{t}}=0.5 \Delta u +20 Hill(u,0.2,2)/v – 80U +1 \]

\[\frac{\partial{v}}{ \partial{t}}=5.0 \Delta u +30 (u-v) \]

(dx=0.1,dt=0.0005)

 

samneiluSAMUNEIRU

二次元流体シミュレーターVer-314を開発しました

こんにちは,UMUです.

今日は定例会ですね.新歓などの関係もあって3月下旬の発表から何も作ってないなと思い,二次元流体シミュレータを作ってみました.

矩形内に流体の速度を固定したりすることで流体感が出ます.

計算誤差とか何も考えずに差分法で実装しましたが,意外に綺麗になりました.以下,いろいろ計算結果の動画です.

グラフの見方:

左図:X方向の流体速度 右図:Y方向の流体速度

色は,赤ほど正,青ほど負となるようになっています.

たとえば,左図の赤い所は,X正方向に強い流れがあることを示します.

(左と右が両方青い時は左上に流れているということですね.)

▽動画

最後のやついいかんじ

 

理論

数値計算の理論は(割とテキトーなので)省きます

この流体シミュレータで用いた式は無次元化されたナビエストークス方程式(以下)です.

\[\frac{\partial v}{\partial t} + (v \cdot  \nabla)v = -\nabla p + Re^{-1} \nabla ^2 v   \]

左から,慣性項,移流項,圧力項,粘性項です(Reはレイノズル数)が,誤差お関係から粘性項がなくても拡散が起こってしまいます(数値拡散).次回があれば,数値拡散を抑える方法をやってみたいと思います.

おしまい.

余談

定例会のために何かひねり出した結果がこれ.22時に帰ってきてその後の作業はきついということがわかった.

追加

人工知能講習会(機械学習講習会)やることまとめ

目標

MNISTデータベースの学習

内容目次

全4回,予備1回

1.Python入門

1.1 環境構築(Anaconda + Mako + PyOpenCL,Device for OpenCL)

1.2 Dive Into Python

2.数値計算入門

2.1 NUMPY+SCIPY入門

3 .ニューラルネットワーク入門

3 .1 機械学習の基礎知識

3 .2 ニューラルネットワークの基礎知識

3.3 NNPropagator入門

4.NNPropagatorを使ったニューラルネットワークの学習

4.1 多層パーセプトロンによる関数近似

4.2 MNIST学習

前提条件

手続き型プログラミング言語を少し習得していること.

※手続き型プログラミング言語:C,C++,Python,Java,C#,Javascript…など.

内容

1.Python入門

人工知能プログラミングで必要となる,Pythonについて入門します.

1.1 環境構築(Anaconda + Mako + PyOpenCL,Device for OpenCL)

Pythonでプログラミングできる環境を作ります.KCS製のニューラルネットワークライブラリ「NNPropagator」を動作させるための環境となっています.

1.2 Python入門

Pythonを知らない人向けに,一通り簡単にPythonを学びます.

2.数値計算入門

数値計算を知らない人向けに,一通り簡単にPythonを学びます

2.1 NUMPY+SCIPY入門

Pythonには,優秀な数値計算用のライブラリが多数存在します.そのうち,汎用性が高く習得が容易な「Numpy」について一通り簡単に学びます.Scipyは必要に応じて学びます.

3 .ニューラルネットワーク入門

機械学習でも最先端で活用されている,ニューラルネットワークについて学びます.

3 .1 機械学習の基礎知識

機械学習とは何かについて簡単に学びます.

3 .2 ニューラルネットワークの基礎知識

ニューラルネットワークとは何かについて簡単に学びます.確率と線形代数の知識が多少必要ですが,わからなければ教えます.

3.3 NNPropagator入門

KCS独自開発のライブラリNNPropagatorについて学びます.独自開発のライブラリと言っても,ニューラルネットワークの実装を補助するものです.ライブラリもナイーブに書いており,実装を理解するのも簡単です.

4.NNPropagatorを使ったニューラルネットワークの学習

NNPropagatorを用いて,こちら側で用意したデータベースを解析して頂きます.

今回の講習会では,2系統(関数近似系/MNIST)用意してあります.また,確認用に自分でデータベースを作り,評価します.

4.1 多層パーセプトロンによる関数近似

こちらで用意した学習用のデータ群をNNPropagatorを用いて学習・評価します.

4.2 MNIST学習(数字文字認識)

MNISTデータ群をNNPropagatorを用いて学習・評価します.

MNISTとは,下の画像の様な数字の画像と教師データのデータベースです.

(画像を与えられた時に何の数値であるか推定する人工知能を作ります.)

img679

教科書

1~3.2インターネット上または書籍で適当なものを引用し,活用します.

3.3~ 独自の教科書を用います.

カハンの加算アルゴリズムが最強な件(修正版)

いつものように,流体や機械学習の数値計算のプログラミングを書いているときの話である.

実行中,数値計算が破綻していることに気がついた.

検証の結果,約数千の数値の総和をとるコードの誤差が異常に大きくなることを発見した.

これは,数値の大きさに幅があり,大きい数値が小さい数値の精度を奪ってしまうことに基づく.

すなわち,

1+0.0000000001+0.0000000001+…(1000000000回)..+0.000000001

が,1になってしまうのと同様の現象である(本来は2が正しいが,精度の問題で1+0.0000000001=1となるため).

この誤差をどうにかできないかと調べた結果,有名な方法があるらしい(知らなかった).

それは,カハンの加算アルゴリズムである.

カハンの方法では,誤差を保持することにより,総和による誤差の累積を防ぐようになっている.

擬似コードは以下.

</p>
<p>ARR = 総和を求めたい数値の列<br />
SUM = 0 //総和<br />
ERR = 0 //誤差<br />
for num in ARR:<br />
    y = num &#8211; ERR //誤差の反映<br />
    t = SUM + y  //試しに加算<br />
    ERR = (t &#8211; SUM) &#8211; y  //誤差の算出<br />
    SUM = t<br />

このように,ある計算を行いたい場合,適切なアルゴリズムを用いるなど,十分な精度を持つ方法で実装を行う必要がある.

そうでないと,数値計算で特異な現象を発見したと思ったら,じつは計算誤差のせいでしたということになりかねない(というか,よくなる).

数値計算についてもっと学ばなければいけないと思った.

緑目はなくなるのか?

インターネットを見ていたら,次のような画像があった.

(The color of your future child’s eyes.)

eyet1

これは,目の色の遺伝確率を示しているらしい.

この図の一番下を見ると,青+青→ほぼ青となっており,これだと長い時間たった後,青しか残らないというコメントがあった.

そこで今回は,この図の記述に沿って世代が変わるごとに目の色の分布の比率がどうなるかシミュレートしてみた.

以下,シミュレーションのための簡略した設定を示す.

・第0世代において,茶・緑・青の比率は等しいとする(1/3).

・人数が十分に多いとし,個人の好みは考えない.

この条件で目の色比率がどうなるかを計算すると,次のようになる.

eyecolor

 緑目はなくなりませんね.めでたしめでたし.

(厳密にいうと,緑+緑→茶の確率がゼロでない限り,

茶色もなくなりません)