Griffin-Limアルゴリズムの実行時間をlibrosaとtorchaudioで比較してみた話

はじめに

音声の振幅スペクトルから位相を推定し、元の音声を復元するためのGriffin-Limアルゴリズムが知られている。
Griffin-Limアルゴリズムはlibrosaパッケージとtorchaudioパッケージの両方に実装されている。

• librosa

librosa.org

• torchaudio

pytorch.org

今回はGriffin-Limアルゴリズムに関してCPUとGPUとで実行時間を計測し比較してみた。

環境

• OS、ハードウェア、ドライバ
  • Ubuntu 18.04
  • GeForce RTX 2070
  • CUDA 11.2
  • CUDA Driver 460.32.03

• ソフトウェア
  • Python 3.6.9
  • librosa 0.8.0
  • torch 1.7.1
  • torchaudio 0.7.2
  • torchlibrosa 0.0.9

スクリプト

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.8.1
# torchaudio 0.8.1
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import numpy as np
import torch
import torchaudio
import torchaudio.transforms as T

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_ITER = 50

if __name__ == "__main__":

    waveform, sample_rate = librosa.load(librosa.ex("trumpet"))

    # librosa
    start = time.time()
    specgram = np.abs(
        librosa.stft(
            waveform, n_fft=N_FFT, hop_length=HOP_LENGTH, win_length=WIN_LENGTH
        )
    )
    y_inv = librosa.griffinlim(
        specgram, n_iter=N_ITER, hop_length=HOP_LENGTH, win_length=WIN_LENGTH
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    # load waveform as torch.Tensor
    waveform, sample_rate = torchaudio.load(librosa.ex("trumpet"))
    waveform_cuda = waveform.cuda()

    # torchaudio (CPU)
    specgram = T.Spectrogram(
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    )
    griffin_lim = T.GriffinLim(
        n_fft=N_FFT,
        n_iter=N_ITER,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    )
    start = time.time()
    reconstructed = griffin_lim(specgram(waveform))
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    specgram = T.Spectrogram(
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    ).cuda()
    griffin_lim = T.GriffinLim(
        n_fft=N_FFT,
        n_iter=N_ITER,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    ).cuda()
    torch.cuda.synchronize()
    start = time.time()
    reconstructed = griffin_lim(specgram(waveform_cuda))
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

実行結果

GPUを使ったtorchaudioが一番早い結果となった。

elapsed_time (librosa): 0.264543[sec]
elapsed_time (torchaudio; CPU): 0.147203[sec]
elapsed_time (torchaudio; GPU): 0.022305[sec]

まとめ

位相復元に関しても、torchaudioのCPU版を使うだけで実行時間の短縮が望めそうである。GPU版も悪くない結果である。今回は内部的に実装されているアルゴリズムを考慮したフェアな比較というよりは、いちユーザが関数を外から触るときの体感を重視した。

音声認識結果を音声合成するPythonスクリプトをSpeechRecognitionとPyOpenJTalkで書いたみた話

はじめに

かつて、音声認識と音声合成を組み合わせて遊んでみるという主旨の記事を書いたことがある。

tam5917.hatenablog.com

音声合成には、コマンドラインから音声合成できるOpenJTalkパッケージを用いたのだった。これをPythonから動かす場合には、専用の関数を定義して呼び出す必要があり、ひと手間かかる。

最近、山本りゅういち氏がPyOpenJTalkなるパッケージを整備し、単体でTTS（Text-to-Speech、テキストを音声に変換する）機能を実現した。

https://t.co/Fi7L4SSdeq アップデート終了。デモノートブック作りました。hts_engine_API, open_jtalk を事前にインストールしなくても良いように内部的に変更しています。以前はテキスト処理だけでしたが、波形生成の仕組みも入れて、単体でTTSとして遊べるようになりました。

— 山本りゅういち / Ryuichi Yamamoto (@r9y9) March 21, 2021

これを使って前回記事のPythonスクリプトを少しきれいにしたというのが本記事の主旨である。

インストール

音声認識にはSpeechRecognitionパッケージを使う。PyOpenJTalkが今回のキモなので忘れずにインストールする。最新版をインストールするという意味で--upgradeオプションもつけておく。

pip3 install SpeechRecognition
pip3 install --upgrade pyopenjtalk

スクリプト

以下のスクリプトを適当な名前で保存し、実行する。 「話しかけてみましょう！」と表示されたのち、適当な言葉をしゃべると、認識された言葉をオウム返しで音声合成して終了するというシンプルなものである。

import subprocess

import numpy as np
import pyopenjtalk
import speech_recognition as sr
from scipy.io import wavfile

OUT_WAV = "/tmp/tts.wav"

if __name__ == "__main__":

    # マイクからの音声入力
    r = sr.Recognizer()
    with sr.Microphone() as source:
        print("話しかけてみましょう！")
        audio = r.listen(source)

    try:
        # 日本語でGoogle音声認識
        text = r.recognize_google(audio, language="ja-JP")
    except sr.UnknownValueError:
        print("Google音声認識は音声を理解できませんでした。")
    except sr.RequestError as e:
        print("Google音声認識サービスからの結果を要求できませんでした;" " {0}".format(e))
    else:
        print(text)

        # 音声認識結果を音声合成する
        x, sr = pyopenjtalk.tts(text)

        # wavファイルに保存
        wavfile.write(OUT_WAV, sr, x.astype(np.int16))

        # 保存したwavファイルを再生
        subprocess.run("afplay " + OUT_WAV, shell=True)

おわりに

PyOpenJTalkのおかげでTTSが手軽に実現でき、楽しい。PyOpenJTalkには話速変換やピッチ変換のオプションもあるので、調べて使ってみるのが良いだろう。

おまけ

話速変換とピッチ変換のオプションを使ったnotebookはこちらから。 nbviewer.jupyter.org

注意点

実行時のエラー

File "pyopenjtalk/htsengine.pyx", line 1, in init pyopenjtalk.htsengine
ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 88 from C header, got 80 from PyObject

このようなエラーが出たら、numpyを再インストールしてみるのが吉。

Google音声認識API

ただし、本記事で利用しているGoogle音声認識APIはあくまでテスト/デモ用であり、いつGoogle側にAPI利用を打ち切られても文句は言えない。Google製の音声認識APIを使いたい場合はCould Speech-to-Textの認証情報を取得して使うのが本筋である。ひと月あたり60分の音声までなら利用無料であるが、それ以上を利用するにはお金がかかる。

cloud.google.com

torchaudioとtorchlibrosaの実行速度に違いはあるのか？

はじめに

PyTorchには音声系データを処理するのに便利なtorchaudioというライブラリが存在する。
pytorch.org

一方、音声系データの処理に便利なlibrosaというパッケージが存在する。
librosa.org

さらにtorchlibrosaという、librosa内部の行列計算まわりをPyTorchで置き換えたパッケージが存在する。
github.com

ここで一つ疑問：
「で、結局どれ使えばいいの？（どれが早いの？）」

この疑問が気になったので、CPU実行とGPU実行における実行時間を比較検証するための簡易的なスクリプトを書いて調べてみたということである。

環境

• OS、ハードウェア、ドライバ
  • Ubuntu 18.04
  • GeForce RTX 2070
  • CUDA 11.2
  • CUDA Driver 460.32.03

• ソフトウェア
  • Python 3.6.9
  • librosa 0.8.0
  • torch 1.7.1
  • torchaudio 0.7.2
  • torchlibrosa 0.0.9

スクリプト

メルスペクトログラムを計算するスクリプトである。

• 3月24日追記：GPUの計測方法に誤りがあったため、コメントに従ってtorch.cuda.synchronize()を用いるように修正した。

スクリプトを表示する

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.7.1
# torchaudio 0.7.2
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import librosa.feature
import torch
import torchaudio
import torchaudio.transforms as T
import torchlibrosa as tl

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_MELS = 128

if __name__ == "__main__":

    waveform, sample_rate = librosa.load(librosa.ex("nutcracker"))
    # librosa
    start = time.time()
    mel_spectrogram = librosa.feature.melspectrogram(
        y=waveform,
        sr=sample_rate,
        n_fft=N_FFT,
        hop_length=HOP_LENGTH,
        n_mels=N_MELS,
        power=2.0,
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    waveform, sample_rate = torchaudio.load(librosa.ex("nutcracker"))
    waveform_cuda = waveform.cuda()

    # torchaudio (CPU)
    start = time.time()
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    melspec = mel_spectrogram(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    torch.cuda.synchronize()
    start = time.time()
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    mel_spectrogram = mel_spectrogram.cuda()
    melspec = mel_spectrogram(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (CPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    start = time.time()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    logmel = feature_extractor(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (GPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    torch.cuda.synchronize()
    start = time.time()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    feature_extractor = feature_extractor.cuda()
    logmel = feature_extractor(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

• 5月8日追記：koukyo1213様のコメントに従い、torch.nn.Sequentialやtorchaudio.transforms.MelSpectrogramのインスタンス化の部分を計測対象から外し、計算処理のみにフォーカスしたスクリプトは以下である。

スクリプトを表示する

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.8.1
# torchaudio 0.8.1
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import librosa.feature
import torch
import torchaudio
import torchaudio.transforms as T
import torchlibrosa as tl

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_MELS = 128

if __name__ == "__main__":

    torchaudio.set_audio_backend("sox_io")

    waveform, sample_rate = librosa.load(librosa.ex("nutcracker"))
    # librosa
    start = time.time()
    mel_spectrogram = librosa.feature.melspectrogram(
        y=waveform,
        sr=sample_rate,
        n_fft=N_FFT,
        hop_length=HOP_LENGTH,
        n_mels=N_MELS,
        power=2.0,
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    waveform, sample_rate = torchaudio.load(librosa.ex("nutcracker"))
    waveform_cuda = waveform.cuda()

    # torchaudio (cpu)
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    start = time.time()
    melspec = mel_spectrogram(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    torch.cuda.synchronize()
    start = time.time()
    mel_spectrogram = mel_spectrogram.cuda()
    melspec = mel_spectrogram(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (CPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    start = time.time()
    logmel = feature_extractor(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (GPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    torch.cuda.synchronize()
    start = time.time()
    feature_extractor = feature_extractor.cuda()
    logmel = feature_extractor(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

実行結果

その1（3月24日版）

GPUを使ったtorchaudioが一番早い結果となった。

elapsed_time (librosa): 0.042655[sec]
elapsed_time (torchaudio; CPU): 0.021293[sec]
elapsed_time (torchaudio; GPU): 0.002648[sec]
elapsed_time (torchlibrosa; CPU): 0.237388[sec]
elapsed_time (torchlibrosa; GPU): 0.209118[sec]

その2（5月8日版）

koukyo1213様のコメントに従って、torch.nn.Sequentialやtorchaudio.transforms.MelSpectrogramのインスタンス化部分を計測対象から外した場合の結果を以下に示す。

elapsed_time (librosa): 0.040965[sec]
elapsed_time (torchaudio; CPU): 0.014024[sec]
elapsed_time (torchaudio; GPU): 0.001927[sec]
elapsed_time (torchlibrosa; CPU): 0.027673[sec]
elapsed_time (torchlibrosa; GPU): 0.006538[sec]

今回もGPUを使ったtorchaudioが一番早い結果となった。

考察とまとめ

最初のスクリプトでtorchlibrosaの実行時間がより多くかかったのは、koukyo1213様ご指摘の通り、インスタンス化のオーバーヘッドがかなり大きいことに起因すると考えられる。その影響を除いた新たなスクリプトでは、GPU使用のtorchlibrosaは、同じくGPU使用のtorchaudioと比較して3倍ほど低速という結果が得られた。とはいえ実行時間の桁は揃い、どちらも十分実用に供すると考える。

black用py-isortの設定

私はblackをpythonコードのformatterに用いているので、blackの仕様に合うようにisortを設定して動かしたいというわけである。isortをEmacsから触るためにはpy-isortを用いる。以下の記事は参考になる。

qiita.com

isortのマニュアルはこちらから： pycqa.github.io

まずisortをインストールする。

pip3 install isort

py-isortはEmacsのパッケージシステムから容易にインストール可能である。その後、Emacs設定ファイルに以下を追記する。こうするとバッファ保存時にisortが走ってくれるので便利である（＝つまりインポート順を自分で編集する必要がなくなる！）。

;; py-sort: isortを使ってライブラリのimport順をソート
(eval-when-compile (require 'py-isort))
;; for black
(setq py-isort-options '("--multi-line=3"
                         "--line-length=88"
                         "--trailing-comma"
                         "--use-parentheses"
                         "--ensure-newline-before-comments"
                         "--force-grid-wrap=0"))
;; バッファ保存時に自動ソート
(add-hook 'before-save-hook 'py-isort-before-save)

なおこれらのオプション値はblack公式のドキュメントを参考にした。 black.readthedocs.io

numpy.arange関数におけるstepに0.1を指定したときの振る舞いとその原因を追ってみた話

はじめに

先日、接点QB氏による以下のツイートを見かけた：

??????
俺の脳みそがバグってるのか？？ pic.twitter.com/zyzzdb3ijS

— 接点QB (@setten_QB) February 10, 2021

興味深いので少し調べてみたというわけである（少々長い記事）。

numpyのarange関数の振る舞い

numpyのarange関数のステップ数に0.1を指定したときの振る舞いについて、

np.arange(1.0, 1.5, 0.1)

とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4までのarrayが得られた。array要素の最大値1.4は引数stopの1.5より小さいというのが、想定される挙動である。

一方、

np.arange(1.0, 1.6, 0.1)

とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6のように、array要素の最大値が1.5を超えて1.6に達したので困惑したということである（「1.0, 1.1, 1.2, 1.3, 1.4, 1.5」かと思ったのに！）。

原因を探るため、numpyのマニュアルのarangeの項を見てみた。 numpy.org

すると、戻り値のarrayについて以下の記述が見つかった：

For floating point arguments, the length of the result is ceil((stop - start)/step).

つまり、浮動小数点数をarange関数の引数に与えた場合、array長がceil（＝小数点以下切り上げ）を伴う式で計算されるということだった。

手元で挙動を試してみると、以下のようになった。

こんな感じ？ pic.twitter.com/iJVVV91oUg

— mat (@ballforest) February 10, 2021

つまり、(1.5 - 1.0) / 0.1の計算結果が5.0となる一方で、(1.6 - 1.0) / 0.1の計算結果が6.000000000000001となり6.0にならない。

したがってceil関数を通したとき、

ceil((1.5 - 1.0) / 0.1) は5になるが、ceil((1.5 - 1.0) / 0.1) は7になった

というわけである。

問題は「なぜこのような計算結果が得られたのか」ということである。原因を一言で言えば「浮動小数点数演算に由来する計算誤差」ということに尽きるが、本記事ではこれを詳細に追ってみた（浮動小数点数の復習問題！）。

浮動小数点数の表現形式と演算

さて、浮動小数点数のIEEE754倍精度形式は次式で表される。 $$(-1)^{sgn} \times 2^{(exp - 1023)} \times (1+frac×2^{-52})$$ $sgn$は符号部、$exp$は指数部、$frac$は仮数部である（ここでは10進数表現）。

この表現形式を元にして、計算を追っていこうという方針である。 ここでは以下の記事を参考にした。 qiita.com

浮動小数点数の倍精度表現

検証にあたり、接点QB氏の具体例に基づいて、対象とする浮動小数点数を

0.1, 1.0, 1.5, 1.6

の4つに限定する。

0.1の倍精度表現

0.1の倍精度表現は

符号部 = 0
指数部 = 01111111011
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1019
仮数部 = 2702159776422298

となる。0.1の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1019 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{-4} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-56} \\ &= 7205759403792794 / 72057594037927936 \end{aligned} $$

1.0の倍精度表現

1.0の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 0000000000000000000000000000000000000000000000000000

であり、10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 0

となる。1.0の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 0 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \\ &= \frac{2^{52}}{2^{52}} \\ &= 4503599627370496 / 4503599627370496 \end{aligned} $$

1.5の倍精度表現

1.5の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 2251799813685248

となる。 1.5の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2251799813685248 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2251799813685248) \times 2^{-52} \\ &= (4503599627370496 + 2251799813685248) \times 2^{-52} \\ &= 6755399441055744 / 4503599627370496 \end{aligned} $$

1.6の倍精度表現

1.6の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 2702159776422298

となる。 1.6の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-52} \\ &= 7205759403792794 / 4503599627370496 \end{aligned} $$

倍精度表現に基づく浮動小数点数演算

以上で倍精度表現が揃ったので、検証に入ることができる（前置き長すぎワロタ）。 表に示しておく。

浮動小数点数	倍精度表現
0.1	7205759403792794 / 72057594037927936
1.0	1 (=4503599627370496 / 4503599627370496)
1.5	6755399441055744 / 4503599627370496
1.6	7205759403792794 / 4503599627370496

以上の表に基づいて、(1.5 - 1.0) / 0.1 を計算すると

((6755399441055744  - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936)
= 5

が得られた。ぴったし5である。

一方、(1.6 - 1.0) / 0.1 については

((7205759403792794 - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936)
= 6.000000000000001

と計算できた。

以上の計算結果に基づいてceilは小数点以下を切り上げるため、冒頭で示した挙動が説明できたわけである。

まとめ

本記事ではarange関数の引数のひとつであるstep数に浮動小数点数を適用した場合に、 戻り値のarrayの要素数がどのように変わりうるかを、具体的な数値例を元に調べてみた。

Numpyのマニュアルによれば、arangeによる戻り値のarrayの要素数はceil((stop - start)/step)で計算され、これらstart, stop, stepはarange関数の引数として与えられる。

接点QB氏のツイートをもとに、arangeの引数が具体的に

• (1) start=1.0, stop=1.5, step=0.1
• (2) start=1.0, stop=1.6, step=0.1

の場合に、浮動小数点数の倍精度表現に基づいて要素数を計算し、arange関数の振る舞いを検証した。 検証の結果、確かにPythonの計算結果は正しかったことが確かめられた。

Numpyマニュアルにも浮動小数点数が引数の場合はlinspace使ったほうがいいよ的な記述もあり、arangeを使うときには思い込みによるバグを埋め込まないように気をつけたい。

参考

• 2進数から10進数への変換 https://www.muryou-tools.com/sinsuu-change.php
• 関数電卓 https://tomari.org/main/java/dentaku_kansuu.html
• 浮動小数点数内部表現シミュレーター https://tools.m-bsys.com/calculators/ieee754.php
• 0.1は浮動小数点数で正確に表せないのに、printしたときに0.1と表示されるのはなぜか - Qiita
• Double-precision floating-point format - Wikipedia

vscodeのメモ

• ターミナルの表示およびカーソル切り替えはCommand + Shift + @
• ターミナルを閉じる場合はControl + Shift + @ （ターミナルにフォーカスされている場合）
• サイドバーを閉じる場合はControl + Option + Space (Toggle Sidebar Visibility)

わかりやすいパターン認識（第2版）の書評

本書は1998年に発売された『わかりやすいパターン認識』の改訂版であり，全体を通して統計的パターン認識の基礎が「わかりやすい」筆致でまとめられている。本書は全9章から構成されている。パターン認識の考え方（第1章）から始まり，識別関数の学習・評価・設計法（第2～4章），特徴抽出・次元削減（第5～7章）など，統計的パターン認識の分野における基礎的かつ重要な話題が精選されている。第8章では期待損失学習，第9章ではベイズ決定則と学習アルゴリズムの関係など，高度な話題にも触れられている。「統計的パターン認識はベイズ決定則を土台とする一学問」というフレーズには、改訂前から強い感銘を受けた。

改訂前後の大きな違いは，各章に数値計算・実験例や演習問題が数多く追加されていることである。例えば第2章には，パーセプトロンの学習過程が具体的な数値例とともに説明されており，第3章には3層ニューラルネットの実験例などが追加されている。演習問題では，具体的な数値例を与えて学習アルゴリズムの動作を確認させるものや，本文中の数式を具体的に導出させるものがある。詳細な解答例を載せたPDFが出版社のサポートページで公開されているのがありがたい。著者らの知見やノウハウを述べたCoffee Breakの欄にも加筆修正がされており，深層学習についても触れられている。「むすび」には，旧版以降に出版された参考書がコメントつきで紹介されており，初学者への良き指針となるだろう。

旧版は近年の深層学習・機械学習ブーム以前に発行されているが，本改訂で新たな章を割いて取り入れた話題はない。ただし深層学習を理解する上で重要な，「浅い」ニューラルネットであるパーセプトロンや３層ニューラルネット，誤差逆伝播法などに多くの紙面を割いて説明しており，本書発行の主旨を考えれば欠点にはあたらない。本書で取り上げられるトピックは技術の進展を経ても陳腐化しないものばかりが精選されているのである。

パターン認識に関心を持つ各位に本書の一読を広く薦めたい。