Author Archives: hamukazu

数値計算で必要な数式変形について

@beam2dさんの記事が面白かったのだが、これを読んでよくわからなかった人もそれなりにいるだろうという想定で、補足説明のつもりでこれを書いている。

いろんなところで話をするのだが、「数学≠数値計算」である。数学の知識と数値計算の知識はちょっと違う。方程式の解が解析的に解けて、やった解けた!と思っても、その値をコンピュータで計算して結果を出すまでにさらに工夫が必要なことがよくある。コンピュータ内部での数値表現の特性を理解した上で、コンピュータフレンドリーな数式にさらに変形しなければいけないことがある。

例1: Sigmoid関数の計算

機械学習でよく出てくるSigmoid関数
$$ \sigma(x) = \frac{1}{1+e^{-x}} $$
であるが、これは\( \lim_{x\to \infty} \sigma(x) =1 \), \( \lim_{x\to -\infty} \sigma(x) =0 \)という性質がある。これは、この式のまま実装すると、\( x \)が小さい時にうまくいかない。

実際にPythonでやってみる。

>>> def sigmoid_(x):
...     return 1/(1+np.exp(-x))
... 
>>> sigmoid_(10)
0.99995460213129761
>>> sigmoid_(-10)
4.5397868702434395e-05
>>> sigmoid_(-1000)
__main__:2: RuntimeWarning: overflow encountered in exp
0.0

これは、np.expのカッコの中の絶対値が大きすぎてオーバーフローしてしまったことが原因である。しかし、expの計算ができなくてもsigmoid関数は0〜1の値なので計算できないとおかしい。

これを解決するにはいろんな方法が考えられるが、最も単純には、\(x\)がある値\(a\)(ただし、aは負で、exp(-\(a\))がオーバーフローしないようなぎりぎりの値)より小さいときは0にしてしまえばよい。つまり次のように定義すればよい。
$$
\sigma(x)=
\left\{
\begin{array}{ll}
0 & (x \lt a)\\
\frac{1}{1+e^{-x}} & (x\geq a)
\end{array}
\right.
$$

(ただし、\(a\)はexp(-\(a\))がオーバーフローしない程度に十分小さい値)

>>> def sigmoid(x):
...     return 0.0 if x<-709 else 1/(1+np.exp(-x))
... 
>>> sigmoid(10)
0.99995460213129761
>>> sigmoid(-10)
4.5397868702434395e-05
>>> sigmoid(-1000)
0.0

ここで、\(a\)として-709を使った。実際sigmoid(-709)を計算すると結果は1.2167807506234229e-308となり、これはほぼ0である。

例2: Softplus関数の計算

Softplus関数
$$ \tau(x)= \log (1+e^x) $$
の計算をしたいとする。これもこの式のまま計算すると\(x\)が大きい時にexpがオーバーフローするが、一方で関数の真の値はほぼ\(x\)に近い値になる。

ここで@beam2dさんのブログに書かれていたとおり、
\[
\begin{align}
\tau(x)&= \left\{
\begin{array}{ll}
x+\log (1+e^{-x}) & (x > 0) \\
\log (1+e^x) & (x \leq 0)
\end{array}\right.\\
&=\max(0,x)+\log (1+e^{-|x|})
\end{align}
\]

と変形してから実装すればいい。実際、式変形する前と後での動作の違いを見てみるとわかるかと思う。

>>> def softplus_(x):
...     return np.log(1+np.exp(x))
... 
>>> softplus_(1000)
inf
>>> def softplus(x):
...     return np.maximum(0,x)+np.log(1+np.exp(-np.abs(-x)))
... 
>>> softplus(1000)
1000.0

ここでsoftplus_は元の定義式をそのまま使ったものだが、これに1000を引数に与えるとinfになってしまう。softplusは式変形したものを利用したもので、1000を与えると1000が返ってくる。数学的にはこっちのほうが正しい。

おまけ: NumPy芸

上記のsoftplusの実装は、すべてユニバーサル関数なので、引数が配列のときも有効である。実際やってみると次のようになる。

>>> softplus(np.array([10,100,1000]))
array([   10.0000454,   100.       ,  1000.       ])

上記のsigmoid関数の方は配列にそのまま適用できないが、numpy.whereを使って書き換えれば簡単にユニバーサル関数化できる(読者の演習問題とする)。

まとめ

* 数学的考察で解析的な解が求まったとしても、そこからコンピュータで数値的に計算できるかまでにはギャップがあることがある。
* 関数値が有限でも、計算過程でとても大きな値が発生するケースが要注意で、その場合に数式の同値変形で回避できることがある。
* Pythonで大量に同じ関数を計算することを考えると、場合分けも含めてできるだけユニバーサル関数で書いておいたほうがいい。なにしろ、PythonではFor文書いたら負けなので。

文献案内

数値計算に興味をもったら、是非これを読むことをお勧めする。


数値計算の常識

古い本だが内容は今でも十分に役に立つ。桁落ちのワナなどの説明がくわしい。ページ数も少ないし読みやすいと思う。

更新履歴:
2015/8/1 softplusが途中からsoftmaxになってたのを直しました。

Python関連のEmacs設定について補足

以前書いた「PythonプログラミングのためのEmacs設定(動画あり)」だが、情報が古くなり、いつのバージョンからかうまく動かなくなってしまった。
(以前のブログは最新情報に合わせて書き換えた)

問題なのはpy-autopep8の設定なのだが、

(add-hook 'before-save-hook 'py-autopep8-before-save)

(add-hook 'python-mode-hook 'py-autopep8-enable-on-save)

に変えなければいけない。

ここで面白いのは、以前はbefore-save-hookへの設定だったものが、python-mode-hookに変わったこと。つまりファイル保存時の起動をそとから明示的にキックしてもらうのではなく、自発的にやるようになった。

最新版ではpy-autopep8-before-saveは消えてしまったので、旧設定のままpy-autopep8のアップグレードをして、pythonファイルを開くとエラーで保存ができなくなってしまいハマった。ということを再現しようと思ったら、いまは再現せず保存できるようになった。Emacs本体側で対応したのかかもしれない。問題に気づいてからしばらく経っているうちに状況が変わってしまったようだ。

ところで、Windowsでpy-autopep8を使うと保存するたびに一番左上にカーソルが動いてすごく不便なんだけど、なんとかならないだろうか。

Remark.jsにMathJaxの数式を入れる方法

色々とハマりつつも、Remark.jsとMathJaxを使ってプレゼン用スライドを作るのがとても便利なことに気づいたので紹介する。

HTML+JavaScriptでプレゼン用スライドを作る仕組みはいくつかあって、色々と試してみたが、最近はRemark.jsが便利だと思っている。Markdown形式を自動で変換してくれるというのと、配布資料用として印刷してもちゃんとレイアウトしてくれるところがよい。あと、デザインも変に凝らずにシンプルなところがよい。といってもRemark.jsでも、CSSをいじれば凝ったことはいくらでもできるようだが私はやらなない。

仕事柄どうしてもスライドに数式が多くなるので、その点MarkdownでTeX形式の数式を挿入できるのがとても嬉しいと思っていたのだが、色々といじってるとハマってしまった。以下それを説明する。

途中をすっ飛ばしてとりあえずRemark.js+MathJaxを使いたいという人は、最後の「まとめ」だけ読めばよい。

本家のページからテンプレートをもらってきて、次の数式を入力してみる。

$$\sum_{a_i\in A_j}f(a_i)$$

するとちゃんと、
$$\sum_{a_i\in A_j}f(a_i)$$
と表示される。

やった、これで便利! 数式も簡単に入れられる! と思って、

$$\sum_{i=1}^n \lVert a_i \Vert_\infty$$

と入れるとなぜがうまく動かない。実際のHTMLを見てみると、「_(アンダースコア)」のペアがMarkdown文法の「強調」の意味に取られて<em> , </em> に変換されてしまうのが問題らしい。前者の例でもアンダースコア2つを使ってるのだがなぜか問題なく、正確なルールはよくわかっていない。

実際にうまくいっていない様子は、こちらを参照してほしい。

ちなみに2つ目の例は、正しくレンダリングされると
$$\sum_{i=1}^n \lVert a_i \Vert_\infty$$
となるはずである。

で、解決策だが、こちらを参考にした。

具体的には、MathJaxのパラメータに「delayStartupUntil=configured」というのを加え、以下のコードを加えればよい。

      MathJax.Hub.Config({
          tex2jax: {
          skipTags: ['script', 'noscript', 'style', 'textarea', 'pre']
          }
      });
      MathJax.Hub.Queue(function() {
          $(MathJax.Hub.getAllJax()).map(function(index, elem) {
              return(elem.SourceElement());
          }).parent().addClass('has-jax');
      });

      MathJax.Hub.Configured();

具体的に何をやっているかまで理解できていないのだが、ここで呼んでいるのはすべてMathJaxのAPIなので、JavaScriptを使った他のシステムでも有効な方法なのかもしれない。

まとめ

以上のことを全部含めた自分用のテンプレートを作ったので公開する。ソースはここにおいてある。

使い方は以下のとおり:

  • レポジトリからindex.htmlをダウンロードする
  • index.htmlと同じディレクト内のCONTENT.mdにMarkdownでスライドのテキストを書く
  • Webサーバを立ち上げる(Python2なら「python -m SimpleHTTPServer」、Python3なら「python -m http.server」とするのが便利)
  • ブラウザでindex.htmlを開く(サーバでPythonを使った場合は「http://localhost:8000」を開く)

そして、ここで数式の記述にも注意が必要で、バッククォートでくくってやる必要がある。最初の例で言うと次のようにする:

`$$\sum_{a_i\in A_j}f(a_i)$$`

また、インラインの数式の場合も同様で、例えば次のようにする:

`\(\sum_{a_i\in A_j}f(a_i)\)`

サーバが必要なのは、index.htmlからCONTENT.mdをインクルードしてるからで、html内の<textarea id=”source”>~</textarea>にMarkdownを直接書けばサーバは必要ない。また、CONTENT.md以外のファイルをインクルードしたければindex.html?<filename>とすれば動くようになっている。

CONTENT.mdのサンプルはレポジトリに含まれていて、ここにそのデモを用意しておいた。

しばらくはこれで、いわゆるコンサル風のプレゼンとは違う方向に突っ走っていきたい。

追記:

  • 数式の記述でバッククォートを入れるという説明を忘れたので加筆した
  • Pythonを使ったHTTPサーバの立ち上げ方を加筆した

コードが書けなくてもOSSに貢献できる3つの方法

世の中オープンソース・ソフトウェア(以下「OSS」)にあふれているし、プログラマじゃなくても、コードが書けなくても、例えばLinuxの上で動いているMySQLの上で動いているWordpressでホームページをデザインしている人は多いはずだ(LinuxもMySQLもWordpressもOSSである)。

OSSはソースコードが公開されていてだれでも自由に改変できることが特徴なのだが、実際にコードを改変する人は少ないし、改変したコードで元のソフトウェアに貢献している人はもっと少ない。だが、OSSに貢献できるのは一部の技術者だけだと思っていたらそれは間違いだし、コードを書けなくても貢献する方法がいくつかある。それを知ってもらおうと思いこの記事を書いた。といっても、恩恵を受けてるのだから必ず貢献しろというつもりは毛頭なく、もしできることがあればやってほしいと思う。自分にはなんの貢献もできないという考えを改めてほしい。

以下貢献の方法を、安易な順に挙げる。

寄付をする

これが一番簡単で時間のかからない貢献のしかたである。例えば、Free Software Foundationのトップページには「donate」というリンクがあるし、Ubuntuをダウンロードしようとする寄付しろと言ってくる(もちろん寄付しなくてもダウンロードはできる)。あるいは、Tシャツやマグカップなどのグッズを買うことが実質的な寄付になっていることもある。

私は少額であるが、FSFにもUbuntuにも寄付している。

バグを報告する

バグを見つけたときに放置せずにしかるべきところに報告するのは、OSSの発展に貢献することになる。また、報告することで自分が困っているバグの修正が早くなる可能性があり、直接のメリットを受けるかもしれない。これは場合によっては英語力が必要とされるかもしれないが、日本語だけで対応できることもありうる。使っているOSSにバグを見つけた場合、まずは(それが意外と難しかったりするが)どこに報告すべきかを見つけ、できるだけ報告しよう。

ある動作がバグが仕様かというのは悩ましいケースもあるが、「仕様としては不自然で、ちょっと調べてみたところどこにも仕様という記載はない」というケースはまず報告していいと思う。バグを報告するのに何時間もドキュメントを読んで確認する必要はないだろう。そのかわり、「仕様です。却下。」と言われても腹を立てたり落ち込んだりしない心が必要。

また、報告したバグが長期間放置されることもあるかもしれないが、それも腹を立てたり落ち込んだりしない心が必要。

私はUbuntuのバグを見つけたらほぼすべて報告するようにしている。

翻訳をする

もしあなたが英語ができて、日本語ドキュメントのないOSSを使っているのなら、それを日本語に翻訳するという貢献の仕方もある。ドキュメント以外にも、メニューやチュートリアルなど、訳してもらうと喜ばれることはいろいろある。もちろん、訳し始めたからには必ず全部やる必要があるわけではなく、途中まででもそれを公開することには価値がある。勉強会などのコミュニティで翻訳するというのもよくある話で、そういう場で翻訳プロジェクトを提案するというのも貢献である。英語と技術を同時に勉強できるということで、勉強としてやるのもいいかもしれない。

ただ、翻訳についてひとつ気をつけて欲しいのは、ライセンスについて。OSSでよく使われるGPL, BSD, MITなどのライセンスについては、勝手に訳して勝手に公開しても問題ないが、一般にはいくら公開情報でも勝手に翻訳して公開することは法的に問題があることがある。ちゃんと確認してからやること。

かなり昔のことではあるが、私はboost(というC++のライブラリ)の翻訳プロジェクトに関わった。

まとめ

  • コードが書けなくてもOSSへの貢献はできる。できるところからやって欲しい。
  • 全部実践している自分偉い

2015年の素数日

今年に入って5日が過ぎてしまったが、ここで毎年恒例の素数日リスト公開。

今年の素数日は以下のとおり。

20150111
20150131
20150227
20150303
20150327
20150411
20150513
20150821
20151011
20151031
20151127
20151221
20151227

最初の素数日に間に合ってよかった。
ソースコードはこちら参照。

PythonプログラミングのためのEmacs設定(動画あり)

EmacsのPython環境周りを最近整えたので自分のためのメモ。似たような情報は色々とネットに転がっているのだが、古いEmacs向けだったり、いろいろとやりすぎだったりというのもあったので、自分なりにシンプルにまとめてみた。

以下、Ubuntuでemacs24を使ってることを前提とする(Linuxならばそんなに変わらないかなと思うが)。python-modeとpackage.elはすでに入っているとする。python-pipがインストールされてて、~/.local/binにパスが通っていることも前提とする。

以下の順で説明する。

  • 準備:package.elの設定
  • 自動インデントの設定
  • yasnippetのインストール
  • py-autopep8のインストール
  • auto-completeのインストール
  • elpyインストール
  • flymake関連のインストール
  • 動作デモ(動画)

それぞれ独立なので、機能を理解して自分がほしい物を入れればいいと思う。

準備:package.elの設定

まず前提としてpackage.elで使うため、~/.emacs.d/init.elに以下を加える。

(require 'package)
(setq package-archives
      (append '(("marmalade" . "http://marmalade-repo.org/packages/")
                ("melpa" . "http://melpa.milkbox.net/packages/"))
              package-archives))
(package-initialize)

自動インデントの設定

python-modeのときにEnterキー(またはC-h)を押すと、自動でいい感じにインデントしてくれるように設定する。

.emacs.d/init.elに以下の行を挿入する。

(add-hook 'python-mode-hook
          (lambda ()
            (define-key python-mode-map (kbd "\C-m") 'newline-and-indent)
            (define-key python-mode-map (kbd "RET") 'newline-and-indent)))

自動インデントにすると、余計な空白が入ったままになることがよくあるが、後述するpy-autopep8で保存時に変換してくれるようにすると、そういう余計な空白は消えてくれる。

yasnippetのインストール

典型的なコードスニペットを入力しやすくするモジュール。これを入れると、例えば「ifm」と入力して[tab]を押すと、

if __name__=='__main__':
    main()

に自動的に変換してくれる。

インストールは、emacsから「M-x package-install」と入力して、次にパッケージを聞かれたら「yasnippet」と入力。

.emacs.d/init.elに以下を加える。

(yas-global-mode t)

py-autopep8のインストール

Pythonのプログラムを保存するときに自動的にPEP8準拠の形式に変換してくれるように設定する。

まずはシェルからautopep8をインストールする。

pip install --user autopep8

次にemacsを起動して「M-x package-install」と入力して、次にパッケージを聞かれたら「py-autopep8」と入力。

.emacs.d/init.elに以下を加える。

(add-hook 'python-mode-hook 'py-autopep8-enable-on-save)

auto-completeのインストール

Pythonに限らない入力自動補完のモジュール。

次にemacsを起動して「M-x package-install」と入力して、次にパッケージを聞かれたら「auto-complete」と入力。

.emacs.d/init.elに以下を加える。

(require 'auto-complete-config)
(ac-config-default)

elpyのインストール

Pythonのオートコンプリートをしてくれるモジュールのインストール。くれる。

まずはPythonにropeとjediを入れる。

pip install --user rope jedi

次にemacsを起動して「M-x package-install」と入力して、次にパッケージを聞かれたら「elpy」と入力。

flymake関連のインストール

入力中にpyflakesを実行して、エラーを指摘するようにする。

まずはpyflakesをインストールする

pip install --user pyflakes

次に、emacsに「M-x list-packages」を使って、flymake, flymake-cursor, flymake-easy, flymake-python-pyflakesの4つのパッケージをインストールする。

最後に.emacs.d/init.elに以下を追加

(require 'tramp-cmds)
(when (load "flymake" t)
  (defun flymake-pyflakes-init ()
     ; Make sure it's not a remote buffer or flymake would not work
     (when (not (subsetp (list (current-buffer)) (tramp-list-remote-buffers)))
      (let* ((temp-file (flymake-init-create-temp-buffer-copy
                         'flymake-create-temp-inplace))
             (local-file (file-relative-name
                          temp-file
                          (file-name-directory buffer-file-name))))
        (list "pyflakes" (list local-file)))))
  (add-to-list 'flymake-allowed-file-name-masks
               '("\\.py\\'" flymake-pyflakes-init)))

(add-hook 'python-mode-hook
          (lambda ()
            (flymake-mode t)))

動作デモ(動画)

上記モジュールをすべて入れた状態で、実際に動かしたときの動作の動画を公開する。(みんな大好き)フィボナッチ数列を出力するプログラムを書いてから実行するまでの様子を示している。

参考サイト

EmacsWiki: Programming Python in Emacs

更新履歴

2015/7/29 py-autopep8についての記述が古くなったので、一部変更した(こちら参照)

scipy.sparseの内部データ構造


前回の記事で予告した記事が遅くなってしまった。申し訳ない。

ここでは、Pythonの疎行列ライブラリscipy.sparseについて、データの内部構造を説明する。内部構造を知らずにブラックボックスとして使いたい場合は前回の記事を参照のこと。

はじめに(あるいは前回記事のまとめのまとめ)

scipy.sparseで疎行列を表すクラスは主に以下のものがある(他にもあるが、ここではとりあげない。これだけ知ってれば大体のことはできる。)

  • lil_matrix
  • csr_matrix
  • csc_matrix

通常は、lil_matrixを用意してそこに値を詰めてから、csr_matrixまたはcsc_matrixに変換してから計算を行う。csr_matrixとcsc_matrixのどちらを使うかは以下のことから判断する。

  • 演算はcsr同士またはcsc同士が高速
  • csrは行を取り出すのが高速
  • cscは列を取り出すのが高速
  • csrは転置をとるとcscになり、cscは転置をとるとcsrになる

内部構造

以下csr_matrixとcsc_matrixのデータ内部構造について説明する。このCSRとCSCという形式は、Pythonに限定された話ではなく、数値計算で疎行列を扱うときに一般的に使われるものである。

用語

以下、行と列の名前は、行0, 行1, 行2, …, 列0, 列1, 列2, …という風に0からナンバリングして呼ぶことにする。「i行目」「j列目」というときも、0を起点として考える。

csr_matrix

内部データは、3つのメンバ変数data, indices, indptrで構成される。

  • dataには方向に非ゼロ要素の値が格納される。(dataのサイズは非ゼロ要素の数)
  • indices[j]はdata[j]が何目の要素かを示す。(indicesのサイズはdataのサイズに等しい)
  • jがindptr[i]以上indptr[i+1]未満の範囲にあるときに、data[j]およびindices[j]はi目の要素である。つまりindptrはdataおよびindicesにおいて、どこからどこまでが同じかを示す。(indptrのサイズは行列の数+1)

実際の例を以下に示す。
csr_matrix

上図の行列で、行0、行1、行2にはそれぞれ、2個、1個、2個の非ゼロ要素がある。したがって、dataおよびindicesは両方とも2,1,2のサイズに区切り、それぞれの行の情報を格納する。行0の非ゼロ要素の値は1., 2.だからdata[0], data[1]には1., 2.が格納され、行0の非ゼロ要素の列インデックスは1,2だから、indices[0], indices[1]には1,2が格納される。行0の情報はdata[0]とindices[0]から始まるので、indptr[0]には0が入り、行1の情報はdata[2]とindices[2]から始まるので、indptr[1]には2が入る。以下同様に格納される。

実際の実行セッションは以下のようになる。

>>> from scipy.sparse import lil_matrix, csr_matrix
>>> a=lil_matrix((3,3))
>>> a[0,1]=1.; a[0,2]=2.; a[1,2]=3.; a[2,0]=4.; a[2,1]=5.
>>> b=a.tocsr()
>>> b.todense()
matrix([[ 0.,  1.,  2.],
        [ 0.,  0.,  3.],
        [ 4.,  5.,  0.]])
>>> b.indices
array([1, 2, 2, 0, 1], dtype=int32)
>>> b.data
array([ 1.,  2.,  3.,  4.,  5.])
>>> b.indptr
array([0, 2, 3, 5], dtype=int32)

csc_matrix

csr_matrix同様に内部データは、3つのメンバ変数data, indices, indptrで構成される。csrの説明で、「行」と「列」を入れ替えるとcsc_matrixの内部構造の説明になる。つまり以下のようなことになる。

  • dataには方向に非ゼロ要素の値が格納される。(dataのサイズは非ゼロ要素の数)
  • indices[j]はdata[j]が何目の要素かを示す。(indicesのサイズはdataのサイズに等しい)
  • jがindptr[i]以上indptr[i+1]未満の範囲にあるときに、data[j]およびindices[j]はiの要素である。つまりindptrはdataおよびindicesにおいて、どこからどこまでが同じかを示す。(indptrのサイズは行列の数+1)

具体例は図だけを示す。
csc_matrix

内部構造を知っていると何がうれしいか

内部構造を理解することで、高速な操作ができるケースがある。非ゼロ要素に同一関数を作用させるケースについては以前のブログで書いたが、ここではコンストラクタの高速化の話をする。

例えば新しく、行列をつくるとき、(公式ドキュメントのユースケースにかかれているように)lil_matrixに詰めてからcsr_matrix(またはcsc_matrix)に変換するよりも、csr_matrix(またはcsc_matrix)の内部構造を直接用意してコンストラクトしたほうが処理速度が速い。

実際に、次のような行列について、初期化にどのくらいかかるかのベンチマークを取ってみる。
\[
\begin{pmatrix}
2&1& & &\cr
&2&1& &\cr
& &\ddots&\ddots&\cr
& & &\ddots&1 \cr
& & & &2
\end{pmatrix}
\]

ベンチマークに使ったコードは以下の通り。ここではサイズ10000×10000の行列で実験している。

from scipy.sparse import lil_matrix, csr_matrix
import numpy as np
from timeit import timeit

def set_lil(n):
    a=lil_matrix((n,n))
    for i in xrange(n):
        a[i,i]=2.
        if i+1<n:
            a[i,i+1]=1.
    return a

def set_csr(n):
    data=np.empty(2*n-1)
    indices=np.empty(2*n-1,dtype=np.int32)
    indptr=np.empty(n+1,dtype=np.int32)
    # to be fair, for-sentence is intentionally used
    # (using indexing technique is faster)
    for i in xrange(n):
        indices[2*i]=i
        data[2*i]=2.
        if i<n-1:
            indices[2*i+1]=i+1
            data[2*i+1]=1.
        indptr[i]=2*i
    indptr[n]=2*n-1
    a=csr_matrix((data,indices,indptr),shape=(n,n))
    return a

print "lil:",timeit("set_lil(10000)",
                    number=10,setup="from __main__ import set_lil")
print "csr:",timeit("set_csr(10000)",
                    number=10,setup="from __main__ import set_csr")

そして、ある環境における実行結果が以下の通り。

lil: 11.6730761528
csr: 0.0562081336975

csr_matrixに直接詰めたほうが圧倒的に速い。とはいえ、csr_matrixに直接詰めるコードの方が複雑になるので、開発効率の面からlil_matrixを利用する選択肢も捨てがたいかと思う。csr_matrix(またはcsc_matrix)に直接詰めるのが有効なのは以下のような場合かと考えられる。

  • 疎行列を何度も用意しなければいけない場合(一度用意するだけならあまりパフォーマンスの違いは気にならないかもしれない)
  • データがあらかじめある程度ソートされた形で用意されている場合(オンラインアルゴリズムのような、次にどの非ゼロ要素がくるかわからないような状況だと、lil_matrixを使ったほうが手っ取り早いし、そこでわざわざソートするのならパフォーマンスの優位性は薄れる)

まとめ

  • scipy.sparseを使うとき、通常はlil_matrixを用意して値を入れて、csr_matrixかcsc_matrixに変換してから計算する。そのように使っている限りは内部構造を知る必要はなく、ブラックボックスとして使える。
  • しかし場合によっては、csr_matrixやcsc_matrixに直接値を詰めたり、中身をいじったりするほうが、計算速度が速くなる。そのためには内部データ構造を理解する必要がある。

追記:
2015/4/1 csc_matrixの図が間違っていたので修正しました。指摘していただいた@KeWata09さんに感謝します。

Pythonの疎行列ライブラリscipy.sparseの基本

今回は珍しく(?)初心者フレンドリーに書こうと思う。先日のPyCon JPでの発表で、Pythonで使える疎行列のライブラリ「scipy.sparse」について反響があったので、改めてここで使い方の基本をまとめてみることとする。そもそも、便利で使い方もそれほど難しくないのに、日本語の情報がほとんどないようなので敬遠されていることもあるかもしれない。ここでは、使い方の基本を解説する。

基本: 疎行列とは?

要素のほとんどが0であるような行列を疎行列と呼ぶ。そうではない普通の行列を、区別したいときには密行列と呼ぶ。疎行列では、非ゼロ要素だけを覚えておけば良いので、メモリ、計算時間ともに大幅な節約になる。特に非ゼロ要素の割合が小さい場合は、その節約は大きくなる。

特に機械学習の分野だと、計算過程で大規模な疎行列が必要になることがよくある。

使い方

scipy.sparseには、疎行列を表すクラスがいくつかあるが、一般の行列について演算をしたいのならlil_matrix, csc_matrix, csr_matrixだけを覚えれば良い。

典型的な使い方としては
1) lil_matrixを用意する
2) 値をセットする
3) csr_matrix(またはcsc_matrix)に変換する(どちらを使うべきかの判断は後述)
4) 演算をする
という流れになる。

つまりまとめると、

  • 値をセットするときはlil_matrixを使う
  • 計算するときはcsr_matrixまたはcsc_matrixを使う

lil_matrixでも演算はできるが、計算効率を考えるとそれは避けた方がいい。

コード例は以下のようになる。

from scipy.sparse import lil_matrix, csr_matrix

# 疎行列aを用意する
a=lil_matrix((3,3))

# 非ゼロ要素を設定する
a[0,0]=1.; a[0,2]=2.

# lil_matrixをcsr_matrixに変換する
a=a.tocsr()

# 疎行列bを用意する
b=lil_matrix((3,3))

# 非ゼロ要素を設定する
b[1,1]=3.; b[2,0]=4.; b[2,2]=5.

# lil_matrixをcsr_matrixに変換する
b=b.tocsr()

# aとbの積を計算する
c=a.dot(b)

# aとbの和を計算する
d=a+b

このコードでは、
\[
A=
\begin{pmatrix}
1&0&2\cr
0&0&0\cr
0&0&0
\end{pmatrix},\quad
B=
\begin{pmatrix}
0&0&0\cr
0&3&0\cr
4&0&5
\end{pmatrix}
\]
という行列について
\[
C=AB,\quad D=A+B
\]

という計算をしている。もちろん、これはユースケースを示しているだけで、実際に疎行列が威力を発揮するのは大規模行列のときである。3×3程度の大きさだと疎行列を使うメリットはない。

lil_matrixのコンストラクタでは、引数として(m,n)のように行列のサイズだけをしていすれば、まず全要素が0の行列ができる。そうしてできた行列をaとすると、a[i,j]=xの形で非ゼロ要素だけセットする。その後、a.tocsr()またはa.tocsc()で、csr_matrixまたはcsc_matrixに変換してから、和・積などの計算をする。

CSRかCSCか?

では、計算するときにはcsr_matrixにすべきか、csc_matrixにすべきか、という問題が残るが、それは以下のルールを考慮した上、適宜よい方を選択すればよい。

  • csr_matrixは〜行目を取り出す操作、csc_matrixは〜列目を取り出す操作が高速である
  • 同じ型同士の和・積は高速である。つまりcsr_matrix同士の和・積、csc_matrix同士の和・積はともに高速である
  • csr_matrix, csc_matrixは転置行列を取ると型が移り合う。つまりcsr_matrixの転置をとるとcsc_matrixになり、csr_matrixの転置をとるとcsc_matrixになる。

以上の説明で、「高速になる」というのは、大規模データの場合はすべて「それ以外はやってはいけない」と読み替えても間違いではない。例えば、csr_matrixとcsc_matrixの和・積の計算はやってはいけない。

ではなぜこうなっているのかというと、内部構造の理解が必要になるが、それはPyCon JPでもしゃべったが次回の記事で書こうと思う。

更新履歴:
2016/10/24 最初の例で行列Aの説明が間違っているとの指摘を受けましたので、訂正しました。

PyCon JP 2014で発表してきた

先日PyCon JP(Python開発者のお祭り)で、登壇発表してきました。
そのときの資料がこちら

内容としては、NumPyを使って高速に計算するコツと、SciPyの中の特に疎行列の扱いについてです。

そしてビデオがこちら。
https://www.youtube.com/watch?v=VJwjzY6jdBY#t=293

Togetterでの反応
http://togetter.com/li/719865?page=23
を見ると、やはり難しかったという人もいるようですが、こういう内容の発表だとやはり全員にその場で理解してもらうというのは難しいと思うので、今後ブログ等で補完的な内容を発信していこうかと思ってます。

人間が読めない数字

先日の記事「コンピュータが読めない数字」で、「コンピュータさん悪くない」って言ったが、どう悪くないかイマイチ伝わってない気がしたので、MNISTデータで個人的に納得できないものワーストを書き出してみようと思う。

worst
これすべて数字なのだが、読めるだろうか。
正解は下をスクロールして参照。

答:
左から
9 7 3 2 9 5 4
(MNISTデータに付随している正解データによる)