2021-02-24

Axial-Attentionについて調べたことのメモ

機械学習メモ

Axial-DeepLab (ECCV 2020, Spotlight、arXiv:2003.07853)の論文¹を読んでいて、Axial-Attentionが具体的にどんなものなのかよく分からなかったので、色々調べたり著者実装²を自分で実行してみたりして「こんな感じかなぁ」という認識に至ったので内容をメモしておく。

Attention？

そもそも「Attentionが何か」というのは、こちらの記事³が網羅的で分かりやすいと思う。

自分の理解では入力featureからqueryとkeyとvalueを生成し、queryとkeyの行列積で類似度表のようなものを作り、類似度表とvalueの行列積で値を取り出すと思っている。

類似度表がいわゆるAttentionなのだと認識している。（ただし、理解が正しいという自信はない）

Axial-Attentionモジュールの構造

論文より図を引用する。

$W_q$ 、 $W_k$ 、 $W_v$ はquery、key、value。r^q、r^k、r^vはquery、key、valueのpositional encodingベクトルとのこと。

処理内容を箇条書きすると以下のようになる。

入力featureを全結合層に通してquery、key、valueを得る
チャネル次元をグループに分割する
queryとkeyの行列積で類似度表をグループの数だけ生成する
query、keyにそれぞれpositional encodingベクトルを行列積して類似度表に加算する
類似度表とvalueの行列積、類似度表とpositional encodingベクトルの行列積を求めて合算する

※Fig.1中の「softmax」とある辺りが類似度表である

以下に補足する。

query、key、valueの生成、グループ分割

入力featureのshapeを[N, C, H, W]とする。（Nはミニバッチサイズ、Cはチャネル数、HとWのpixel数）

これを[N*W, C, H]に形状変更する。（Height-Axis版の場合、Width-Axis版はHとWが入れ替わる。以下同様）

さらに全結合層に通す。著者実装では全結合層は Conv2d 1x1 で実装されている。queryとkeyはチャネル数がAttentionモジュールの出力チャネル数の半分になる。さらにチャネル次元はグループに分割される。

query: [N*W, C, H] → Conv → [N*W, outs/2, H] → 分割 → [N*W, group, plane/2, H]
key:   [N*W, C, H] → Conv → [N*W, outs/2, H] → 分割 → [N*W, group, plane/2, H]
value: [N*W, C, H] → Conv → [N*W, outs,   H] → 分割 → [N*W, group, plane,   H]

out
sは出力チャネル数
groupはグループ数
planeはouts/group

類似度表の生成

query×keyで類似度表を生成する。

einsum('bgci, bgcj -> bgij', query, key)
[N*W, group, plane/2, H]×[N*W, group, plane/2, H] -> [N*W, group, H, H]

positional encodingベクトルがそれぞれquery、keyに適用されて類似度表に合算される。（要素ごとの和）

einsum('bgci, cij -> bgij', query, embedding)
[N*W, group, plane/2, H]×[plane/2, kernel, kernel] -> [N*W, group, kernel, H]
kernelはHと同じサイズ

kernelというのが、具体的に何を指しているのかは不明。著者実装を動かしてわかったことは、とにかくkernelサイズとHが一致する、ということだけだった。

類似度表はH次元でsoftmaxを取るので以下のようなH vs Hの表になる。

H vs Hのsimilarity表

これはある画素(i, j)の出力値を決める時に、W（行）方向は固定してH（列）方向の各pixelの出力を何％ずつ持ってくるか、という成分表のようなイメージになる。

上図を例に取ると、j=0の出力はH=0が10%、H=1が80%、H=2が10%でミックスした値となる。

なお、similarity表のshapeは[N*W, group, kernel, H]であるので、上図の表は実際には(N*W)×group数個存在している。（W方向はpixelごとに別の表になる）

類似度表から出力feature生成

valueとpositional encodingベクトルをそれぞれ類似度表と行列積して要素ごとの和を取る。

einsum('bgij, bgcj -> bgci', similarity, value)
[N*W, group, H, H]×[N*W, group, plane, H] -> [N*W, group, plane, H]
einsum('bgij, cij -> bgci', similarity, embedding)
[N*W, group, H, H]×[plane, kernel, kernel] -> [N*W, group, plane, H]

行列積の最内ループは内積なので H dot H になっていて、類似度表は0〜1の成分表なので成分表に従ってミックスする操作になっている。

下図は画像の中心点（だいだい色に塗ったpixel）の出力を決める時に縦方向のpixelを類似度表で参照している、という図である。 similarityベクトルのイメージ

類似度表との行列積が終われば、あとは形状変更で[N, C, H, W]にして出力となる。

なお、成分表がグループごとに別々なのでチャネル次元でみた時にグループごとに別々の領域にAttention（注目）してくれる可能性がある。（そのように学習が進めば）

positional encodingベクトルについて

自分はpositional encodingについて正直よくわかっていない。位置ごとに別の値を足した状態で正しく識別できるように学習するのだから、何らかの位置を加味した重みに学習されるに違いないくらいのふんわりした認識しかない。

学習対象のパラメータ

著者実装を見る限り、以下が学習パラメータ（重み）になりそうだった。

全結合層（Conv 1x1）の重み
positional encodingベクトル
Batch Normalizationのパラメータ

Axial-Attentionの（乱暴な）イメージ

自分の勝手なイメージだとBatch Normalizationは正則化項のような働きをしてくれる何かと認識しているので、実質Axial-Attentionは入力featureのどの部分を取り出すかを positional encodingベクトル、Conv 1x1 で決めているという理解であながち間違っていないのではないかと思う。

類似度表を直接重みにしていないのでとっつきにくいが、入力とConvとpositional encodingベクトルで類似度表が生成される、という点が理解できれば十分な気がした。

Axial-Attentionブロックの構造

Height-Axis版とWidth-Axis版のAxial-Attentionモジュールを組み合わせて縦横全域をカバーするAxial-Attentionブロックになる。

論文の図でMulti-Head、Concatと記載されている箇所は、チャネルをグループに分割してグループごとに別々の類似度表を適用していることに相当する。著者実装ではshapeの形状変更で対応できているのでConcatしていない。（要素ごとの和を取る際にConcatしているが、Multi-Headの件とは別の実装都合と思われる）

Attentionの可視化

論文よりAttentionの可視化画像を引用する。

最初にこの画像をみた際に、どこに着目して見ればよいのかが分からなかった。

column headとある図は、下図のように青いpixelの出力が縦方向に見て特に赤色になっているpixelの値を使っている、ということである。

row headとある図は、青いpixelの出力が横方向で赤い箇所の出力を多く使っているということになる。

以上。

Wang, Huiyu and Zhu, Yukun and Green, Bradley and Adam, Hartwig and Yuille, Alan and Chen, Liang-Chieh. Axial-DeepLab: Stand-Alone Axial-Attention for Panoptic Segmentation. European Conference on Computer Vision (ECCV). 2020↩
csrhddlam. 2020. Axial-DeepLab (ECCV 2020, Spotlight). https://github.com/csrhddlam/axial-deeplab , 2021/02/24閲覧↩
@halhorn. 2018. 作って理解する Transformer / Attention. https://qiita.com/halhorn/items/c91497522be27bde17ce , 2021/02/24閲覧↩

2020-07-04

ROCm3.5でPyTorchビルド

Python PyTorch ROCm

前回の続き？のような何か。

maminus.hatenadiary.org

今回やったこと

ROCm3.5 Docker環境
PyTorchをmasterブランチでビルド
torchvisionをmasterブランチでビルド
- 上記に伴い公式にバグ報告
前回同様PyTorch版EfficientDet¹で動作確認

前回からの変化点を中心にメモ。

ホスト

$ xhost +local:
$ docker pull rocm/pytorch:rocm3.5_ubuntu18.04_py3.6
$ docker run -it --name rocm35_pytorch -v ~/:/data -v /tmp/.X11-unix/:/tmp/.X11-unix/ --shm-size=4g --privileged --device /dev/video0 --device=/dev/kfd --device=/dev/dri --group-add video -e DISPLAY=$DISPLAY -e QT_X11_NO_MITSHM=1 rocm/pytorch:rocm3.5_ubuntu18.04_py3.6

Dockerのshmを4GB指定するオプションを追加。

コンテナ(PyTorch、torchvisionのインストール)

今回使ったROCmのdockerイメージは最初からAnaconda環境になっていて、virtualenvを使おうとするとビルドしたPyTorchがimportでエラーになったりするので、Anaconda環境にインストールするようにした。

あと、地味にユーザがjenkins（ID 1014）になっているので、ホストをマウントしてアクセスするのが面倒になった・・・

$ sudo apt-get update && sudo apt-get install kmod x11-apps

$ git clone https://github.com/pytorch/pytorch.git
$ cd pytorch/
$ git submodule sync
$ git submodule update --init --recursive
$ cd pytorch
$ .jenkins/pytorch/build.sh

$ git clone https://github.com/pytorch/vision.git
$ cd vision
$ vi setup.py
$ python setup.py install

ここは今回手順をごっそり入れ替えてソースコードからのビルド。前回はdockerにPyTorchが入っていたが、今回のバージョンは入っていないように見えた。

今回指定しなかったけど、export HCC_AMDGPU_TARGET=gfxXXXを指定したほうが良いと思う。

あと、torchvisionはpip installしようとするとPyTorch1.5.1が上書きインストールされるし、PyTorchだけビルド版を使うとtorchvision.nmsを呼び出すとエラーになるので、torchvisionもビルドする形にした。（ビルドしたらhipifyもやってくれてびっくり）

torchvisionの編集内容

libpng-configのオプションがUbuntuではdisableされているのでsetup.pyを書き換え。

diff --git a/setup.py b/setup.py
index 01276f1..e95137f 100644
--- a/setup.py
+++ b/setup.py
@@ -262,14 +262,21 @@ def get_extensions():
             png_version = parse_version(png_version)
             if png_version >= parse_version("1.6.0"):
                 print('Building torchvision with PNG image support')
-                png_lib = subprocess.run([libpng, '--libdir'],
+                help_result = subprocess.run([libpng, '--help'],
+                                         stdout=subprocess.PIPE)
+                lib_opt = '--libdir'
+                if '--L_opts' in help_result.stdout.decode('utf-8'):
+                    lib_opt = '--L_opts'
+                png_lib = subprocess.run([libpng, lib_opt],
                                          stdout=subprocess.PIPE)
                 png_include = subprocess.run([libpng, '--I_opts'],
                                              stdout=subprocess.PIPE)
                 png_include = png_include.stdout.strip().decode('utf-8')
                 _, png_include = png_include.split('-I')
                 print('libpng include path: {0}'.format(png_include))
-                image_library += [png_lib.stdout.strip().decode('utf-8')]
+                png_lib = png_lib.stdout.strip().decode('utf-8')
+                if png_lib:
+                    image_library += [png_lib]
                 image_include += [png_include]
                 image_link_flags.append('png')
             else:

~~本件、公式にIssue¹立てておきました。~~ 2020/07/07追記 masterに修正が入りました。

ここまでで、PyTorch・torchvisionはなんとなく動くはず。

コンテナ(EfficientDetの環境準備)

$ pip install opencv-python
$ cd EfficientDet.Pytorch
$ pip install -r requirements.txt
$ pip install git+https://github.com/cocodataset/cocoapi.git#subdirectory=PythonAPI

$ mkdir -p ~/.cache/torch/hub/checkpoints
$ cp efficientnet-b0-355c32eb.pth ~/.cache/torch/hub/checkpoints/

以上。

A Pytorch Implementation of EfficientDet Object Detection, https://github.com/toandaominh1997/EfficientDet.Pytorch, 2020/07/04閲覧↩
build error in ubuntu with libpng, https://github.com/pytorch/vision/issues/2392, 2020/07/04閲覧↩

2020-05-04

cv::fisheye::initUndistortRectifyMapのPについて

OpenCV メモ

自分用のメモ

OpenCV公式のIssue¹にある通り、cv::fisheye::initUndistortRectifyMap()の実装にバグがあるような気がする。x、yの値が0〜size.width、0〜size.heightになっているが、uvの計算式は-1.0〜+1.0？の範囲を想定しているように見える。

    for( int i = 0; i < size.height; ++i)
    {
// snip
        double _x = i*iR(0, 1) + iR(0, 2),
               _y = i*iR(1, 1) + iR(1, 2),
               _w = i*iR(2, 1) + iR(2, 2);

        for( int j = 0; j < size.width; ++j)
        {
// snip
                u = f[0]*x*scale + c[0];
                v = f[1]*y*scale + c[1];

2020/05/07追記

fisheyeモデルでは引数Pがデフォルトで単位行列なのが良くない？ような気もしてきた。
pinholeモデルの場合はgetDefaultNewCameraMatrix()を内部で呼んでる。fisheyeの場合、どのような値が望ましいのだろうか。とりあえずPにKをそのまま同じ値でわたすとそれらしい結果になった。

※バグと言い切れるか微妙なのでタイトルも変更

https://github.com/opencv/opencv/issues/5389#issuecomment-158017331 ↩

2020-04-16

PyTorch版EfficientDetをROCm環境で実行してみた

ROCm PyTorch

PyTorch版EfficientDet¹のカメラで撮影した画像に推論結果を重畳してくれるデモをROCm環境で動かしたので手順などをメモ。

環境

Radeon RX 580
Ubuntu 18.04.4 LTS（ホスト側）

dockerイメージ

Docker Hubでrocm/pytorchの中からそれっぽいのを選択。 Ubuntu18.04はROCm2.10が一番新しかった。

docker pull rocm/pytorch:rocm2.10_ubuntu18.04_py3.6_pytorch_profiling

ホスト側の準備[2020/05/12 追記]

以下のエラーが出てしまうのでxhostコマンドで許可を出しておく。

No protocol specified
: cannot connect to X server :0

ホスト側（コンテナの外）で実行する。

xhost local:

コンテナの起動コマンド

ホームディレクトリを/dataにマウント
X11アプリを動かすため/tmp/.X11-unixをコンテナに見せる
カメラをデモで使うので/dev/video0をコンテナに見せる
cv2.imshow()、cv2.waitKey()を動かすためDISPLAY変数、QT_X11 _NO_MITSHMを設定する
他の設定はコンテナ内でGPUを使うため

docker run -it --name pytroch_rocm -v ~/:/data -v /tmp/.X11-unix/:/tmp/.X11-unix/ --privileged --device /dev/video0 --device=/dev/kfd --device=/dev/dri --group-add video -e DISPLAY=$DISPLAY -e QT_X11_NO_MITSHM=1 rocm/pytorch:rocm2.10_ubuntu18.04_py3.6_pytorch_profiling

※1：コンテナ起動時にカメラは接続済みの状態にしておく
※2：QT_X11 _NO_MITSHMが無いとcv2.waitKey()で以下のエラーになる

X Error: BadAccess (attempt to access private resource denied) 10
  Extension:    130 (MIT-SHM)
  Minor opcode: 1 (X_ShmAttach)
  Resource id:  0x4eb

Docker上で追加作業

X11の準備

デモの中でOpenCVのウィンドウを表示するのでx11-appsを追加した。

apt-get update
apt-get install x11-apps

Pythonパッケージ

pip installでエラーになるケースがあるので、いつもpipとsetuptoolsを最新にしている。システムのpip3を更新すると環境を壊してしまうことがあるのでvenv環境を作ってそちらに入れている。

apt-get install python3-venv
python3 -m venv --system-site-packages env_rocm
source env_rocm/bin/activate
pip install --upgrade pip
pip install --upgrade setuptools

EfficientDetの依存パッケージ

デモでOpenCVを使っているのと、COCO APIを使っているので一緒にインストールする。

pip install -r requirements.txt
pip install opencv-python
pip install git+https://github.com/cocodataset/cocoapi.git#subdirectory=PythonAPI

EfficientNetの重みデータをキャッシュディレクトリに置く

EfficientDetのGitHubサイトにはEfficintDet部分の重みデータへのリンクが貼られているが、どうやらバックボーンネットワークのEfficientNetの重みを別途実行時にダウンロードするらしい。

現状ソースコードに埋め込まれているダウンロードURLがリンク切れしているようなので、EfficientNetのPyTorch実装GitHub²サイトからダウンロードして、キャッシュディレクトリに置く。

# wgetなどでダウンロードしておく

mkdir -p /root/.cache/torch/checkpoints
cd /root/.cache/torch/checkpoints
ln -s /data/python3.6/efficientnet-b0-355c32eb.pth efficientnet-b0-355c32eb.pth

ソースコード修正

dockerイメージに含まれているtorchvisionはCUDA実装がビルドされていないようでNMS呼び出しで以下のエラーが出る。

  File "demo.py", line 177, in <module>
    detect.camera()
  File "demo.py", line 155, in camera
    show_image = self.process(img=img)
  File "demo.py", line 81, in process
    scores, classification, transformed_anchors = self.model(img)
  File "/root/.local/lib/python3.6/site-packages/torch/nn/modules/module.py", line 539, in __call__
    result = self.forward(*input, **kwargs)
  File "/data/python3.6/EfficientDet.Pytorch/models/efficientdet.py", line 83, in forward
    transformed_anchors[0, :, :], scores[0, :, 0], iou_threshold=self.iou_threshold)
  File "/root/.local/lib/python3.6/site-packages/torchvision/ops/boxes.py", line 31, in nms
    return torch.ops.torchvision.nms(boxes, scores, iou_threshold)
RuntimeError: Not compiled with GPU support (nms at /tmp/pip-wnkd089r-build/torchvision/csrc/nms.h:20)

torchvisionのHIPIFYのやり方がよくわかってないのでビルドは諦めて呼び出し側でtorch.tensorをいったんCPUに戻すようにしたらエラー回避できた。

diff --git a/models/efficientdet.py b/models/efficientdet.py
index 43357c9..83d60d9 100644
--- a/models/efficientdet.py
+++ b/models/efficientdet.py
@@ -80,7 +80,7 @@ class EfficientDet(nn.Module):
             transformed_anchors = transformed_anchors[:, scores_over_thresh, :]
             scores = scores[:, scores_over_thresh, :]
             anchors_nms_idx = nms(
-                transformed_anchors[0, :, :], scores[0, :, 0], iou_threshold=se
lf.iou_threshold)
+                transformed_anchors[0, :, :].cpu(), scores[0, :, 0].cpu(), iou_threshold=self.iou_threshold).cuda()
             nms_scores, nms_class = classification[0, anchors_nms_idx, :].max(
                 dim=1)
             return [nms_scores, nms_class, transformed_anchors[0, anchors_nms_idx, :]]

実行

実行は作者様のGitHubサイトの通りに実行。EfficientDetの重みデータは、あらかじめGitHubに貼られているリンクからダウンロードしておく。

python demo.py --threshold 0.5 --iou_threshold 0.5 --cam --score --weight ../checkpoint_VOC_efficientdet-d0_268.pth

A Pytorch Implementation of EfficientDet Object Detection, https://github.com/toandaominh1997/EfficientDet.Pytorch ↩
EfficientNet PyTorch, https://github.com/lukemelas/EfficientNet-PyTorch/releases ↩

2020-03-12

pip installできるようにしました

Python ONNX

前回作ったなんちゃってライブラリを改良？したので報告。 maminus.hatenadiary.org

前回はライブラリをpip install .すらできなかった。

こんなエラーになる。

$ pip install .
Processing /home/maminus/python3.6/ommle
    ERROR: Command errored out with exit status 1:

<snip>

File "/tmp/pip-req-build-n0by817_/ommle/__init__.py", line 1, in <module>
    from extractor import add_middle_outputs
ModuleNotFoundError: No module named 'extractor'

__init__.pyからextractor.pyを指定してるところでモジュールが見つからないと怒られている。

ディレクトリ構成はこんな感じ

ommle
  +-- ommle
         +-- __init__.py
         +-- extractor.py

パッケージ名から指定するとエラーにならないし、pip install git+...もできた。

-from extractor import add_middle_outputs
+from ommle.extractor import add_middle_outputs

よく考えたらfrom .extractor import add_middle_outputsでもいけるかもしれない気がする。

いまだにimportの挙動が難しすぎて理解できてない・・・

2020-02-29

ONNXモデルの中間ノード出力を抽出するライブラリを作ってみた

Python ONNX

先日記事にしたこれを実装してみた。 maminus.hatenadiary.org

ソースコード

GitHubに上げました。まともな？ライブラリを作るの初めてなので何か間違ってるかも・・・（特にPipenv自信なし。pip install も直接gitを指定できなかった。たぶんどこか間違ってる）

やってること

ONNXモデル（onnx.ModelProto）を引数に受け取り、編集後のONNXモデルを戻り値で返す関数 add_middle_output() を実装
各ノードのoutput[]からIdentityノードを生やす
model.graph.outputに中間層出力用のValueInfoProtoを追加
- Identityノードの出力がValueInfoProto
- 出力名はOMMLE.node_output_name.middle_0の形式
  - node_output_nameがノードの出力名
  - 上記の名前で既存の名前と衝突する場合は最後の0をインクリメントする
引数cast_typeにonnx.TensorProtoのデータ型を指定した場合、Identityノードの代わりにCastノードを使う
- 指定したデータ型にキャストした中間ノード出力が得られる
- 事前に中間ノード出力のデータ型が不明な場合に使う想定
- cast_typeを省略した場合、onnx.TensorProto.FLOAT（numpy.float32相当）である前提で動くので注意
引数exclude_op_types
リストにノードのop_type名を入れると該当するノードから中間出力を出さない
- op_type名は完全一致で比較する
引数exclude_output_namesリストにノードのoutput名を入れると該当するoutput名から中間出力を出さない
- 同じく完全一致で比較
- output[]名はNetronなり何なりで調べて指定する

メモ

Identityノードを挟まずに直接ノードの出力をONNX出力に指定すると推論時に以下のようにエラーになる。（ONNX仕様なのか、ONNXRuntineの実装都合なのかは未調査）

>       self._sess.load_model(providers)
E       onnxruntime.capi.onnxruntime_pybind11_state.Fail: [ONNXRuntimeError] : 1 : FAIL : This is an invalid model. Graph output (conv_1) does not exist in the graph.

2020-02-28

浮動小数点誤差の自分メモ

メモ

普段適当にfloat32やfloat16使ってるけど、適当すぎてたまに困るので勉強のために机上での考察とプログラムを動かして検証してみた。

オーバーフロー、アンダーフロー

指数部が表現可能な範囲外を超えるケース
float32なら指数が+127を超えるとオーバーフローしてinfになる
下限の-126を下回ると非正規化数
- 正規化数は ±1.xxxx × 2^yyy で表現している
- 非正規化数は ±0.xxxx × 2^-126
  - 一番小さな数値になるのは ±0.0000...0001 × 2^-126
  - float32の場合、仮数部は23bitなので 0.00...01 は 2^-23
  - 仮数部 2^-23 × 指数部 2^-126 で ±2^-149 が非正規化数の最小値
  - 上記を一般化すると ±2^（仮数部bit数＋指数部最小値）が最小値
演算結果が非正規化数の最小値より小さな値になるとアンダーフローして0になる

確認コード

    // 16進数浮動小数リテラルはC++17以降で使える
    std::cout << 0x1p+127f * 0x1p+127f << std::endl;     // オーバーフロー（生の数）
    std::cout << -0x1p+127f * 0x1p+127f << std::endl;    // オーバーフロー（負の数）
    std::cout << 0x1p-126f * 0x1p-25f << std::endl;      // アンダーフロー

実行結果

inf
-inf
0

float16

float16は指数部が5bitしかない
- 指数部は-14〜+15
- 15乗と聞くと大きな数値に聞こえそうだが２進数なので2¹⁵==32768
- 仮数部がall 1でも65504が最大値
- 普通にuint16_tの最大値の方が大きい
- 255×255はセーフだが、256×256はオーバーフローする
正規化数の最小値は ±2^-14
- こちらも 1/16384 と正直小さな値と言い難い
非正規化数の最小値は ±2^-24
- 指数部2^-14 仮数部 2^-10
- 1/16777216。微妙な数値。うっかりアンダーフローしそう
以前float16でディープラーニングモデルを学習させようとしたら速攻でNaNになって学習できなかった
- おそらくアンダーフローが発生して0になり0除算でNaNという流れ

実験用ルーチン

NaN、INFなどの特殊値や非正規化数、なんなら±0にすら非対応のなんちゃって実験用ルーチンのコード。文字列のオレオレバイナリ表現とfloat値の相互変換ができる。（実用には耐えられない）

#include <cstdint>
#include <string>
#include <algorithm>
#include <iostream>

std::string remove_char(const std::string &str, char c)
{
        std::string result = str;
        auto new_end = std::remove(std::begin(result), std::end(result), c);
        result.erase(new_end, std::end(result));
        return result;
}
std::string cleansing(const std::string &str)
{
        auto s = remove_char(str, ' ');  // 半角空白を削除
        return remove_char(s, ',');  // カンマを削除
}
float str_to_fp32(const std::string &str)  // なんちゃってバイナリ表現文字列 → float変換
{
        auto s = cleansing(str);  // 区切り文字と空白を削除
        std::uint32_t raw = 0;     // バイナリ表現

        // 符号ビット
        if (s[0] == '-') {
                raw = 0x80000000;
        }

        // 仮数部
        std::uint32_t mantissa = 0;
        for (auto i=0;i<23;i++) {
                auto bit = (s[3 + i] == '1')? 1: 0;
                mantissa <<= 1;
                mantissa |= bit;
        }

        // 指数部
        auto exponent = std::stoi( std::string(std::begin(s) + 3 + 23 + 4, std::end(s)) );
        exponent += 127;

        raw += exponent << 23;
        raw += mantissa;

        return *reinterpret_cast<float *>(&raw);
}
std::string fp32_to_str(float f)  // float → バイナリ表現文字列変換
{
        std::uint32_t v = *reinterpret_cast<std::uint32_t *>(&f);
        std::string result;

        result = (f > 0)? "+": "-";
        result += "1.";

        // 仮数部
        std::uint32_t mantissa = v & 0x7fffff;
        for (auto i=0;i<23;i++) {
                auto bit = (v >> (22 -i)) & 0x01;
                result += (bit == 1)? "1": "0";
                if((i % 4) == 2 && i!=22) {
                        result += ",";
                }
        }

        // 指数部
        auto exponent = (v >> 23) & 0xff;
        result += " x 2 ** " + std::to_string(static_cast<int>(exponent) - 127);
        return result;
}

情報落ち

加減算する時に２数の差が大きいと計算結果が変化しない
- 2^{仮数部のbit数} だけ差があると情報落ちする
計算結果は指数部の大きい方になる
例えば 0.5f は 1.000...×2^-1 で指数部-1なので仮数部は 2^-2 〜 2^-24 の 0/1 を表現できる
- 仮数部の分解能？が2^-24という解釈ができそう
  - 仮数部が 2^-24 の倍数と解釈できる（仮数部の２進数×2^-24）
  - 上記の解釈を受けてこの記事では仮に分解能と呼称する
- 分解能が2^-24なので、2^-25以下の数値は表現できない
- よって2^-25以下の数値を加算しても結果が変化しない
指数部の数値によって分解能は変化する
- 2.0f は 1.000...×2¹ 指数部1 で仮数部は 2⁰〜2^-22
- 1.0f は 2^-1〜2^-23
- 0.25f は 2^-3〜2^-25
指数部が小さい(-126に近づく方向)ほど仮数部の分解能が小さくなる
- 指数部-126 なら 2^-127〜2^-149
- イメージとしては指数部の値によって物差しの大きさが変化するイメージ
  - 仮数部が23bitの固定なので物差しの目盛り数が23個の固定イメージ
  - この記事で分解能と呼んでいるものが物差し自体の大きさ
- 仮数部目盛り値×2^{指数部−仮数部bit数}
- 指数部が小さいと分解能が上がる件を解説してくれている書籍¹
  お試しコード

// 情報落ちしないギリギリのライン
// fp32の場合、仮数部23bitなので2^23までの差は情報落ちしない
std::cout << "1.0 + 2**-23 -> " << fp32_to_str(1.f + std::pow(2, -23)) << std::endl;
std::cout << "0.5 + 2**-24 -> " << fp32_to_str(.5f + std::pow(2, -24)) << std::endl;

// 情報落ちする。2^23を超える差がある
std::cout << "1.0 + 2**-24 -> " << fp32_to_str(1.f + std::pow(2, -24)) << std::endl;
std::cout << "0.5 + 2**-25 -> " << fp32_to_str(.5f + std::pow(2, -25)) << std::endl;

実行結果

1.0 + 2**-23 -> +1.000,0000,0000,0000,0000,0001 x 2 ** 0
0.5 + 2**-24 -> +1.000,0000,0000,0000,0000,0001 x 2 ** -1

1.0 + 2**-24 -> +1.000,0000,0000,0000,0000,0000 x 2 ** 0
0.5 + 2**-25 -> +1.000,0000,0000,0000,0000,0000 x 2 ** -1

情報落ちと桁落ちのコンボ

桁落ちそのものは「有効数字減って困るんだっけ？」くらいの認識だった
情報落ちと桁落ちが組み合わさって誤差を観測するケースを体験したので共有しておく
①で小さな値(指数部-8)と大きな値(指数部0)を加算する
- 情報落ちによって小さな値の仮数部後半が消失する
➁で先程加算した値を減算する
- 理想状態なら元の小さな値に戻る
- 実際には桁落ちが発生して仮数部後半がなくなる
  - 感覚的にはむしろ最初の加算で仮数部後半が消失しているので元に戻らないという感じ

①
+1.000,1000,1000,1000,1000,1000 x 2 ** -8
 + 
+1.000,0000,0000,0000,0000,0000 x 2 ** 0

結果１
+1.000,0000,1000,1000,1000,1001 x 2 ** 0

➁
+1.000,0000,1000,1000,1000,1001 x 2 ** 0
 - 
+1.000,0000,0000,0000,0000,0000 x 2 ** 0

結果２
+1.000,1000,1000,1001,0000,0000 x 2 ** -8

乗算の誤差

個人的には指数部は単純な指数同士の加減算イメージ
仮数部も乗算
- この時 1.00...001 みたいに1で囲んだbit数（最後の1が出現するまでのbit数）は足し算で増える
- 例えば 1.001 (3bit) × 1.00001 (5bit) は 1.00101001 (3+5=8bit)
- 上記の数値が仮数部bit数を超えると仮数部で表現しきれず誤差が発生する

除算の誤差

正直はっきりした「この条件で誤差が出る」というイメージがない（自分が知らないだけかも）
4bit相当の除算を試した感じだと割り切れないケースで仮数部を超えてしまう(丸め誤差発生)という印象

お試しコード

for (auto i=3;i<15;i++) {
    auto lhs = 15.f;
    std::cout << "15 / " << i << " -> " << fp32_to_str(lhs / (i * 16)) << std::endl;
}

実行結果

15 / 3 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -2
15 / 4 -> +1.111,0000,0000,0000,0000,0000 x 2 ** -3
15 / 5 -> +1.100,0000,0000,0000,0000,0000 x 2 ** -3
15 / 6 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -3
15 / 7 -> +1.000,1001,0010,0100,1001,0010 x 2 ** -3
15 / 8 -> +1.111,0000,0000,0000,0000,0000 x 2 ** -4
15 / 9 -> +1.101,0101,0101,0101,0101,0101 x 2 ** -4
15 / 10 -> +1.100,0000,0000,0000,0000,0000 x 2 ** -4
15 / 11 -> +1.010,1110,1000,1011,1010,0011 x 2 ** -4
15 / 12 -> +1.010,0000,0000,0000,0000,0000 x 2 ** -4
15 / 13 -> +1.001,0011,1011,0001,0011,1011 x 2 ** -4
15 / 14 -> +1.000,1001,0010,0100,1001,0010 x 2 ** -4

平山尚　（2008）　ゲームプログラマになる前に覚えておきたい技術　秀和システム　pp.745-746．↩