[iOS 8] Accelerate.Frameworkに加わった新たな線形代数ライブラリ

2014.09.18

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

New BLAS

行列やベクトルを処理するためのライブラリ(Basic Linear Algebra Subprograms - BLAS)はこれまでもAccelerate.framework内のvecLib.frameworkで提供されていました。今回のOSアップデートではvecLib.frameworkにLinearAlgebraのフォルダが新しく加わり、行列やベクトルを扱いやすくなりました。

行列、ベクトルを表すla_object_t

この新しいAPI群の中での主役はla_object_tです。これはベクトルや行列一般を表す型です。Objective-Cではオブジェクトとして扱われる為に、ARC環境下ではメモリの煩わしい管理を避けることができます。

Objective-Cでオブジェクトとして扱われている為に、Swiftへの移行もやりやすいです。

本記事ではla_object_tのログを出力するために以下のla_object_t拡張を用意します。

extension la_object_t {
    final public var rows: UInt {
        return UInt(la_matrix_rows(self))
    }

    final public var cols: UInt {
        return UInt(la_matrix_cols(self))
    }

    final public var toArray: [Double] {
        var arrayBuf = [Double](count: Int(rows * cols), repeatedValue: 0.0)
        let status = la_matrix_to_double_buffer(&arrayBuf, cols, self)
        return arrayBuf
    }

    final public var matrix: [[Double]] {
        var matrix = [[Double]]()
        for row in 1...rows {
            let firstCol = Int(cols * (row - 1))
            let lastCol = Int(cols * row - 1)
            let partCols = Array(toArray[firstCol...lastCol])
            matrix.append(partCols)
        }
        return matrix
    }

    final public var matrixDescription: String {
        return matrix.reduce("")
            {(acc, rowVals) in
                acc +
                rowVals.reduce(""){(ac, colVal) in ac + "\(colVal) "} +
                "\n"
            }
    }
}

この拡張を既存のコードに追加することでmatrixDescriptionプロパティ経由でログの為の文字列が取得できます。本記事ではDoubleを行列の要素として扱っていますが、ライブラリではFloatを行列の要素としても扱えます。関数などでdoubleという文字が出てくる場合は適宜float用のAPIも用意されていると読み替えていってください。

遅延評価

新しいライブラリが持つ優れた特徴としては遅延評価があります。

la_object_tの各演算のための関数ではオブジェクト同士で実際に計算が行われるわけではなく、la_object_t同士の関係性が保存されていくだけです。

メモリのアロケーションも各演算では行われず、実際に内部で計算が行われるのは、バッファに値を書き込む次の関数を呼び出した時だけになります。

  • la_matrix_to_float_buffer
  • la_matrix_to_double_buffer
  • la_vector_to_float_buffer
  • la_vector_to_double_buffer

計算結果が知りたい時だけこの関数を呼び出すことにより計算コストの大幅な削減が見込まれます。

APIの機能別にAPI関数とコード例の詳細を合わせて見ていきましょう。

生成のための関数群

まず、la_object_tを生成するための関数についてです。

メモリバッファからの生成

doubleのメモリバッファからの行列(la_object_t)の生成は次の関数でおこないます。

func la_matrix_from_double_buffer(buffer: UnsafePointer<Double>, // メモリバッファ
                                  matrix_rows: la_count_t, // 行数
                                  matrix_cols: la_count_t, // 列数
                                  matrix_row_stride: la_count_t, // 改行ごとにメモリバッファをスライドさせる幅
                                  matrix_hint: la_hint_t, // 計算効率化の為の行列に関する情報指定
                                  attributes: la_attribute_t) // API操作に関する属性
                                  -> la_object_t!

使用例

var matABuf: [Double] =
[1, 2, 3,
4, 5, 6,
7, 8, 9]
let matA = la_matrix_from_double_buffer(matABuf, 3, 3, 3, 
la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))
println(matA.matrixDescription)
// 1.0 2.0 3.0
// 4.0 5.0 6.0
// 7.0 8.0 9.0

la_count_tはUIntのタイプエイリアスです。第5引数と第6引数のla_hint_t型とla_attribute_t型は内部的にはUIntのタイプエイリアスで、それぞれ次のような値で初期化を行います。

la_hint_t

行列が対角行列であるとか上三角行列であるといったヒントに合わせて実際の行列計算の際に最適化を行ってくれ、計算に必要なコストを減らすことが可能になっています。各ヒントを見てみましょう。

  • LA_NO_HINT : 特に特徴を持たない行列を表します。
  • LA_SHAPE_DIAGONAL : 対角行列を表します。
  • LA_SHAPE_LOWER_TRIANGULAR : 下三角行列を表します。
  • LA_SHAPE_UPPER_TRIANGULAR : 上三角行列を表します。
  • LA_FEATURE_SYMMETRIC : 対称行列を表します。
  • LA_FEATURE_POSITIVE_DEFINITE : 正定値行列(二次形式が非負)を表します。
  • LA_FEATURE_DIAGONALLY_DOMINANT : 対角優位行列(対角成分以外の行の絶対値の総和が対角成分の絶対値よりも小さい)を表します。

la_attribute_t

計算に伴う属性の指定をできます。いまのところ、LA_ENABLE_LOGGINGのみがあります。

  • LA_DEFAULT_ATTRIBUTES : デフォルト属性です
  • LA_ATTRIBUTE_ENABLE_LOGGING : デバグ用のログを出力する属性です。エラーが出た時にどこで躓いたのかを探るのに便利です。

単位行列の生成

単位行列は次の関数で取得することができます。

func la_identity_matrix(matrix_size: la_count_t, // 行列のサイズ
                        scalar_type: la_scalar_type_t, // 要素のタイプ、FloatかDouble用の値を指定
                        attributes: la_attribute_t) 
                        -> la_object_t!

使用例

let matE = la_identity_matrix(3, la_scalar_type_t(LA_SCALAR_TYPE_DOUBLE), la_attribute_t(LA_DEFAULT_ATTRIBUTES))
println(matE.matrixDescription)

第二引数のla_scalar_type_tには次の2つの値を指定して初期化します。

  • LA_SCALAR_TYPE_FLOAT : 各要素をfloatの型として初期化します。
  • LA_SCALAR_TYPE_DOUBLE : 各要素をdoubleの型として初期化します。

対角行列と対角行列を上下にずらした行列の生成

対角行列と、対角行列を右上、左下にずらして行列のサイズを大きくしたものは次の関数で取得することができます。

func la_diagonal_matrix_from_vector(vector: la_object_t!, // 列ベクトルから対角成分を決定する
                                    matrix_diagonal: la_index_t) // どのくらいずらすか 
                                    -> la_object_t!

使用例

var vecBBuf: [Double] = [1, 2, 3]
let vecB = la_matrix_from_double_buffer(vecBBuf, 3, 1, 1, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))
let matC = la_diagonal_matrix_from_vector(vecB, 0)
println(matC.matrixDescription)
// 1.0 0.0 0.0 
// 0.0 2.0 0.0 
// 0.0 0.0 3.0 
let matD = la_diagonal_matrix_from_vector(vecB, 1)
println(matD.matrixDescription)
// 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 matF = la_diagonal_matrix_from_vector(vecB, -2)
println(matF.matrixDescription)
// 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_index_tはIntのタイプエイリアスです。

取得のための関数群

メモリバッファへの書き込み

la_object_tの内容をメモリバッファに読み込むために、la_matrix_to_double_buffer関数を使用できます

func la_matrix_to_double_buffer(buffer: UnsafeMutablePointer<Double>, // バッファポインタ
                                buffer_row_stride: la_count_t, // 改行によるバッファのずらし幅
                                matrix: la_object_t!) // バッファに出力するオブジェクト
                                -> la_status_t // 読み込みがうまくいったかどうかを表すステータス

メモリバッファへの書き込みについては先ほどのログを取るための拡張からプロパティを例にして見てみましょう。

使用例

final public var toArray: [Double] {
    var arrayBuf = [Double](count: Int(rows * cols), repeatedValue: 0.0)
    let status = la_matrix_to_double_buffer(&arrayBuf, cols, self)
    return arrayBuf
}

ベクトルに関しても同じような機能をもったla_vector_to_float_buffer, la_vector_to_double_buffer関数が用意されています。

関数の戻り値にla_status_tという型がでてきています。これはInt型のタイプエイリアスです。

どのようなものがあるかをまとめました。

la_status_t

  • LA_SUCCESS : 計算が成功したことを表します。
  • LA_WARNING_POORLY_CONDITIONED : 多くの計算に数値的に条件の不備があることを示しています。計算結果は正確ではありません。
  • LA_INTERNAL_ERROR : LinearAlgebraライブラリ内部のエラーを表しています。例えばメモリの確保が失敗したり、メモリが書き換わったなどの現象です。
  • LA_INVALID_PARAMETER_ERROR : オブジェクトか、計算に用いたオブジェクトを作成する際の関数に渡した引数が不正であったことをしめします。
  • LA_DIMENSION_MISMATCH_ERROR : オブジェクトか、計算に用いたオブジェクトを作成する際の関数に渡された引数の次元に不整合が発生したことを示します。例えば異なる次元のベクトルを足した時などにこのエラーが生成されます。
  • LA_PRECISION_MISMATCH_ERROR : Double型かFloat型のla_object_tがオブジェクトの作成の際に交じり合ってミスマッチを起こしたことを示しています。
  • LA_SINGULAR_ERROR : 行列関係を解く際に逆行列が存在しないために正しい結果を得られない事をしめしています。
  • LA_SLICE_OUT_OF_BOUNDS_ERROR : 行列、ベクトルの一部を取得する際に、取得範囲が行列の領域外まで達してしまったことを示しています。

行列、ベクトルサイズの取得

行列に対しては以下の関数で行数と列数を取得することができます。

func la_matrix_rows(matrix: la_object_t!) -> la_count_t // 行数を取得
func la_matrix_cols(matrix: la_object_t!) -> la_count_t // 列数を取得

また先ほどのログの為の拡張から例を見てみましょう。

使用例

final public var rows: UInt {
    return UInt(la_matrix_rows(self))
}
final public var cols: UInt {
    return UInt(la_matrix_cols(self))
}

また、ベクトルを表すla_object_t(行数、列数のどちらかが1な行列)に対しては次の関数で長さを取得することができます。

func la_vector_length(vector: la_object_t!) -> la_count_t // ベクトル長さを取得

行列、ベクトルの一部を取得

la_object_tで表せる行列・ベクトルの一部を取得するためにはla_matrix_sliceを用います。

func la_matrix_slice(matrix: la_object_t!, // 取得対象
                     matrix_first_row: la_index_t, // 取得初めの行
                     matrix_first_col: la_index_t, // 取得初めの列
                     matrix_row_stride: la_index_t, // 行間隔
                     matrix_col_stride: la_index_t, // 列間隔
                     slice_rows: la_count_t, // 取得する行数
                     slice_cols: la_count_t) // 取得する列数
                     -> la_object_t!

例をla_object_tに対するsubscriptで見てみましょう。

使用例

extension la_object_t {
    final public subscript(rowRange: Range<UInt>, colRange: Range<UInt>) -> la_object_t {
        return la_matrix_slice(self,
            la_index_t(rowRange.startIndex), la_index_t(colRange.startIndex), 1, 1,
            rowRange.endIndex - rowRange.startIndex,
            colRange.endIndex - colRange.startIndex)
    }
    final public subscript(row: UInt, col:UInt) -> Double {
        return la_matrix_slice(self, la_index_t(row), la_index_t(col), 0, 0, 1, 1).toArray[0]
    }
}
let matG = matA[1...2, 1...2]
println(matG.matrixDescription)
// 5.0 6.0 
// 8.0 9.0
println(matA[1, 1])
// 5.0

1つ目のsubscriptは行列・ベクトルの範囲を指定して新たなla_object_tを取得します。

2つ目のsubscriptは行列・ベクトルの要素の位置を指定して中身を取得します。

また以下のようにして間隔を空けて要素を取得することも可能です。

let vecH = la_matrix_from_double_buffer([1, 2, 3, 4, 5, 6, 7], 7, 1, 1, 
la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))
let vecI = la_matrix_slice(vecH, 0, 0, 2, 0, 4, 1)
println(vecI.matrixDescription)
// 1.0 
// 3.0 
// 5.0 
// 7.0

行間隔、列間隔は負の整数を指定することもでき、この場合は上から、左からではなく逆向きに要素を取得するようになります。上の例で第2引数を6、第4引数を-2にすると取得結果は

// 7.0 
// 5.0 
// 3.0 
// 1.0

となります。

行列からベクトルを取得

次の関数で行列の行、列を指定して行ベクトル、列ベクトルを取得することができます。

// 行ベクトルの取得
func la_vector_from_matrix_row(matrix: la_object_t!, // 取得対象行列
                               matrix_row: la_count_t) // 取得する行
                               -> la_object_t! // 行ベクトル
// 列ベクトルの取得
func la_vector_from_matrix_col(matrix: la_object_t!, // 取得対象行列
                               matrix_col: la_count_t) // 取得する列
                               -> la_object_t! // 列ベクトル

使用例

let vecK = la_vector_from_matrix_col(matA, 1)
let vecL = la_vector_from_matrix_row(matA, 1)
println(vecK.matrixDescription)
// 2.0 
// 5.0 
// 8.0 
println(vecL.matrixDescription)
// 4.0 5.0 6.0

対角行列などからベクトルを取得

la_diagonal_matrix_from_vectorの逆の操作で、対角成分やそれをずらした成分からベクトルを取得する関数が用意されています。

func la_vector_from_matrix_diagonal(matrix: la_object_t!, // 取得対象行列
                                    matrix_diagonal: la_index_t) // 対角成分からどのくらいずらすか
                                    -> la_object_t! 指定された成分から取得したベクトル

使用例

let vecM = la_vector_from_matrix_diagonal(matC, 0)
println(vecM.matrixDescription)
// 1.0 
// 2.0 
// 3.0 
let vecN = la_vector_from_matrix_diagonal(matD, 1)
println(vecN.matrixDescription)
// 1.0 
// 2.0 
// 3.0

ベクトルのノルム取得

ベクトルに対しては次の関数でノルム(長さという概念を一般化したもの)を取得可能です。

func la_norm_as_double(vector: la_object_t!, // 取得対象ベクトル
                       vector_norm: la_norm_t) // ノルムの指定
                       -> Double

この中でla_norm_tはUIntのタイプエイリアスで、次の値を用いて初期化することができます。

la_norm_t

  • LA_L1_NORM : L1ノルム(各要素絶対値の総和)を表します
  • LA_L2_NORM : L2ノルム(ユークリッド距離)を表します
  • LA_LINF_NORM : infinityノルム(各要素絶対値の最大値)を表します

使用例

let l1norm = la_norm_as_double(vecB, la_norm_t(LA_L1_NORM))
let l2norm = la_norm_as_double(vecB, la_norm_t(LA_L2_NORM))
let linorm = la_norm_as_double(vecB, la_norm_t(LA_LINF_NORM))
println(l1norm)
// 6.0
println(l2norm)
// 3.74165738677394
println(linorm)
// 3.0

単位ベクトルの取得

ノルムを指定して単位ベクトルを取得することが可能です。

func la_normalized_vector(vector: la_object_t!, // 取得対象ベクトル
                          vector_norm: la_norm_t) // ノルムを指定
                          -> la_object_t! // 単位ベクトル

使用例

let vecO = la_normalized_vector(vecB, la_norm_t(LA_L2_NORM))
// 0.267261241912424 
// 0.534522483824849 
// 0.801783725737273

各計算のための関数群

転置

行列に対して転置操作を行うことができます。行ベクトルは転置で列ベクトルになります。

func la_transpose(matrix: la_object_t!) -> la_object_t! // 転置操作

使用例

extension la_object_t {
    final public var trans : la_object_t {
        return la_transpose(self)
    }
}
let matP = matA.trans
println(matP.matrixDescription)
// 1.0 4.0 7.0 
// 2.0 5.0 8.0 
// 3.0 6.0 9.0

スカラー倍

次の関数で行うことが可能です。

func la_scale_with_double(matrix: la_object_t!, // 操作対象
                          scalar: Double) // 何倍したいか
                          -> la_object_t! // 操作結果

使用例

public func *(left: Double, right: la_object_t) -> la_object_t {
    return la_scale_with_double(right, left)
}
let matQ = 3 * matE
println(matQ.matrixDescription)
// 3.0 0.0 0.0 
// 0.0 3.0 0.0 
// 0.0 0.0 3.0

和と差

次の関数で行列、ベクトル同士の和、差を得られます。

func la_sum(obj_left: la_object_t!, obj_right: la_object_t!) -> la_object_t! // 和
func la_difference(obj_left: la_object_t!, obj_right: la_object_t!) -> la_object_t! // 差

使用例

public func +(left: la_object_t, right: la_object_t) -> la_object_t {
    return la_sum(left, right)
}
public func -(left: la_object_t, right: la_object_t) -> la_object_t {
    return la_difference(left, right)
}
let matR = matA + matE
let matS = matA - matE
println(matR.matrixDescription) 
// 2.0 2.0 3.0 
// 4.0 6.0 6.0 
// 7.0 8.0 10.0 
println(matS.matrixDescription)
// 0.0 2.0 3.0 
// 4.0 4.0 6.0 
// 7.0 8.0 8.0

要素同士の積

行列やベクトルの各要素同士をかけた要素をもつ同じサイズのla_object_tを生成できます。

func la_elementwise_product(obj_left: la_object_t!, obj_right: la_object_t!) -> la_object_t!

使用例

let matT = la_elementwise_product(matA, matE)
println(matT.matrixDescription)
// 1.0 0.0 0.0 
// 0.0 5.0 0.0 
// 0.0 0.0 9.0

内積

ベクトル同士の内積を行い、1x1のla_object_tを生成できます。

func la_inner_product(vector_left: la_object_t!, vector_right: la_object_t!) -> la_object_t!

使用例

infix operator ** {}
public func ** (left: la_object_t, right: la_object_t) -> Double {
    return la_inner_product(left, right).toArray[0]
}
let innerProductB = vecB ** vecB
println(innerProductB)
// 14.0

直積

ベクトル同士の直積を行い、各要素同士の積の結果をもった行列を生成します。

func la_outer_product(vector_left: la_object_t!, vector_right: la_object_t!) -> la_object_t!

使用例

infix operator *** {}
public func *** (left: la_object_t, right: la_object_t) -> la_object_t {
    return la_outer_product(left, right)
}
let outerProductB = vecB *** vecB
println(outerProductB.matrixDescription)
// 1.0 2.0 3.0  
// 2.0 4.0 6.0 
// 3.0 6.0 9.0

行列の積

行列やベクトルの積の結果を得ることも可能です。

func la_matrix_product(matrix_left: la_object_t!, matrix_right: la_object_t!) -> la_object_t!

使用例

public func *(left: la_object_t, right: la_object_t) -> la_object_t {
    return la_matrix_product(left, right)
}
let vecU = matA * vecB
println(vecU.matrixDescription)
// 14.0 
// 32.0 
// 50.0

関係AX=BからXを求める

A、Bを既知の行列やベクトルであるとして、それが関係: AX=B (AXは積)を満たしている時にXを求めたいことがあります。

このXを次の関数を用いて簡単に求めることが可能です。

func la_solve(matrix_system: la_object_t!, // AX = BのAに当たる行列やベクトル 
              obj_rhs: la_object_t!) // AX = BのBに当たる行列やベクトル
              -> la_object_t! // AX = B のXに当たる行列やベクトル

Aが逆行列を持たない場合など、計算がうまくいかない場合は最終的にla_status_tの取得結果にエラーとして現れます。

使用例

let matV = la_solve(matC, matA)
println(matV.matrixDescription)
// 1.0 2.0 3.0 
// 2.0 2.5 3.0 
// 2.33333333333333 2.66666666666667 3.0

splat

計算のための便利なla_object_tとしてLinearAlgebraライブラリにはsplatという概念があります。

splatオブジェクトはすべての要素が等しい行列やベクトルを表し、各演算の対象によってサイズを自動的に調整してくれます。

サイズの推測は各演算のもう一つのla_object_tによって次のように定まります。

演算  推測されるサイズ (行数、列数)
la_sum(A, splat) (Aの行数、Aの列数)
la_sum(splat, A) (Aの行数、Aの列数)
la_difference(A, splat) (Aの行数、Aの列数)
la_difference(splat, A) (Aの行数、Aの列数)
la_elementwise_product(A, splat) (Aの行数、Aの列数)
la_elementwise_product(splat, A) (Aの行数、Aの列数)
la_inner_product(A, splat) (Aの長さ, 1)
la_inner_product(splat, A) (1, Aの長さ)
la_matrix_product(A, splat) (Aの列数, 1)
la_matrix_product(splat, A) (1, Aの行数)

この仕組の為に、演算について複数の対象がsplatであることは認められません。

splatの生成

次の関数を用いてsplatを生成できます。

func la_splat_from_float(scalar_value: Float, // splat要素のFloat値
                         attributes: la_attribute_t) 
                         -> la_object_t!

func la_splat_from_double(scalar_value: Double, // splat要素のDouble値
                          attributes: la_attribute_t) 
                          -> la_object_t!

func la_splat_from_vector_element(vector: la_object_t!, // splatの生成に用いるベクトル
                                  vector_index: la_index_t) // どの要素をsplatの要素として用いるか指定
                                  -> la_object_t!

func la_splat_from_matrix_element(matrix: la_object_t!, // splatの生成に用いる行列
                                  matrix_row: la_index_t, // どの行、列の要素を 
                                  matrix_col: la_index_t) // splatの要素として用いるか指定
                                  -> la_object_t!

使用例

let splatA = la_splat_from_double(2.0, la_attribute_t(LA_DEFAULT_ATTRIBUTES))
let splatB = la_splat_from_vector_element(vecB, 1)
let vecW = vecB + splatA + splatB
println(vecW.matrixDescription)
// 5.0 
// 6.0 
// 7.0     
let splatC = la_splat_from_matrix_element(matA, 1, 1)
let matX = matA + splatB
println(matX.matrixDescription)
// 6.0 7.0 8.0 
// 9.0 10.0 11.0 
// 12.0 13.0 14.0

splatから指定サイズのベクトル、行列を生成

時にはsplatから一定の大きさのベクトルや行列を取得したい場合もあると思います。

そのような場合に合わせた関数も提供されています。

func la_vector_from_splat(splat: la_object_t!, // splat
                          vector_length: la_count_t) // 取得したいベクトルの長さ
                          -> la_object_t!

func la_matrix_from_splat(splat: la_object_t!, // splat
                          matrix_rows: la_count_t, // 取得したい行列の行数
                          matrix_cols: la_count_t) // 取得したい行列の列数
                          -> la_object_t!

使用例

let splatA = la_splat_from_double(2.0, la_attribute_t(LA_DEFAULT_ATTRIBUTES))
let vecY = la_vector_from_splat(splatA, 6)
let matZ = la_matrix_from_splat(splatA, 2, 3)
println(vecY.matrixDescription)
// 2.0 
// 2.0 
// 2.0 
// 2.0 
// 2.0 
// 2.0
println(matZ.matrixDescription)
// 2.0 2.0 2.0 
// 2.0 2.0 2.0

参考サイト