このブログの更新は Twitterアカウント @m_hiyama で通知されます。
Follow @m_hiyama

メールでのご連絡は hiyama{at}chimaira{dot}org まで。

はじめてのメールはスパムと判定されることがあります。最初は、信頼されているドメインから差し障りのない文面を送っていただけると、スパムと判定されにくいと思います。

［参照用 記事］

Scilabで行列計算をしてみる

雑記／備忘

今までR言語で行列計算してたんですが、Scilabに替えてみました。行列計算に限れば、Scilabのほうがずっと使い勝手がいいです。

内容：

R言語には不満が募ってきた

数値的行列計算が必要なとき僕は、R言語を使っていました。Rコンソールを立ち上げれば、その場ですぐ計算できるので便利ではあるのですが、いくつか不満もあります。

  1. 行列を直接表現するリテラル表記がなくて、matrix(c(1, 3, 2, 4), 2, 2) のような関数呼び出しが必要でうっとうしい。
  2. 行列の掛け算記号が %*% 。基本演算子が3打鍵必要なのがイヤ。
  3. 行列のブロック操作がサクサクできないのでフラストレーションを感じる。

ここでRの悪口を言っても生産的ではないので詳細は書きませんが、3番目の不満が意味不明でしょうから説明します。「行列のブロック操作」と書きましたが、正式になんと言うのか知りません。例示するなら、

A = \begin{bmatrix}1 & 2 \\ 3 & 4\end{bmatrix} \\B = \begin{bmatrix}5 & 6 \\ 7 & 8\end{bmatrix}

という行列A, Bがあったとき、例えばAとBを対角位置に並べてCを作る操作とかです。

C = \begin{bmatrix}A & O \\ O & B\end{bmatrix} = \begin{bmatrix}1 & 2 & 0 & 0 \\3 & 4 & 0 & 0 \\0 & 0 & 5 & 6 \\0 & 0 & 7 & 8 \end{bmatrix}

もっと簡単な例だと：

X = \begin{bmatrix}1 \\ 2 \end{bmatrix} \\Y = \begin{bmatrix}3 \\ 4\end{bmatrix} \\Z = \begin{bmatrix}X & Y \end{bmatrix} = \begin{bmatrix} 1 && 3 \\ 2 && 4 \end{bmatrix}

Scilab

ちょっと調べてみて、Scilab（サイラボ）がヨサゲな印象でした。Scilabのホームページは：

https://www.scilab.org/download/latestから各プラットフォームごとのインストール用アーカイブ（バイナリ配布物）をダウンロードできます。現時点（2017年春）での最新バージョンは Scilab 6.0.0 です。インストールは簡単です、特に説明することはないです。

https://www.scilab.org/resources/documentation/tutorials から、日本語PDF文書Scilab_beginners_jp.pdf（はじめてのScilab -- 33ページ）、Xcos_beginners_jp.pdf（はじめてのXcos -- 15ページ）をダウンロードできます。もう少し詳しい解説はintroscilab.pdf（INTRODUCTION TO SCILAB -- 87ページ）。文書は若干古めで、最新バージョンとはズレがあります。が、深刻な問題にはならないでしょう。

ソフトウェアに付属のヘルプも日本語化されていて（見出しは英語のまま）、かなり充実しています。ブラウザでマニュアルを見たいなら、インターネット上のページがあります。

ソフトウェア付属のヘルプとWebオンラインマニュアルは、最新バージョンの情報だと期待してよさそうです。


どうでもいいことですが、Scilabのlabを「ラボ」と読むのは「ラボラトリー（laboratory）」の略だという前提があるからで、labだけなら「ラブ」に近いと思うのだけど…。僕は「サイラブ」と読んで（呼んで）しまいますね。

[追記]次のサイトに、日本語を含む多言語による情報があります。

機械翻訳を使っているのか？ 日本語の品質はよくありません。何言ってるのかワケワカンナイことがあるので注意してください。
[/追記]

Scilabコンソール

Scilabをインストールしたら起動します。初期状態では、メインウィンドウ内にいくつかのサブウィンドウが詰め込まれた最近のIDEっぽい外観です。ゴチャゴチャしていてい僕は好きじゃないのですが、カスタマイズとかは後回しにしてスーパー電卓（a super desktop calculator）*1としてScilabを使ってみます。

「Scilab 6.0.0 コンソール」（日本語UI）とタイトルされたペイン（サブウィンドウ）がScilabコンソールです。このコンソールが対話的インタプリタのコマンドラインインターフェイスです。今回はScilabコンソールしか使いません。

「-->」がインタプリタのプロンプトマークで、この後にキーボード入力して、Enter/Returnキーを押せばScilabインタプリタが応答を返します。

--> disp("Hello, world")

 Hello, world


--> 1 + 2
 ans  =

   3.


--> x = 1 + 2
 x  =

   3.


--> y = x*2 + 1
 y  =

   7.


--> sqrt(y) // yの平方根
 ans  =

   2.6457513

1行目のdisp()はプリント用関数です。1 + 2 のような計算式を入力すれば結果が表示されますが、その結果はansという変数に代入されます。ただし、明示的な代入文を書けば（x = 1 + 2）、ansへの自動代入はされません。計算は倍精度浮動小数点数（64ビット実数）で行われ、値がちょうど整数でも小数点は表示されます*2。関数呼び出しは、たいていのプログラミング言語と同じ func(arg) の形です。スラッシュ2個から行末まではコメントです。

行列計算をはじめる

Scilabでは、行列を非常に分かりやすい形式で入力できます。1行分のデータはカンマか空白で区切り、次の行に移るにはセミコロンを使います。

--> [1, 2; 3, 4]
 ans  =

   1.   2.
   3.   4.


--> [1 2; 3 4] // カンマの代わりに空白でもよい
 ans  =

   1.   2.
   3.   4.


--> a = [1 2; 3 4] // 行列を変数に代入
 a  =

   1.   2.
   3.   4.

行列の掛け算はアスタリスク1文字です。加減はもちろんプラスとマイナス。

--> b = [5 6; 7 8]
 b  =

   5.   6.
   7.   8.


--> a*b
 ans  =

   19.   22.
   43.   50.


--> a + b
 ans  =

   6.    8.
   10.   12.


--> a - b
 ans  =

  -4.  -4.
  -4.  -4.

スカラー倍もアスタリスクです。

--> 2*a
 ans  =

   2.   4.
   6.   8.


--> a*-1
 ans  =

  -1.  -2.
  -3.  -4.

Scilabにベクトル型データはありません*3が、1行n列-行列かn行1列-行列でベクトルを表現できます。どっちにするかは「決め」（約束）の問題です。例えば、ベクトルをn行1列-行列＝列ベクトルで表現するとして、行列とベクトルの積は次のようです。

--> a*[2; -1]
 ans  =

   0.
   2.


--> b*[2; -1]
 ans  =

   4.
   6.

行列の転置（縦と横の入れ替え）は、'（シングルクォート）を後置します*4。

--> a'
 ans  =

   1.   3.
   2.   4.


--> [2; -1]'
 ans  =

   2.  -1.

ベクトルの内積は転置して掛ければいいので：

--> x = [2; -1]; y = [5; 3]; // 2つの代入文を行って、結果の表示は行わない

--> x'*y // ベクトルxとyの内積
 ans  =

   7.

行列のブロック操作と寸法

特別な行列を表現する関数が用意されています。対角成分が1で他は0である行列は eye(行数, 列数) で、すべての成分が0である行列は zeros(行数, 列数) で作れます。

--> eye(2, 2)
 ans  =

   1.   0.
   0.   1.


--> eye(3, 2)
 ans  =

   1.   0.
   0.   1.
   0.   0.


--> zeros(3, 2)
 ans  =

   0.   0.
   0.   0.
   0.   0.

冒頭に挙げたブロック操作をやってみます。

--> a
 a  =

   1.   2.
   3.   4.


--> b
 b  =

   5.   6.
   7.   8.


--> o = zeros(2, 2)
 o  =

   0.   0.
   0.   0.


--> c = [a o; o b]
 c  =

   1.   2.   0.   0.
   3.   4.   0.   0.
   0.   0.   5.   6.
   0.   0.   7.   8.


--> x = [1; 2]; y = [3; 4];

--> z = [x y]
 z  =

   1.   3.
   2.   4.

簡単でしょ。

行列の行数と列数が欲しいときはsize()関数、成分の個数が欲しいときはlength()関数を使います。

--> size(c)
 ans  =

   4.   4.


--> length(c)
 ans  =

   16.


--> size([1; 2])
 ans  =

   2.   1.


--> size([1; 2]')
 ans  =

   1.   2.


--> length([1; 2])
 ans  =

   2.

成分と部分行列の取り出し

多くのプログラミング言語では、配列、リスト、タプル、ベクトル、行列などの成分（要素、項目）を取り出すにはブラケットを使って、x[2], a[1, 2] のように書きます。Scilibでは、丸括弧で x(2), a(1, 2) のように書きます。

--> x
 x  =

   1.
   2.


--> x(1)
 ans  =

   1.


--> a
 a  =

   1.   2.
   3.   4.


--> a(1, 2)
 ans  =

   2.

関数呼び出しと区別が付かなくて困らないか？ 困りません。それどころか、区別が付かないのでかえって便利です、後で例を挙げます。JVM上で動く関数型言語Scalaでは、やはりインデックス・アクセスと関数呼び出しが同じ形式ですね。

次に、行列の単一成分に限らず、いくつかの成分をまとめて取り出してみます。そのために範囲記法を説明しておくと、n:m で「nからmまでの整数範囲」*5を意味します。

--> 0:3
 ans  =

   0.   1.   2.   3.


--> 1:10
 ans  =

   1.   2.   3.   4.   5.   6.   7.   8.   9.   10.

範囲の実体は横ベクトル（1行n列の行列）です。この範囲記法を使うと、行列の一部を素早く抜き出せます。

--> c
 c  =

   1.   2.   0.   0.
   3.   4.   0.   0.
   0.   0.   5.   6.
   0.   0.   7.   8.


--> c(1:4, 2:3)
 ans  =

   2.   0.
   4.   0.
   0.   5.
   0.   7.


--> c([1 2 3 4], [2 3])
 ans  =

   2.   0.
   4.   0.
   0.   5.
   0.   7.


--> c([2 3], [2 3])
 ans  =

   4.   0.
   0.   5.

範囲は横ベクトルなので、直接横ベクトルを書いても同じです -- 長い範囲だと直接書くのはシンドイですがね。範囲を示すコロンの特別な用法として、次のような書き方ができます。

--> c(:, [2, 3])
 ans  =

   2.   0.
   4.   0.
   0.   5.
   0.   7.


--> c(:, 2)
 ans  =

   2.
   4.
   0.
   0.


--> c(1, :)
 ans  =

   1.   2.   0.   0.

コロンだけだと、その行列が持つ行範囲、列範囲になるのです。行列のサイズを意識せずに「行全部／列全部」を指定できるので便利です。

抜き出す部分は、範囲では記述できないパランパランな部分でもいいです。

--> c
 c  =

   1.   2.   0.   0.
   3.   4.   0.   0.
   0.   0.   5.   6.
   0.   0.   7.   8.


--> c(:, [1 4])
 ans  =

   1.   0.
   3.   0.
   0.   6.
   0.   8.

インデックス・アクセスと関数呼び出し

前節で、「インデックス・アクセスと関数呼び出しを区別しないほうが便利だ」と書きました。実例を挙げましょう。

整数引数の関数fと、整数n, m（n≦m）に関して次の総和を考えます。

\sum_{i = n}^{m} f(i)

これを、関数f, 整数n, mを引数とする関数（高階関数）と考えて、

total(f,\, n,\, m) \: := \sum_{i = n}^{m} f(i)

としましょう。これをScilabの関数で書くと次のようになります。プログラムとして見やすいように、f→fun、n→from、m→to とリネームしています。

function ret = total(fun, from, to)
    ret = 0
    for i = from:to do
        ret = ret + fun(i)
    end
endfunction

Scilabの関数定義構文の説明はしませんが、何となくは分かるでしょう。

さて、奇数の総和は平方数になるという事実があります。

\sum_{i = 1}^{n} 2n - 1 \: = \: n^2

この事実を関数total()を使って確認してみます。

--> function ret = odd(n)
  >     ret = 2*n - 1
  > endfunction

--> odd(1)
 ans  =

   1.


--> total(odd, 1, 4)
 ans  =

   16.


--> total(odd, 1, 5)
 ans  =

   25.

関数odd()を定義して、関数（高階関数）total()に渡したのですが、奇数の列は別に関数でなくてもかまいません。

--> 1:2:9
 ans  =

   1.   3.   5.   7.   9.


--> total(1:2:9, 1, 4)
 ans  =

   16.


--> total(1:2:9, 1, 5)
 ans  =

   25.

ここで、1:2:9は範囲（実体は行ベクトル）ですが、1から2ずつ増加して9に至る整数並びです。

関数も行ベクトル（あるいは列ベクトル）も同じ形式でアクセスできるので、total()の第1引数は関数でもベクトル（特殊な行列）でもいいのです。計算手順とデータを統一的に扱えて便利です。

逆行列と行列式

行列Aの逆行列はA-1と書きます。この書き方にならって、Sailabでも a^-1 という記法で逆行列を表現できます。また、逆数と同じ 1/a でも逆行列です。

--> a
 a  =

   1.   2.
   3.   4.


--> a^-1
 ans  =

  -2.    1.
   1.5  -0.5


--> 1/a
 ans  =

  -2.    1.
   1.5  -0.5


--> a*1/a
 ans  =

   1.   0.
   0.   1.

ウーン、まーいいんですけど、僕は inv(a) という書き方が落ち着きます。

--> inv(a)
 ans  =

  -2.    1.
   1.5  -0.5


--> a*inv(a)
 ans  =

   1.          0.
   8.882D-16   1.


--> inv(a)*a
 ans  =

   1.          0.
   2.220D-16   1.

アレッ？ a*inv(a)が単位行列になってませんね。2.220D-16はものすごく小さい数なので、浮動小数点数の誤差でしょうが、1/aだと正確に0だったのは不可解。ちょっと事情が分かりません。

逆行列の存在・非存在は行列式を見れば分かりますが、aの行列式はdet(a)です。

--> a
 a  =

   1.   2.
   3.   4.


--> det(a)
 ans  =

  -2.


--> inv(a)
 ans  =

  -2.    1.
   1.5  -0.5


--> det(inv(a))
 ans  =

  -0.5


--> det(a*inv(a))
 ans  =

   1.

おわりに

というわけで、Scilabだと行列計算が気持ちよく出来るよ というお話でした。

*1:"INTRODUCTION TO SCILAB"の22ページに、Scilabは"a super desktop calculator"であり、またそれ以上のものである、と書いてあります。

*2:デフォルトが実数（浮動小数点数）ですが、n = int32(10) のような型変換を使うと、適当なビット長（この場合は32ビット）の整数データを作れます。整数の10は、10. ではなくて 10 と表示されます。

*3:実は、厳密に言えばスカラー型もありません。1行1列の行列がスカラーとして扱われています。

*4:'（シングルクォート）は実際はエルミート転置（Hermitian transpose, conjugate transpose）です。単なる転置は.'（ドットとクォート）という演算子ですが、実数成分の場合は単なる転置もエルミート転置も変わりはありません。

*5:x:y という書き方の意味は、xを初項とする等差数列のyを越えない部分です。例えば、1.2:5 の値は [1.2, 2.2, 3.2, 4.2] となります。今の文脈では、初項と上限に整数値（データ型は実数）を取っているので整数範囲を表すことになっています。