いものやま。

雑多な知識の寄せ集め

Swiftでの行列計算について調べてみた。（その1）

技術 Swift

ニューラルネットワークの実装をするときに必要となるのが、行列計算。
もちろん、Arrayを使ってベタに実装してもいいんだけど、それだと速度が気になるところ。
調べてみると、Accelerateフレームワークを利用するといいようだったので、少し調べてみた。

なお、参考にしたのは以下のサイト：

Accelerateフレームワーク

行列計算のAPIとしてデファクトスタンダードとなっているのが、BLASというAPIらしい。
そして、その上位のパッケージとして、LAPACKがあるらしい。

Accelerateフレームワークにはこれらの実装が含まれていて、高速な行列計算が出来るようになっている。

ただ、これらのAPIは命名規則にかなりクセがあって、分かりにくい。。。

そこで、iOS 8.0（もしくはMac OS X 10.10）以降では、分かりやすいインタフェースが用意されている。
ここでは、その分かりやすいインタフェースを使っていくことにする。

なお、Accelerateフレームワークを使えるようにするには、リンクするライブラリにAccelerateフレームワークを追加し、ソースコードで以下のようにインポートを行う：

import Accelerate

la_object_t

行列（やベクトル、および、後述のSplat）を表す型として、la_object_tという型が用意されている。
ちなみに、接頭辞のlaはLinearAlgebraの略。

メモリをARCで管理している場合、このオブジェクトのメモリ管理は他のオブジェクトと同じようにすればいい。
（なお、ARCで管理していなかったり、C言語で使う場合には、la_retain()およびla_release()でメモリ管理を行う必要がある）

今見ていこうとしているインタフェースは、C言語でも使えるようにC言語の関数の形で書かれているので、C言語でオブジェクト指向プログラミングを行うときと同様に、このオブジェクトを関数の第1引数で取るようになっているものが多い。

行列の生成

まずは行列の生成から。

配列から生成する

行列を生成する基本的な方法として、Floatの配列、もしくはDoubleの配列から、行列を生成する関数が用意されている：

// Float/Doubleのバッファ（配列）から行列を生成する。
// （バッファの内容はコピーされる）
la_object_t
la_matrix_from_float_buffer(
    const float *buffer,           // Floatの配列
    la_count_t matrix_rows,        // 行数
    la_count_t matrix_cols,        // 列数
    la_count_t matrix_row_stride,  // バッファの各行（2行目以降）のオフセット
    la_hint_t matrix_hint,         // 行列の特徴
    la_attribute_t attributes);    // 属性
la_object_t
la_matrix_from_double_buffer(
    const double *buffer,          // Doubleの配列
    la_count_t matrix_rows,        // 行数
    la_count_t matrix_cols,        // 列数
    la_count_t matrix_row_stride,  // バッファの各行（2行目以降）のオフセット
    la_hint_t matrix_hint,         // 行列の特徴
    la_attribute_t attributes);    // 属性

なお、la_count_tはUIntの別名。

la_hint_tは次のように定義されていて、行列の特徴を指定でき、正しく指定すれば計算の効率を上げることが出来る：

typedef unsigned long la_hint_t;
#define LA_NO_HINT                       (0U)        // ヒントなし
#define LA_SHAPE_DIAGONAL                (1U << 0)   // 対角行列
#define LA_SHAPE_LOWER_TRIANGULAR        (1U << 1)   // 下三角行列
#define LA_SHAPE_UPPER_TRIANGULAR        (1U << 2)   // 上三角行列
#define LA_FEATURE_SYMMETRIC             (1U << 16)  // 対称行列
#define LA_FEATURE_POSITIVE_DEFINITE     (1U << 17)  // 正定値行列
#define LA_FEATURE_DIAGONALLY_DOMINANT   (1U << 18)  // 対角優位行列

la_attribute_tは次のように定義されていて、普通はLA_DEFAULT_ATTRIBUTESを指定しておけばいい：

typedef unsigned long la_attribute_t;
#define LA_DEFAULT_ATTRIBUTES           (0)
#define LA_ATTRIBUTE_ENABLE_LOGGING     (1U << 0)

ちょっと分かりにくいのが、matrix_row_strideという引数。
これは、1次元で並んでいる配列を2次元に並べていくときに、2行目、3行目、・・・のデータの開始位置を決定するための値。

例えば、この値が5なら、行列の1行目のデータは配列のインデックスの0から、行列の2行目のデータは配列のインデックスの5から、行列の3行目のデータは配列のインデックスの10から、取ってくることになる。
（ただし、重なってしまうので、列数より小さい値は指定できない）

サンプルを書くと、以下のようになる：

import Foundation
import Accelerate

let originalBuffer: [Double] = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]

let matrix1 = la_matrix_from_double_buffer(originalBuffer, 3, 2, 2,
                                           la_hint_t(LA_NO_HINT),
                                           la_attribute_t(LA_DEFAULT_ATTRIBUTES))
// matrix1は次のようになる：
// 1.0 2.0
// 3.0 4.0
// 5.0 6.0

let matrix2 = la_matrix_from_double_buffer(originalBuffer, 3, 2, 3,
                                           la_hint_t(LA_NO_HINT),
                                           la_attribute_t(LA_DEFAULT_ATTRIBUTES))
// matrix2は次のようになる：
// 1.0 2.0
// 4.0 5.0
// 7.0 8.0

なお、行列の生成に成功したかどうかというのは、la_object_tのオブジェクト自体が保持している。
（関数の戻り口が1つしかないC言語では、常套手段）
それを知るには、次の関数を使う：

// 状態を取得する。
la_status_t
la_status(
    la_object_t object);  // 行列

la_status_tは次のように定義されている：

typedef long la_status_t;
#define LA_SUCCESS                       (0)
#define LA_WARNING_POORLY_CONDITIONED    (1000)
#define LA_INTERNAL_ERROR                (-1000)
#define LA_INVALID_PARAMETER_ERROR       (-1001)
#define LA_DIMENSION_MISMATCH_ERROR      (-1002)
#define LA_PRECISION_MISMATCH_ERROR      (-1003)
#define LA_SINGULAR_ERROR                (-1004)
#define LA_SLICE_OUT_OF_BOUNDS_ERROR     (-1005)

行列の生成に限らず、他の処理についても、処理が正常に終了したかどうかの情報はla_object_t自身が保持するようになっているので、適宜la_status()を呼んで状態をチェックしないといけない。

単位行列を生成する

単位行列を生成するために、次の関数が用意されている：

// 単位行列を生成する。
la_object_t
la_identity_matrix(
    la_count_t matrix_size,        // 行列のサイズ
    la_scalar_type_t scalar_type,  // 精度
    la_attribute_t attributes);    // 属性

la_scalar_type_tは行列の精度を指定するためのもので、次のように定義されている：

typedef unsigned int la_scalar_type_t;
#define LA_SCALAR_TYPE_FLOAT  (0x8000)
#define LA_SCALAR_TYPE_DOUBLE (0x4000)

対角行列を生成する

対角行列を生成するために、次の関数が用意されている：

// 対角行列を生成する。
la_object_t
la_diagonal_matrix_from_vector(
    la_object_t vector,           // 対角成分のベクトル（1行もしくは1列の行列）
    la_index_t matrix_diagonal);  // 右上へいくつズラすか

la_index_tはIntの別名。
matrix_diagonalに正の値を指定すると右上へズレ、負の値を指定すると左下へズレることになる。

サンプルを書くと、以下のようになる：

import Foundation
import Accelerate

let vectorBuffer: [Double] = [1.0, 2.0, 3.0]
let vector = la_matrix_from_double_buffer(vectorBuffer, 3, 1, 1,
                                          la_hint_t(LA_NO_HINT),
                                          la_attribute_t(LA_DEFAULT_ATTRIBUTES))
// vectorは次の列ベクトル：
// 1.0
// 2.0
// 3.0

let matrix1 = la_diagonal_matrix_from_vector(vector, 0)
// matrix1は次のようになる：
// 1.0 0.0 0.0
// 0.0 2.0 0.0
// 0.0 0.0 3.0

let matrix2 = la_diagonal_matrix_from_vector(vector, 1)
// matrix2は次のようになる：
// 0.0 1.0 0.0 0.0
// 0.0 0.0 2.0 0.0
// 0.0 0.0 0.0 3.0
// 0.0 0.0 0.0 0.0

let matrix3 = la_diagonal_matrix_from_vector(vector, -2)
// matrix3は次のようになる：
// 0.0 0.0 0.0 0.0 0.0
// 0.0 0.0 0.0 0.0 0.0
// 1.0 0.0 0.0 0.0 0.0
// 0.0 2.0 0.0 0.0 0.0
// 0.0 0.0 3.0 0.0 0.0

部分的な行列を生成する

すでにある行列から行と列を部分的にピックアップして新しい行列を生成するために、次の関数が用意されている：

// 部分的な行列を生成する。
la_object_t
la_matrix_slice(
    la_object_t matrix,            // 元の行列
    la_index_t matrix_first_row,   // 元の行列からピックアップする最初の行
    la_index_t matrix_first_col,   // 元の行列からピックアップする最初の列
    la_index_t matrix_row_stride,  // 行をピックアップしていくときにズラす幅
    la_index_t matrix_col_stride,  // 列をピックアップしていくときにズラす幅
    la_count_t slice_rows,         // ピックアップする行数
    la_count_t slice_cols);        // ピックアップする列数

例えば、matrix_first_rowに0、matrix_row_strideに2、slice_rowsに3を指定すると、元の行列の0, 2, 4行目がピックアップされることになる。

サンプルを書くと、以下のようになる：

import Foundation
import Accelerate

let originalBuffer: [Double] = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]

let originalMatrix = la_matrix_from_double_buffer(originalBuffer, 3, 3, 3,
                                                  la_hint_t(LA_NO_HINT),
                                                  la_attribute_t(LA_DEFAULT_ATTRIBUTES))
// originalMatrixは次のようになる：
// 1.0 2.0 3.0
// 4.0 5.0 6.0
// 7.0 8.0 9.0

let sliceMatrix1 = la_matrix_slice(originalMatrix,
                                   0, 0,
                                   1, 1,
                                   2, 2)
// sliceMatrix1は次のようになる：
// 1.0 2.0
// 4.0 5.0

let sliceMatrix2 = la_matrix_slice(originalMatrix,
                                   0, 1,
                                   1, 1,
                                   3, 2)
// sliceMatrix2は次のようになる：
// 2.0 3.0
// 5.0 6.0
// 8.0 9.0

// strideには負の値も指定できる。
let sliceMatrix3 = la_matrix_slice(originalMatrix,
                                   0, 2,
                                   1, -1,
                                   3, 3)
// sliceMatrix3は次のようになる：
// 3.0 2.0 1.0
// 6.0 5.0 4.0
// 9.0 8.0 7.0

行列の行（もしくは列）を取得する

行列の行（もしくは列）を取得して、1行の行列（もしくは1列の行列）を返す関数が用意されている：

// 行列から行ベクトル/列ベクトルを取得する。
la_object_t
la_vector_from_matrix_row(
    la_object_t matrix,      // 行列
    la_count_t matrix_row);  // 取得する行
la_object_t
la_vector_from_matrix_col(
    la_object_t matrix,      // 行列
    la_count_t matrix_col);  // 取得する列

行列の対角成分を取得する

行列の対角成分を取得して、1列の行列を返す関数が用意されている：

// 行列から対角成分を取得して列ベクトルとして返す。
la_object_t
la_vector_from_matrix_diagonal(
    la_object_t matrix,           // 行列
    la_index_t matrix_diagonal);  // 右上へいくつズラすか

matrix_diagonalはla_diagonal_matrix_from_vector()と同じ。

行列の行数、列数

行列の行数、列数は、以下の関数で取得できる：

// 行列の行数を返す。
la_count_t 
la_matrix_rows(
    la_object_t matrix);  // 行列

// 行列の列数を返す。
la_count_t
la_matrix_cols(
    la_object_t matrix);  // 行列

行列の内容の取得

行列の内容を見るには、Floatの配列、もしくはDoubleの配列をバッファとして用意して、そこに行列の内容を書き込む必要がある。

配列に行列の内容を書き込む関数は、以下のようになっている：

// バッファとして用意した配列に行列の内容を書き込む。
// 配列は十分なサイズを用意し、inout引数になるので&をつけて渡す。
la_status_t
la_matrix_to_float_buffer(
    float *buffer,                 // バッファとなるFloatの配列
    la_count_t buffer_row_stride,  // バッファの各行（2行目以降）のオフセット
    la_object_t matrix);           // 行列
la_status_t
la_matrix_to_double_buffer(
    double *buffer,                // バッファとなるDoubleの配列
    la_count_t buffer_row_stride,  // バッファの各行（2行目以降）のオフセット
    la_object_t matrix);           // 行列

buffer_row_strideはla_matrix_from_float_buffer()、la_matrix_from_double_buffer()と同様に、行列の2行目、3行目、・・・のデータを書き込むときの開始位置を決定するための値。

サンプルを書くと、以下のようになる：

import Foundation
import Accelerate

let originalBuffer: [Double] = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]

let matrix = la_matrix_from_double_buffer(originalBuffer, 3, 2, 2,
                                          la_hint_t(LA_NO_HINT),
                                          la_attribute_t(LA_DEFAULT_ATTRIBUTES))
// matrixは次のようになる：
// 1.0 2.0
// 3.0 4.0
// 5.0 6.0

var buffer: [Double] = [Double](count: 9, repeatedValue: 0.0)
la_matrix_to_double_buffer(&buffer, 2, matrix)
// bufferは次のようになる：
// 1.0 2.0 3.0 4.0 5.0 6.0 0.0 0.0 0.0

buffer = [Double](count: 9, repeatedValue: 0.0)
la_matrix_to_double_buffer(&buffer, 3, matrix)
// bufferは次のようになる：
// 1.0 2.0 0.0 3.0 4.0 0.0 5.0 6.0 0.0

今日はここまで！