強化学習用のニューラルネットワークをSwiftで書いてみた。（その3）

昨日は強化学習用のニューラルネットワークの計算を行列で表現した。

今日はそれを使って実際に実装していく。

なお、Swiftでの行列演算については、以下を参照：

ここで定義したMatrixクラス、Vectorクラスを使っていくことになる。

また、Rubyでの実装は以下を参照：
（行列計算を使ってないので、読むのがかなり大変だけど・・・）

ニューラルネットワークの仕様は、上記のRubyのものと同じにする。

ValueNetworkプロトコル

まず、強化学習の関数近似に使うのがニューラルネットワークでもHMEでもいいように、インタフェースとなるプロトコルを定義する：

//==============================
// ValueNetwork
//------------------------------
// ValueNetwork.swift
//==============================

import Foundation

protocol ValueNetwork {
  func getValue(input: Vector) -> Double
  func getValueAndWeightGradient(input: Vector) -> (Double, Weight)
  func getValue(input: Vector, withWeightDiff weightDiff: Weight, scale: Double) -> Double
  func addWeight(weightDiff: Weight)
}

定義の中でWeightという型が使われてるけど、これは次に定義するプロトコルで、重みの内部構成を隠蔽するためのインタフェース。

Weightプロトコル

ということで、Weightプロトコルの定義：

//==============================
// ValueNetwork
//------------------------------
// Weight.swift
//==============================

import Foundation

protocol Weight {
  func scale(scalar: Double) -> Weight
  func add(other: Weight) -> Weight
  func subtract(other: Weight) -> Weight
}

func *(left: Weight, right: Double) -> Weight {
  return left.scale(right)
}

func *(left: Double, right: Weight) -> Weight {
  return right.scale(left)
}

func +(left: Weight, right: Weight) -> Weight {
  return left.add(right)
}

func -(left: Weight, right: Weight) -> Weight {
  return left.subtract(right)
}

重みに対して行いたいことは、足し算（と引き算）と掛け算なので、それらをインタフェースとして定義している。
また、記述が簡単になるように、演算子のオーバーロードも行っている。

ValueNNクラス

そして、強化学習用のニューラルネットワークの実装。
ValueNNクラスとして実装した：

//==============================
// ValueNetwork
//------------------------------
// ValueNN.swift
//==============================

import Foundation

class ValueNN: ValueNetwork {
  class NNWeight: Weight {
    private let hiddenLayerWeight: Matrix
    private let hiddenLayerBias: Vector
    private let outputLayerWeight: Vector
    private let outputLayerBias: Double
    
    private init(hiddenLayerWeight: Matrix,
                 hiddenLayerBias: Vector,
                 outputLayerWeight: Vector,
                 outputLayerBias: Double) {
      self.hiddenLayerWeight = hiddenLayerWeight
      self.hiddenLayerBias = hiddenLayerBias
      self.outputLayerWeight = outputLayerWeight
      self.outputLayerBias = outputLayerBias
    }
    
    func scale(scalar: Double) -> Weight {
      let hiddenLayerWeight = self.hiddenLayerWeight * scalar
      let hiddenLayerBias = self.hiddenLayerBias * scalar
      let outputLayerWeight = self.outputLayerWeight * scalar
      let outputLayerBias = self.outputLayerBias * scalar
      return NNWeight(hiddenLayerWeight: hiddenLayerWeight,
                      hiddenLayerBias: hiddenLayerBias,
                      outputLayerWeight: outputLayerWeight,
                      outputLayerBias: outputLayerBias)
    }
    
    func add(other: Weight) -> Weight {
      let otherNNWeight = other as! NNWeight
      let hiddenLayerWeight = self.hiddenLayerWeight + otherNNWeight.hiddenLayerWeight
      let hiddenLayerBias = self.hiddenLayerBias + otherNNWeight.hiddenLayerBias
      let outputLayerWeight = self.outputLayerWeight + otherNNWeight.outputLayerWeight
      let outputLayerBias = self.outputLayerBias + otherNNWeight.outputLayerBias
      return NNWeight(hiddenLayerWeight: hiddenLayerWeight,
                      hiddenLayerBias: hiddenLayerBias,
                      outputLayerWeight: outputLayerWeight,
                      outputLayerBias: outputLayerBias)
    }
    
    func subtract(other: Weight) -> Weight {
      let otherNNWeight = other as! NNWeight
      let hiddenLayerWeight = self.hiddenLayerWeight - otherNNWeight.hiddenLayerWeight
      let hiddenLayerBias = self.hiddenLayerBias - otherNNWeight.hiddenLayerBias
      let outputLayerWeight = self.outputLayerWeight - otherNNWeight.outputLayerWeight
      let outputLayerBias = self.outputLayerBias - otherNNWeight.outputLayerBias
      return NNWeight(hiddenLayerWeight: hiddenLayerWeight,
                      hiddenLayerBias: hiddenLayerBias,
                      outputLayerWeight: outputLayerWeight,
                      outputLayerBias: outputLayerBias)
    }
  }
  
  private static let activationNormalGradient: Double = 1.0
  private static let activationLesserGradient: Double = 0.1
  
  private var weight: NNWeight
  let outputMin: Double
  let outputMax: Double
  
  init(inputSize: Int, hiddenUnitSize: Int, outputMin: Double, outputMax: Double) {
    self.outputMin = outputMin
    self.outputMax = outputMax
    
    let hiddenLayerWeightVariance = 1.0 / (Double(inputSize) + 1.0)
    let hiddenLayerWeightGenerator = NormalDistRandom.init(expected: 0.0, variance: hiddenLayerWeightVariance)
    
    let hiddenLayerWeightBuffer = (0..<hiddenUnitSize).map { _ in
      return (0..<inputSize).map { _ in
        return hiddenLayerWeightGenerator.getRandom()
      }
    }
    let hiddenLayerWeight = Matrix.fromArray(hiddenLayerWeightBuffer)
    
    let hiddenLayerBiasBuffer = (0..<hiddenUnitSize).map { _ in
      return hiddenLayerWeightGenerator.getRandom()
    }
    let hiddenLayerBias = Vector.fromArray(hiddenLayerBiasBuffer)
    
    let outputLayerWeightVariance = 1.0 / (Double(hiddenUnitSize) + 1.0)
    let outputLayerWeightGenerator = NormalDistRandom(expected: 0.0, variance: outputLayerWeightVariance)
    
    let outputLayerWeightBuffer = (0..<hiddenUnitSize).map { _ in
      return outputLayerWeightGenerator.getRandom()
    }
    let outputLayerWeight = Vector.fromArray(outputLayerWeightBuffer)
    
    let outputLayerBias = outputLayerWeightGenerator.getRandom()
    
    self.weight = NNWeight(hiddenLayerWeight: hiddenLayerWeight,
                           hiddenLayerBias: hiddenLayerBias,
                           outputLayerWeight: outputLayerWeight,
                           outputLayerBias: outputLayerBias)
  }
  
  func getValue(input: Vector) -> Double {
    let hiddenLayerWeightedInput = ((self.weight.hiddenLayerWeight * input) as! Vector) + self.weight.hiddenLayerBias
    let hiddenLayerOutput = hiddenLayerWeightedInput.map {
      [unowned self] (weightedInput: Double) in
      return self.hiddenLayerOutputForWeightedInput(weightedInput)
    }
    
    let outputLayerWeightedInput = self.weight.outputLayerWeight +* hiddenLayerOutput + self.weight.outputLayerBias
    let outputLayerOutput = self.outputLayerOutputForWeightedInput(outputLayerWeightedInput)

    return outputLayerOutput
  }
  
  func getValueAndWeightGradient(input: Vector) -> (Double, Weight) {
    let hiddenLayerWeightedInput = ((self.weight.hiddenLayerWeight * input) as! Vector) + self.weight.hiddenLayerBias
    let hiddenLayerOutput = hiddenLayerWeightedInput.map {
      [unowned self] (weightedInput: Double) in
      return self.hiddenLayerOutputForWeightedInput(weightedInput)
    }
    let hiddenLayerGradient = hiddenLayerWeightedInput.map {
      [unowned self] (weightedInput: Double) in
      return self.hiddenLayerGradientForWeightedInput(weightedInput)
    }
    
    let outputLayerWeightedInput = self.weight.outputLayerWeight +* hiddenLayerOutput + self.weight.outputLayerBias
    let outputLayerOutput = self.outputLayerOutputForWeightedInput(outputLayerWeightedInput)
    let outputLayerGradient = self.outputLayerGradientForWeightedInput(outputLayerWeightedInput)
    
    let outputLayerDelta = outputLayerGradient
    let hiddenLayerDelta = outputLayerDelta * self.weight.outputLayerWeight <*> hiddenLayerGradient
    
    let hiddenLayerWeightGradient = hiddenLayerDelta *+ input
    let hiddenLayerBiasGradient = hiddenLayerDelta
    let outputLayerWeightGradient = outputLayerDelta * hiddenLayerOutput
    let outputLayerBiasGradint = outputLayerDelta
    let weightGradient = NNWeight(hiddenLayerWeight: hiddenLayerWeightGradient,
                                  hiddenLayerBias: hiddenLayerBiasGradient,
                                  outputLayerWeight: outputLayerWeightGradient,
                                  outputLayerBias: outputLayerBiasGradint)
    
    return (outputLayerOutput, weightGradient)
  }
  
  func getValue(input: Vector, withWeightDiff weightDiff: Weight, scale: Double) -> Double {
    let newWeight = (self.weight + weightDiff * scale) as! NNWeight
    
    let hiddenLayerWeightedInput = ((newWeight.hiddenLayerWeight * input) as! Vector) + newWeight.hiddenLayerBias
    let hiddenLayerOutput = hiddenLayerWeightedInput.map {
      [unowned self] (weightedInput: Double) in
      return self.hiddenLayerOutputForWeightedInput(weightedInput)
    }
    
    let outputLayerWeightedInput = newWeight.outputLayerWeight +* hiddenLayerOutput + newWeight.outputLayerBias
    let outputLayerOutput = self.outputLayerOutputForWeightedInput(outputLayerWeightedInput)
    
    return outputLayerOutput
  }
  
  func addWeight(weightDiff: Weight) {
    self.weight = (self.weight + weightDiff) as! NNWeight
  }
  
  private func hiddenLayerOutputForWeightedInput(weightedInput: Double) -> Double {
    if weightedInput >= 0.0 {
      return ValueNN.activationNormalGradient * weightedInput
    } else {
      return ValueNN.activationLesserGradient * weightedInput
    }
  }
  
  private func hiddenLayerGradientForWeightedInput(weightedInput: Double) -> Double {
    if weightedInput >= 0.0 {
      return ValueNN.activationNormalGradient
    } else {
      return ValueNN.activationLesserGradient
    }
  }
  
  private func outputLayerOutputForWeightedInput(weightedInput: Double) -> Double {
    if weightedInput < self.outputMin {
      return (ValueNN.activationLesserGradient * weightedInput
                + (ValueNN.activationNormalGradient - ValueNN.activationLesserGradient) * self.outputMin)
    } else if weightedInput < self.outputMax {
      return ValueNN.activationNormalGradient * weightedInput
    } else {
      return (ValueNN.activationLesserGradient * weightedInput
                + (ValueNN.activationNormalGradient - ValueNN.activationLesserGradient) * self.outputMax)
    }
  }
  
  private func outputLayerGradientForWeightedInput(weightedInput: Double) -> Double {
    if (weightedInput < self.outputMin) || (self.outputMax < weightedInput) {
      return ValueNN.activationLesserGradient
    } else {
      return ValueNN.activationNormalGradient
    }
  }
}

細かく説明しないけど、強化学習用のニューラルネットワークをSwiftで書いてみた。（その2） - いものやま。に書いた行列計算をそのまま実装している。
やっている計算はRubyのものと一緒なんだけど、Matrixクラス、Vectorクラスを定義してあるので、（比較的）スッキリした実装になっている。

動作確認

動作確認として、Rubyのときと同様に、以下のコードを書いた：

//==============================
// ValueNetwork
//------------------------------
// main.swift
//
// Test code for ValueNetwork
//==============================

import Foundation

// ValueNN

let valueNN = ValueNN(inputSize: 3, hiddenUnitSize: 10, outputMin: -1.0, outputMax: 1.0)

let inputMatrix = Matrix.fromArray([[1.0, 1.0, 1.0],
                                    [1.0, 1.0, 0.0],
                                    [1.0, 0.0, 1.0],
                                    [0.0, 1.0, 1.0],
                                    [1.0, 0.0, 0.0],
                                    [0.0, 1.0, 0.0],
                                    [0.0, 0.0, 1.0]]).transpose()

for col in (0..<inputMatrix.col) {
  let input = inputMatrix.colVector(col)
  let (output, weightGradient) = valueNN.getValueAndWeightGradient(input)
  let outputWithWeightGradient = valueNN.getValue(input, withWeightDiff: weightGradient, scale: 1.0)
  var diff = outputWithWeightGradient - output
  
  var scale = 1.0
  var upperBound: Double! = nil
  var lowerBound: Double! = nil
  var last = 100
  for i in (0..<100) {
    if (diff < 0.0) || (1.1 < diff) {
      upperBound =  scale
      scale = (lowerBound == nil) ? scale / 2.0 : (upperBound + lowerBound) / 2.0
    } else if diff < 0.9 {
      lowerBound = scale
      scale = (upperBound == nil) ? scale * 2.0 : (upperBound + lowerBound) / 2.0
    } else {
      last = i
      break
    }
    
    let outputWithScaledWeightGradient = valueNN.getValue(input, withWeightDiff: weightGradient, scale: scale)
    diff = outputWithScaledWeightGradient - output
  }
  
  let outputWith01Scaled = valueNN.getValue(input, withWeightDiff: weightGradient, scale: 0.1 * scale)
  let diffWith01Scaled = outputWith01Scaled - output
  
  print("input: \(input.transpose()), output: \(output)")
  print("  scale: \(scale), iterations: \(last)")
  print("  diff (scaled): \(diff), diff (0.1*scaled): \(diffWith01Scaled)")
  
  let weightDiff = weightGradient * 0.1 * scale
  valueNN.addWeight(weightDiff)
  let newOutput = valueNN.getValue(input)
  let newDiff = newOutput - output
  print("  new output: \(newOutput), diff: \(newDiff)")
}

// 以下略

実行例は、以下：

----------
ValueNN
----------
input: [1.0, 1.0, 1.0], output: 0.131021133192784
  scale: 0.125, iterations: 3
  diff (scaled): 0.939252246852243, diff (0.1*scaled): 0.155481391989248
  new output: 0.286502525182032, diff: 0.155481391989248
input: [1.0, 1.0, 0.0], output: 0.363438002291251
  scale: 0.5, iterations: 1
  diff (scaled): 0.92473632645851, diff (0.1*scaled): 0.331184714810343
  new output: 0.694622717101594, diff: 0.331184714810343
input: [1.0, 0.0, 1.0], output: 0.0041779249095818
  scale: 0.125, iterations: 3
  diff (scaled): 1.03338785763323, diff (0.1*scaled): 0.13662619759527
  new output: 0.140804122504852, diff: 0.13662619759527
input: [0.0, 1.0, 1.0], output: 0.257684625662284
  scale: 0.5, iterations: 1
  diff (scaled): 0.992301758621281, diff (0.1*scaled): 0.335977958448971
  new output: 0.593662584111254, diff: 0.335977958448971
input: [1.0, 0.0, 0.0], output: 0.184531558718531
  scale: 0.5, iterations: 1
  diff (scaled): 1.01382568012975, diff (0.1*scaled): 0.27650775017151
  new output: 0.461039308890041, diff: 0.27650775017151
input: [0.0, 1.0, 0.0], output: 0.166668276265283
  scale: 0.5, iterations: 1
  diff (scaled): 0.974321711025463, diff (0.1*scaled): 0.217514514547388
  new output: 0.384182790812671, diff: 0.217514514547388
input: [0.0, 0.0, 1.0], output: 0.180236823017496
  scale: 0.5, iterations: 1
  diff (scaled): 0.945552060789391, diff (0.1*scaled): 0.200679893784258
  new output: 0.380916716801753, diff: 0.200679893784258
# 以下略

今日はここまで！

強化学習用のニューラルネットワークをSwiftで書いてみた。（その2）

昨日は乱数生成器の実装を行った。

今日は強化学習用のニューラルネットワークの計算を行列で表現する。

強化学習用のニューラルネットワークの計算

説明を簡単にするために、ここでは次のようなニューラルネットワークを考える：

• 3層ニューラルネットワーク
• 入力は 、出力は
• 中間層のユニット数は
  • 各ユニットの重みは 、バイアスは
  • 活性化関数は
• 出力層のユニット数は1
  • ユニットの重みは 、バイアスは
  • 活性化関数は

行列で表現しない計算

まず、おさらいとして、行列で表現していない計算を書いておく：

1. 入力 から、中間層の出力 を求める：
   1. 
   2. 
2. 中間層の出力 から、出力層の出力 を求める：
   1. 
   2. 
3. 出力層のデルタ を求める：
   1. 
4. 中間層のデルタ を求める：
   1. 
5. 各偏微分を求める：
   1. 
   2. 
   3. 
   4. 

行列で表現した計算

ここで、中間層の重みを で次のように表すことにする：

そして、表記をもう少しスマートにすると、上の計算は次のように書き直すことが出来る：

1. 入力 から、中間層の出力 を求める：
   1. 
   2. 
2. 中間層の出力 から、出力層の出力 を求める：
   1. 
   2. 
3. 出力層のデルタ を求める：
   1. 
4. 中間層のデルタ を求める：
   1. 
5. 各偏微分を求める：
   1. 
   2. 
   3. 
   4. 

なお、表記法として、関数 の引数にベクトル（や行列）が来ている場合、各要素に関数 を適用したベクトル（や行列）を返すものとする。
また、 は行列（ベクトルを含む）のアダマール積（要素ごとの積）、 はベクトルの直積（外積）を表すものとする。

ここまでくれば、Swiftでの行列計算について調べてみた。（その4） - いものやま。で作ったクラスを使って計算が出来るようになる。

中間層や出力の次元が増えた場合

おまけとして、中間層が増えた場合や、出力の次元が増えた場合の計算も書いておく。

以下では、各層 の重みが 、バイアスが 、活性化関数が とする。
（第1層は入力層、第 層は出力層）

1. 入力層の出力 を とする。
2. 第 層 の出力 を求める：
   1. 
   2. 
3. 出力 を とする。
4. 出力層のデルタ を求める：
   1. 
5. 第 層 のデルタ を求める：
   1. 
6. 各偏微分 を求める：
   1. 
   2. 

今日はここまで！

強化学習用のニューラルネットワークをSwiftで書いてみた。（その1）

強化学習の関数近似にニューラルネットワークを組合せるということをやってきていた。

強化学習については以下から：

ニューラルネットワークについては以下から：

複数のニューラルネットワークを組合せるHME（Hierarchical Mixtures of Experts）については以下から：

強化学習とニューラルネットワーク（HME含む）の組合せについては以下から：

この強化学習用のニューラルネットワークはRubyで書いてたけど、iOSアプリを作るとなると、Swiftに移植しないといけない。
ということで、Swiftに移植してみた。

乱数生成器

まず必要となるのが、乱数生成器。

一様乱数の生成器は以前も書いたけど、再掲。

//==============================
// Random
//------------------------------
// Random.swift
//==============================

import Foundation
import Security

class Random {
  static func getUniformedRandom(upperBound: Int) -> Int {
    let upper = UInt32(upperBound)
    let ignoreRegion = 0xffff_ffff - 0xffff_ffff % upper

    while true {
      let buffer = UnsafeMutablePointer<UInt32>.alloc(1)
      SecRandomCopyBytes(kSecRandomDefault, 4, UnsafeMutablePointer<UInt8>(buffer))
      let randomValue = buffer.memory
      buffer.dealloc(1)
      if randomValue < ignoreRegion {
        return Int(randomValue % upper)
      }
    }
  }

  static func getRandomProbability() -> Double {
    return Double(Random.getUniformedRandom(0x1_0000)) / Double(0xffff)
  }

  private init() {}
}

extension Array {
  func sample() -> Element {
    let index = Random.getUniformedRandom(self.count)
    return self[index]
  }
}

この一様乱数の生成器の説明は、以下から：

そして、ニューラルネットワークの重みの初期化で正規分布に従う乱数生成器も必要になるので、移植。

正規分布に従う乱数生成器の説明は、以下から：

実装は以下：

//==============================
// Random
//------------------------------
// NormalDistRandom.swift
//==============================

import Foundation

class NormalDistRandom {
  private let expected: Double
  private let variance: Double
  private var values: [Double]
  
  init(expected: Double, variance: Double) {
    self.expected = expected
    self.variance = variance
    self.values = []
  }
  
  func getRandom() -> Double {
    if self.values.count == 0 {
      var a = 0.0
      var b = 0.0
      repeat {
        a = Random.getRandomProbability()
        b = Random.getRandomProbability()
      } while (a == 0.0) || (a == 1.0) || (b == 0.0) || (b == 1.0)
      
      let z1 = sqrt(-2.0 * log(a)) * cos(2.0 * M_PI * b)
      let z2 = sqrt(-2.0 * log(a)) * sin(2.0 * M_PI * b)
      
      let rand1 = z1 * sqrt(self.variance) + self.expected
      let rand2 = z2 * sqrt(self.variance) + self.expected
      
      self.values = [rand1, rand2]
    }
    
    return self.values.removeLast()
  }
}

乱数生成器の動作確認

これらの乱数生成器の動作確認として、以下のコードを書いた：

//==============================
// Random
//------------------------------
// main.swift
//
// Test code for Random
//==============================

import Foundation

let count = 1000

// Random.getUniformedRandom()

var valueCount = [
  0: 0, 1: 0, 2: 0, 3: 0, 4: 0,
  5: 0, 6: 0, 7: 0, 8: 0, 9: 0,
]
for _ in 0..<count {
  let randomValue = Random.getUniformedRandom(10)
  valueCount[randomValue] = valueCount[randomValue]! + 1
}

print("----------")
print("Random.getUniformedRandom()")
print("----------")
for value in valueCount.keys.sort() {
  print("\(value): \(valueCount[value]!)")
}

// Random.getRandomProbability()

var sum = 0.0
var squareSum = 0.0
var minValue = 100.0
var maxValue = -100.0
for i in 0..<count {
  let randomValue = Random.getRandomProbability()
  sum += randomValue
  squareSum += randomValue * randomValue
  if randomValue < minValue {
    minValue = randomValue
  }
  if maxValue < randomValue {
    maxValue = randomValue
  }
}
var mean = sum / Double(count)
var variance = squareSum / Double(count) - mean * mean

print("----------")
print("Random.getRandomProbability()")
print("----------")
print("min: \(minValue)")
print("max: \(maxValue)")
print("mean: \(mean)")
print("variance: \(variance)")

// Array#sample()

let array = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

valueCount = [
  0: 0, 1: 0, 2: 0, 3: 0, 4: 0,
  5: 0, 6: 0, 7: 0, 8: 0, 9: 0,
]
for _ in 0..<count {
  let randomValue = array.sample()
  valueCount[randomValue] = valueCount[randomValue]! + 1
}

print("----------")
print("Array#sample()")
print("----------")
for value in valueCount.keys.sort() {
  print("\(value): \(valueCount[value]!)")
}

// NormalDistRandom

let normalDistRandom = NormalDistRandom(expected: 1.0, variance: 2.0)

sum = 0.0
squareSum = 0.0
minValue = 100.0
maxValue = -100.0
for i in 0..<count {
  let randomValue = normalDistRandom.getRandom()
  sum += randomValue
  squareSum += randomValue * randomValue
  if randomValue < minValue {
    minValue = randomValue
  }
  if maxValue < randomValue {
    maxValue = randomValue
  }
}
mean = sum / Double(count)
variance = squareSum / Double(count) - mean * mean

print("----------")
print("normalDistRandom.getRanodm()")
print("----------")
print("min: \(minValue)")
print("max: \(maxValue)")
print("mean: \(mean)")
print("variance: \(variance)")

実行例は、以下：

----------
Random.getUniformedRandom()
----------
0: 98
1: 104
2: 100
3: 85
4: 105
5: 112
6: 117
7: 86
8: 103
9: 90
----------
Random.getRandomProbability()
----------
min: 9.15541313801785e-05
max: 0.998825055313954
mean: 0.472300785839628
variance: 0.0830635419408054
----------
Array#sample()
----------
0: 90
1: 110
2: 100
3: 112
4: 98
5: 91
6: 101
7: 102
8: 106
9: 90
----------
normalDistRandom.getRanodm()
----------
min: -3.76411577670718
max: 5.42937254648993
mean: 0.989138466448474
variance: 2.1180087175273

今日はここまで！

Swiftでの行列計算について調べてみた。（まとめ）

これまでの各記事は以下から。

• 行列、ベクトル、Splatの生成
  • Swiftでの行列計算について調べてみた。（その1） - いものやま。
  • Swiftでの行列計算について調べてみた。（その2） - いものやま。
• 行列、ベクトルの演算
  • Swiftでの行列計算について調べてみた。（その3） - いものやま。
• ラッパークラスの実装
  • Swiftでの行列計算について調べてみた。（その4） - いものやま。

特に書くことがないw

とりあえず、これでSwiftで行列計算が出来るようになったので、関数近似のためのニューラルネットワークの実装をSwiftでもやってみたい。

今日はここまで！

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

昨日は行列とベクトルの演算について説明した。

今日はもうちょっとスマートな書き方が出来るようにするために、ラッパークラスを作っていく。

Matrixクラス

まず、行列を表すMatrixクラスの実装。

なお、途中に出てくるVectorクラスはMatrixクラスを継承したもので、ベクトルを表す。
詳細は後述。

それと、Matrixクラスはホントはla_object_tを継承して作りたかったんだけど、la_object_tがNSObjectProtocolに適合していなかったので、プロパティとして持つようにしている。

//==============================
// LinearAlgebra
//------------------------------
// Matrix.swift
//==============================

import Foundation
import Accelerate

class Matrix: CustomStringConvertible {
  class func fromArray(array: [[Double]]) -> Matrix {
    let row = array.count
    let col = array[0].count
    let buffer = Array(array.flatten())
    let linearAlgebraObject = la_matrix_from_double_buffer(buffer,
                                                           la_count_t(row), la_count_t(col),
                                                           la_count_t(col),
                                                           la_hint_t(LA_NO_HINT),
                                                           la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  class func identityMatrix(size: Int) -> Matrix {
    let linearAlgebraObject = la_identity_matrix(la_count_t(size),
                                                 la_scalar_type_t(LA_SCALAR_TYPE_DOUBLE),
                                                 la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  class func diagonalMatrixFromVector(vector: Vector, offset: Int = 0) -> Matrix {
    let linearAlgebraObject = la_diagonal_matrix_from_vector(vector.linearAlgebraObject, offset)
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  class func filledWith(value: Double, row: Int, col: Int) -> Matrix {
    let splat = la_splat_from_double(value, la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    let linearAlgebraObject = la_matrix_from_splat(splat, la_count_t(row), la_count_t(col))
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  private class func createMatrixOrVector(linearAlgebraObject: la_object_t) -> Matrix {
    var object = Matrix(linearAlgebraObject: linearAlgebraObject)
    if object.isVector {
      object = Vector(linearAlgebraObject: linearAlgebraObject)
    }
    return object
  }
  
  let linearAlgebraObject: la_object_t
  
  init(linearAlgebraObject: la_object_t) {
    self.linearAlgebraObject = linearAlgebraObject
  }
  
  private(set) lazy var row: Int = {
    [unowned self] in
    return Int(la_matrix_rows(self.linearAlgebraObject))
  }()
  
  private(set) lazy var col: Int = {
    [unowned self] in
    return Int(la_matrix_cols(self.linearAlgebraObject))
  }()
  
  var isVector: Bool {
    return (self.row == 1) || (self.col == 1)
  }
  
  private lazy var buffer: [Double] = {
    [unowned self] in
    var buffer = [Double](count: self.row * self.col, repeatedValue: 0.0)
    la_matrix_to_double_buffer(&buffer, la_count_t(self.col), self.linearAlgebraObject)
    return buffer
  }()
  
  subscript(row: Int, col: Int) -> Double {
    return self.buffer[row * self.col + col]
  }
  
  func toArray() -> [[Double]] {
    return (0..<self.row).map { (row: Int) -> [Double] in
      let from = self.col * row
      let to = self.col * (row + 1)
      return Array(self.buffer[from..<to])
    }
  }
  
  var description: String {
    return self.toArray().map({$0.description}).joinWithSeparator("\n")
  }
  
  var errorCode: Int {
    return la_status(self.linearAlgebraObject)
  }
  
  var hasError: Bool {
    return self.errorCode < 0
  }
  
  func slice(firstRow: Int, _ firstCol: Int,
             _ rowStep: Int, _ colStep: Int,
             _ row: Int, _ col: Int) -> Matrix {
    let linearAlgebraObject = la_matrix_slice(self.linearAlgebraObject,
                                              firstRow, firstCol,
                                              rowStep, colStep,
                                              la_count_t(row), la_count_t(col))
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  func rowVector(row: Int) -> Vector {
    let linearAlgebraObject = la_vector_from_matrix_row(self.linearAlgebraObject, la_count_t(row))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func colVector(col: Int) -> Vector {
    let linearAlgebraObject = la_vector_from_matrix_col(self.linearAlgebraObject, la_count_t(col))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func diagonalVector(offset: Int = 0) -> Vector {
    let linearAlgebraObject = la_vector_from_matrix_diagonal(self.linearAlgebraObject, offset)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func transpose() -> Matrix {
    let linearAlgebraObject = la_transpose(self.linearAlgebraObject)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  func scale(scalar: Double) -> Matrix {
    let linearAlgebraObject = la_scale_with_double(self.linearAlgebraObject, scalar)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  func sum(other: Matrix) -> Matrix {
    let linearAlgebraObject = la_sum(self.linearAlgebraObject, other.linearAlgebraObject)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  func difference(other: Matrix) -> Matrix {
    let linearAlgebraObject = la_difference(self.linearAlgebraObject, other.linearAlgebraObject)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  func matrixProduct(other: Matrix) -> Matrix {
    let linearAlgebraObject = la_matrix_product(self.linearAlgebraObject, other.linearAlgebraObject)
    return Matrix.createMatrixOrVector(linearAlgebraObject)
  }
  
  func elementwiseProduct(other: Matrix) -> Matrix {
    let linearAlgebraObject = la_elementwise_product(self.linearAlgebraObject, other.linearAlgebraObject)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  func map(transform: (Double) -> Double) -> Matrix {
    let buffer = self.buffer.map(transform)
    let linearAlgebraObject = la_matrix_from_double_buffer(buffer,
                                                           la_count_t(self.row), la_count_t(self.col),
                                                           la_count_t(self.col),
                                                           la_hint_t(LA_NO_HINT),
                                                           la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
}

func +(left: Matrix, right: Matrix) -> Matrix {
  return left.sum(right)
}

func -(left: Matrix, right: Matrix) -> Matrix {
  return left.difference(right)
}

func *(left: Matrix, right: Double) -> Matrix {
  return left.scale(right)
}

func *(left: Double, right: Matrix) -> Matrix {
  return right.scale(left)
}

func *(left: Matrix, right: Matrix) -> Matrix {
  return left.matrixProduct(right)
}

infix operator <*> {
precedence 150
associativity left
}
func <*>(left: Matrix, right: Matrix) -> Matrix {
  return left.elementwiseProduct(right)
}

なお、新たに定義している演算子<*>は、行列のアダマール積（要素ごとの積）のための演算子。

Vectorクラス

次に、ベクトルを表すVectorクラスの実装。

//==============================
// LinearAlgebra
//------------------------------
// Vector.swift
//
// by ouka-do
//==============================

import Foundation
import Accelerate

class Vector: Matrix {
  class func fromArray(array: [Double]) -> Vector {
    let linearAlgebraObject = la_matrix_from_double_buffer(array,
                                                           la_count_t(array.count), la_count_t(1),
                                                           la_count_t(1),
                                                           la_hint_t(LA_NO_HINT),
                                                           la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  class func filledWith(value: Double, size: Int) -> Vector {
    let splat = la_splat_from_double(value, la_attribute_t(LA_DEFAULT_ATTRIBUTES))
    let linearAlgebraObject = la_vector_from_splat(splat, la_count_t(size))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  private(set) lazy var size: Int = {
    [unowned self] in
    return Int(la_vector_length(self.linearAlgebraObject))
  }()
  
  subscript(i: Int) -> Double {
    return self[i, 0]
  }
  
  func toArray() -> [Double] {
    let array = super.toArray()
    return Array(array.flatten())
  }
  
  func slice(first: Int, _ step: Int, _ size: Int) -> Vector {
    let linearAlgebraObject = la_vector_slice(self.linearAlgebraObject, first, step, la_count_t(size))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func norm(type: Int) -> Double {
    return la_norm_as_double(self.linearAlgebraObject, la_norm_t(type))
  }
  
  func normalizedVector(type: Int) -> Vector {
    let linearAlgebraObject = la_normalized_vector(self.linearAlgebraObject, la_norm_t(type))
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  override func transpose() -> Vector {
    let linearAlgebraObject = la_transpose(self.linearAlgebraObject)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  override func scale(scalar: Double) -> Vector {
    let linearAlgebraObject = la_scale_with_double(self.linearAlgebraObject, scalar)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func sum(other: Vector) -> Vector {
    let linearAlgebraObject = la_sum(self.linearAlgebraObject, other.linearAlgebraObject)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func difference(other: Vector) -> Vector {
    let linearAlgebraObject = la_difference(self.linearAlgebraObject, other.linearAlgebraObject)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func elementwiseProduct(other: Vector) -> Vector {
    let linearAlgebraObject = la_elementwise_product(self.linearAlgebraObject, other.linearAlgebraObject)
    return Vector(linearAlgebraObject: linearAlgebraObject)
  }
  
  func innerProduct(other: Vector) -> Double {
    let linearAlgebraObject = la_inner_product(self.linearAlgebraObject, other.linearAlgebraObject)
    let vector = Vector(linearAlgebraObject: linearAlgebraObject)
    return vector[0]
  }
  
  func outerProduct(other: Vector) -> Matrix {
    let linearAlgebraObject = la_outer_product(self.linearAlgebraObject, other.linearAlgebraObject)
    return Matrix(linearAlgebraObject: linearAlgebraObject)
  }
  
  override func map(transform: (Double) -> Double) -> Vector {
    let matrix = super.map(transform)
    return Vector(linearAlgebraObject: matrix.linearAlgebraObject)
  }
}

func +(left: Vector, right: Vector) -> Vector {
  return left.sum(right)
}

func -(left: Vector, right: Vector) -> Vector {
  return left.difference(right)
}

func *(left: Vector, right: Double) -> Vector {
  return left.scale(right)
}

func *(left: Double, right: Vector) -> Vector {
  return right.scale(left)
}

func <*>(left: Vector, right: Vector) -> Vector {
  return left.elementwiseProduct(right)
}

// a^t * b -> t* -> +*
infix operator +* {
  precedence 150
  associativity left
}
func +*(left: Vector, right: Vector) -> Double {
  return left.innerProduct(right)
}

// a * b^t -> *t -> *+
infix operator *+ {
  precedence 150
  associativity left
}
func *+(left: Vector, right: Vector) -> Matrix {
  return left.outerProduct(right)
}

なお、新たに定義している演算子+*と*+は、それぞれベクトルの内積と直積（外積）のための演算子。
内積は左項の転置をとって掛け算、直積（外積）は右項の転置をとって掛け算なので、このような記号にしてみた。（ここでは+が転置の意味）

動作確認

動作確認をするために書いたコードが以下。
このラッパークラスを使うことで、どんな感じに書けるのかが伝わると思う。

//==============================
// LinearAlgebra
//------------------------------
// main.swift
//
// Test code for Linear Algebra
//==============================

import Foundation
import Accelerate

// Generate matrices and vectors with class methods

// matrix1 = [1.0 2.0 3.0]
//           [4.0 5.0 6.0]
let matrix1 = Matrix.fromArray([[1.0, 2.0, 3.0],
                                [4.0, 5.0, 6.0]])

// matrix2 = [2.0 2.0 2.0]
//           [2.0 2.0 2.0]
let matrix2 = Matrix.filledWith(2.0, row: 2, col: 3)

// matrix3 = [1.0 0.0 0.0]
//           [0.0 1.0 0.0]
//           [0.0 0.0 1.0]
let matrix3 = Matrix.identityMatrix(3)

// vector1 = [1.0]
//           [2.0]
//           [3.0]
let vector1 = Vector.fromArray([1.0, 2.0, 3.0])

// matrix4 = [1.0 0.0 0.0]
//           [0.0 2.0 0.0]
//           [0.0 0.0 3.0]
let matrix4 = Matrix.diagonalMatrixFromVector(vector1)

// vector2 = [2.0]
//           [2.0]
//           [2.0]
let vector2 = Vector.filledWith(2.0, size: 3)

// Properties

let objects = [
  "matrix1": matrix1,
  "matrix2": matrix2,
  "matrix3": matrix3,
  "matrix4": matrix4,
  "vector1": vector1,
  "vector2": vector2,
]

for (name, object) in objects {
  print("\(name):\n\(object)")
  print("- row: \(object.row)")
  print("- col: \(object.col)")
  print("- vector?: \(object.isVector)")
  if let vector = object as? Vector {
    print("- size: \(vector.size)")
  }
}

// Index access

for row in 0..<matrix1.row {
  for col in 0..<matrix1.col {
    print("matrix1[\(row), \(col)]: \(matrix1[row, col])")
  }
}

for index in 0..<vector1.size {
  print("vector1[\(index)]: \(vector1[index])")
}

// Generate vectors from matrix

print("matrix1 row0:\n\(matrix1.rowVector(0))")
print("matrix1 row1:\n\(matrix1.rowVector(1))")
print("matrix1 col0:\n\(matrix1.colVector(0))")
print("matrix1 col1:\n\(matrix1.colVector(1))")
print("matrix1 col2:\n\(matrix1.colVector(2))")
print("matrix1 diagonal:\n\(matrix1.diagonalVector())")

// Slice

// slice1 = [1.0 2.0]
//          [4.0 5.0]
let slice1 = matrix1.slice(0, 0, 1, 1, 2, 2)
print("slice1:\n\(slice1)")

// slice2 = [6.0 5.0]
//          [3.0 2.0]
let slice2 = matrix1.slice(1, 2, -1, -1, 2, 2)
print("slice2:\n\(slice2)")

// slice3 = [2.0] (Vector)
//          [5.0]
let slice3 = matrix1.slice(0, 1, 1, 0, 2, 1)
print("slice3:\n\(slice3)")
print("- vector?: \(slice3 is Vector)")

// slice4 = [1.0]
//          [3.0]
let slice4 = vector1.slice(0, 2, 2)
print("slice4:\n\(slice4)")

// slice5 = [3.0]
//          [2.0]
//          [1.0]
let slice5 = vector1.slice(2, -1, 3)
print("slice5:\n\(slice5)")

// Transpose

// matrix1^T = [1.0 4.0]
//             [2.0 5.0]
//             [3.0 6.0]
let matrix1transpose = matrix1.transpose()
print("matrix1 transpose:\n\(matrix1transpose)")

// vector1^T = [1.0 2.0 3.0]
let vector1transpose = vector1.transpose()
print("vector1 transpose:\n\(vector1transpose)")

// Operators

// 2.0 * matrix1
let product1 = 2.0 * matrix1
print("2.0 * matrix1:\n\(product1)")

// matrix1 * 2.0
let product2 = matrix1 * 2.0
print("matrix1 * 2.0:\n\(product2)")

// matrix1 + matrix2
let sum1 = matrix1 + matrix2
print("matrix1 + matrix2:\n\(sum1)")

// matrix1 - matrix2
let difference1 = matrix1 - matrix2
print("matrix1 - matrix2:\n\(difference1)")

// matrix1 <*> matrix2
let product3 = matrix1 <*> matrix2
print("matrix1 <*> matrix2:\n\(product3)")

// 2.0 * vector1
let product4 = 2.0 * vector1
print("2.0 * vector1:\n\(product4)")

// vector1 * 2.0
let product5 = vector1 * 2.0
print("vector1 * 2.0:\n\(product5)")

// vector1 + vector2
let sum2 = vector1 + vector2
print("vector1 + vector2:\n\(sum2)")

// vector1 - vector2
let difference2 = vector1 - vector2
print("vector1 - vector2:\n\(difference2)")

// vector1 <*> vector2
let product6 = vector1 <*> vector2
print("vector1 <*> vector2:\n\(product6)")

// matrix1 * matrix4
//   [1.0 2.0 3.0]   [1.0 0.0 0.0]
// = [4.0 5.0 6.0] * [0.0 2.0 0.0]
//                   [0.0 0.0 3.0]
// = [1.0  4.0  9.0]
//   [4.0 10.0 18.0]
let product7 = matrix1 * matrix4
print("matrix1 * matrix4:\n\(product7)")

// matrix1 * vector1 (Vector)
//   [1.0 2.0 3.0]   [1.0]
// = [4.0 5.0 6.0] * [2.0]
//                   [3.0]
// = [14.0]
//   [32.0]
let product8 = matrix1 * vector1
print("matrix1 * vector1:\n\(product8)")
print("- vector?: \(product8 is Vector)")

// vector1^T * matrix4 (Vector)
//                   [1.0 0.0 0.0]
// = [1.0 2.0 3.0] * [0.0 2.0 0.0]
//                   [0.0 0.0 3.0]
// = [1.0 4.0 9.0]
let product9 = vector1.transpose() * matrix4
print("vector1^T * matrix4:\n\(product9)")
print("- vector?: \(product9 is Vector)")

// vector1 +* vector2
let product10 = vector1 +* vector2
let product11 = (vector1.transpose() * vector2)[0, 0]
print("vector1 +* vector2: \(product10)")
print("vector1^T * vector2: \(product11)")

// vector1 *+ vector2
let product12 = vector1 *+ vector2
let product13 = vector1 * vector2.transpose()
print("vector1 *+ vector2:\n\(product12)")
print("vector1 * vector2^T:\n\(product13)")

// Norm

let norms = [
  "L1": LA_L1_NORM,
  "L2": LA_L2_NORM,
  "LINF": LA_LINF_NORM,
]
for (name, norm) in norms {
  print("vector1 norm (\(name)): \(vector1.norm(Int(norm)))")
  print("- normalized:\n\(vector1.normalizedVector(Int(norm)))")
}

// Map

// exp(matrix1) = [exp(1.0) exp(2.0) exp(3.0)]
//                [exp(4.0) exp(5.0) exp(6.0)]
let mapped1 = matrix1.map(exp)
print("exp(matrix1):\n\(mapped1)")

// exp(vector1) = [exp(1.0)]
//                [exp(2.0)]
//                [exp(3.0)]
let mapped2 = vector1.map(exp)
print("exp(vector1):\n\(mapped2)")

// Error

let badProduct = matrix1 * slice1
print("matrix1 * slice1 is error?: \(badProduct.hasError)")
print("- error code: \(badProduct.errorCode)")

これを実行すると、以下のとおり：

vector1:
[1.0]
[2.0]
[3.0]
- row: 3
- col: 1
- vector?: true
- size: 3
vector2:
[2.0]
[2.0]
[2.0]
- row: 3
- col: 1
- vector?: true
- size: 3
matrix2:
[2.0, 2.0, 2.0]
[2.0, 2.0, 2.0]
- row: 2
- col: 3
- vector?: false
matrix3:
[1.0, 0.0, 0.0]
[0.0, 1.0, 0.0]
[0.0, 0.0, 1.0]
- row: 3
- col: 3
- vector?: false
matrix4:
[1.0, 0.0, 0.0]
[0.0, 2.0, 0.0]
[0.0, 0.0, 3.0]
- row: 3
- col: 3
- vector?: false
matrix1:
[1.0, 2.0, 3.0]
[4.0, 5.0, 6.0]
- row: 2
- col: 3
- vector?: false
matrix1[0, 0]: 1.0
matrix1[0, 1]: 2.0
matrix1[0, 2]: 3.0
matrix1[1, 0]: 4.0
matrix1[1, 1]: 5.0
matrix1[1, 2]: 6.0
vector1[0]: 1.0
vector1[1]: 2.0
vector1[2]: 3.0
matrix1 row0:
[1.0, 2.0, 3.0]
matrix1 row1:
[4.0, 5.0, 6.0]
matrix1 col0:
[1.0]
[4.0]
matrix1 col1:
[2.0]
[5.0]
matrix1 col2:
[3.0]
[6.0]
matrix1 diagonal:
[1.0]
[5.0]
slice1:
[1.0, 2.0]
[4.0, 5.0]
slice2:
[6.0, 5.0]
[3.0, 2.0]
slice3:
[2.0]
[5.0]
- vector?: true
slice4:
[1.0]
[3.0]
slice5:
[3.0]
[2.0]
[1.0]
matrix1 transpose:
[1.0, 4.0]
[2.0, 5.0]
[3.0, 6.0]
vector1 transpose:
[1.0, 2.0, 3.0]
2.0 * matrix1:
[2.0, 4.0, 6.0]
[8.0, 10.0, 12.0]
matrix1 * 2.0:
[2.0, 4.0, 6.0]
[8.0, 10.0, 12.0]
matrix1 + matrix2:
[3.0, 4.0, 5.0]
[6.0, 7.0, 8.0]
matrix1 - matrix2:
[-1.0, 0.0, 1.0]
[2.0, 3.0, 4.0]
matrix1 <*> matrix2:
[2.0, 4.0, 6.0]
[8.0, 10.0, 12.0]
2.0 * vector1:
[2.0]
[4.0]
[6.0]
vector1 * 2.0:
[2.0]
[4.0]
[6.0]
vector1 + vector2:
[3.0]
[4.0]
[5.0]
vector1 - vector2:
[-1.0]
[0.0]
[1.0]
vector1 <*> vector2:
[2.0]
[4.0]
[6.0]
matrix1 * matrix4:
[1.0, 4.0, 9.0]
[4.0, 10.0, 18.0]
matrix1 * vector1:
[14.0]
[32.0]
- vector?: true
vector1^T * matrix4:
[1.0, 4.0, 9.0]
- vector?: true
vector1 +* vector2: 12.0
vector1^T * vector2: 12.0
vector1 *+ vector2:
[2.0, 2.0, 2.0]
[4.0, 4.0, 4.0]
[6.0, 6.0, 6.0]
vector1 * vector2^T:
[2.0, 2.0, 2.0]
[4.0, 4.0, 4.0]
[6.0, 6.0, 6.0]
vector1 norm (L1): 6.0
- normalized:
[0.166666666666667]
[0.333333333333333]
[0.5]
vector1 norm (L2): 3.74165738677394
- normalized:
[0.267261241912424]
[0.534522483824849]
[0.801783725737273]
vector1 norm (LINF): 3.0
- normalized:
[0.333333333333333]
[0.666666666666667]
[1.0]
exp(matrix1):
[2.71828182845905, 7.38905609893065, 20.0855369231877]
[54.5981500331442, 148.413159102577, 403.428793492735]
exp(vector1):
[2.71828182845905]
[7.38905609893065]
[20.0855369231877]
matrix1 * slice1 is error?: true
- error code: -1002

問題なく動作しているのが分かると思う。

今日はここまで！

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

昨日はベクトルとSplatの説明をした。

今日はこれらを使った演算について。

転置

まずは転置。
次の関数が用意されている：

// 転置した行列を生成して返す。
la_object_t
la_transpose(
    la_object_t matrix);  // 行列

スカラー倍

次はスカラー倍。

// スカラー倍した行列を生成して返す。
la_object_t
la_scale_with_float(
    la_object_t matrix,  // 行列
    float scalar);       // 何倍するか
la_object_t
la_scale_with_double(
    la_object_t matrix,  // 行列
    double scalar);      // 何倍するか

2つの行列の和と差

2つの行列の和と差の計算。

// 2つの行列の和を生成して返す。
// （サイズが同じでないとダメ）
la_object_t
la_sum(
    la_object_t obj_left,    // 左項
    la_object_t obj_right);  // 右項

// 2つの行列の差を生成して返す。
// （サイズが同じでないとダメ）
la_object_t
la_difference(
    la_object_t obj_left,    // 左項
    la_object_t obj_right);  // 右項

なお、サイズが違うなどで計算が正しくできなかった場合には返されるオブジェクトの状態がエラーになるので、la_status()で状態をチェックする。
（以降の他の関数も同様）

2つのベクトルの内積、外積（直積）

2つのベクトル の内積 、および外積（直積） は、次の関数で計算できる：

// ベクトルの内積。
// （サイズが同じでないとダメ）
// 結果がスカラーでなく1x1の行列になることに注意。
la_object_t
la_inner_product(
    la_object_t vector_left,    // 左項
    la_object_t vector_right);  // 右項

// ベクトルの外積（直積）。
la_object_t
la_outer_product(
    la_object_t vector_left,    // 左項
    la_object_t vector_right);  // 右項

なお、コメントにも書いた通り、内積はスカラーではなくて1x1の行列になるので、注意が必要。
（確かに、1行の行列と1列の行列の積が内積なので、1x1行列ではあるんだけど・・・）

あと、ベクトルには列ベクトル・行ベクトルのいずれを指定してもいい。

2つの行列の積、アダマール積

2つの行列 の積 、およびアダマール積（要素ごとの積） は、次の関数で計算できる：

// 行列の積。
// （左項の列数と右項の行数が同じでないとダメ）
la_object_t
la_matrix_product(
    la_object_t matrix_left,    // 左項
    la_object_t matrix_right);  // 右項

// 行列のアダマール積（要素ごとの積）
// （サイズが同じでないとダメ）
// ただし、ベクトルについては、列ベクトルか行ベクトルかは問われない。
// （列ベクトルと行ベクトルが指定された場合、列ベクトルになる）
la_object_t
la_elementwise_product(
    la_object_t obj_left,    // 左項
    la_object_t obj_right);  // 右項

ベクトルのノルムと正規化

ベクトルのノルムは、次の関数で計算できる：

// ベクトルのノルム
float
la_norm_as_float(
    la_object_t vector,      // ベクトル
    la_norm_t vector_norm);  // ノルムの種類
double
la_norm_as_double(
    la_object_t vector,      // ベクトル
    la_norm_t vector_norm);  // ノルムの種類

la_norm_tはノルムの種類を指定するためのもので、次のように定義されている：

typedef unsigned long la_norm_t;
#define LA_L1_NORM 1    // L1ノルム
#define LA_L2_NORM 2    // L2ノルム
#define LA_LINF_NORM 3  // 無限大ノルム（最大値ノルム）

また、ベクトルの正規化を行うための関数も用意されている：

// 指定されたノルムで正規化したベクトルを生成する。
la_object_t
la_normalized_vector(
    la_object_t vector,      // ベクトル
    la_norm_t vector_norm);  // ノルムの種類

連立一次方程式の解を求める

連立一次方程式（を複数並べたもの） を満たす行列 を求める関数が用意されている：

// AX = Bの解Xを返す。
la_object_t
la_solve(
    la_object_t matrix_system,  // Aにあたる行列
    la_object_t obj_rhs);       // Bにあたる行列

---

用意されているインタフェースはこれで終わりなんだけど、C言語の関数の形なので、Swiftで使うにはちょっと野暮な感じ。
なので、明日はラッパークラスを作ってみる。

今日はここまで！

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

昨日はAccelerateフレームワークの概要と、行列の生成と内容の取得について説明した。

今日はベクトルとSplatについて。

ベクトルの生成

ベクトルといっても、実際には1行もしくは1列の行列。
ただ、扱いやすいように、いろいろマクロや関数が用意されている。

配列から生成する

配列からベクトルを生成できるようにするために、次のマクロが用意されている：

// Float/Doubleのバッファ（配列）から列ベクトルを生成する。
// （バッファの内容はコピーされる）
#define la_vector_from_float_buffer(buffer, vector_length, buffer_stride, attributes) \
    la_matrix_from_float_buffer(buffer, vector_length, 1, buffer_stride, LA_NO_HINT, attributes)
#define la_vector_from_double_buffer(buffer, vector_length, buffer_stride, attributes) \
    la_matrix_from_double_buffer(buffer, vector_length, 1, buffer_stride, LA_NO_HINT, attributes)

ただ、Swiftだと使えないみたい。

もし関数で使いたい場合、自分で次のように定義するといいかもしれない：
（もっといい方法は後述）

func la_vector_from_float_buffer(buffer: [Float],
                                 _ vector_length: la_count_t,
                                 _ buffer_stride: la_count_t,
                                 _ attributes: la_attribute_t) -> la_object_t {
  return la_matrix_from_float_buffer(buffer,
                                     vector_length, la_count_t(1),
                                     buffer_stride,
                                     la_hint_t(LA_NO_HINT),
                                     attributes)
}

func la_vector_from_double_buffer(buffer: [Double],
                                  _ vector_length: la_count_t,
                                  _ buffer_stride: la_count_t,
                                  _ attributes: la_attribute_t) -> la_object_t {
  return la_matrix_from_double_buffer(buffer,
                                      vector_length, la_count_t(1),
                                      buffer_stride,
                                      la_hint_t(LA_NO_HINT),
                                      attributes)
}

部分的なベクトルを生成する

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

// 部分的なベクトルを生成して返す。
/a_object_t
la_vector_slice(
    la_object_t vector,        // 元のベクトル
    la_index_t vector_first,   // ピックアップする最初のインデックス
    la_index_t vector_stride,  // ピックアップしていくときにズラす幅
    la_count_t slice_length);  // ピックアップする要素数

逆順のベクトルを生成する

すでにあるベクトルの順序を逆順にしたベクトルを生成するマクロが用意されている：

// 逆順にしたベクトルを生成して返す。
#define la_vector_reverse(vector) \
    la_vector_slice(vector, la_vector_length(vector)-1, -1, la_vector_length(vector))

ただ、マクロなので、やはりSwiftでは使えない。
必要なら、自分で関数を定義する必要がある。

ベクトルのサイズ

ベクトルのサイズは次の関数で取得できる：

// ベクトルのサイズを返す。
// （引数には行ベクトルも列ベクトルも指定できる）
la_count_t
la_vector_length(la_object_t vector);

ベクトルの内容の取得

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

もちろん、配列に行列の内容を書き込む関数を使うのでもいいんだけど、次の関数も用意されている：

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

buffer_strideは行列版のbuffer_row_strideとほぼ同じで、ベクトルが列ベクトルの場合、実際同じ。
ただ、この関数では、行ベクトルの場合もうまく働くようになっている。

Splatの生成

Splatはすべての要素の値が等しい行列（もしくはベクトル）。
ただ、サイズは決まっていなくて、演算のもう一方のオブジェクトのサイズに合わせて自動的に調整される。
例えば、2x3の行列と足し算をするときには2x3の行列として振舞うし、5x1のベクトルと内積を取るときには5x1のベクトルとして振舞う。

Splatを生成するための関数として、以下のものが用意されている：

// 指定されたFloat/Doubleの値でSplatを生成する。
la_object_t
la_splat_from_float(
    float scalar_value,          // Floatの値
    la_attribute_t attributes);  // 属性
la_object_t
la_splat_from_double(
    double scalar_value,         // Doubleの値
    la_attribute_t attributes);  // 属性

また、ベクトルや行列の要素の値からSplatを生成することも出来る：

// 指定されたベクトルのインデックスにある値でSplatを生成する。
la_object_t
la_splat_from_vector_element(
    la_object_t vector,        // ベクトル
    la_index_t vector_index);  // ベクトルのインデックス

// 指定された行列のインデックスにある値でSplatを生成する。
la_object_t
la_splat_from_matrix_element(
    la_object_t matrix,      // 行列
    la_index_t matrix_row,   // 行列の行
    la_index_t matrix_col);  // 行列の列

Splatからベクトル、行列を生成する

Splatからサイズを指定してベクトル、行列を生成することも出来る：

// 指定されたサイズの列ベクトルを生成する。
la_object_t
la_vector_from_splat(
    la_object_t splat,          // Splat
    la_count_t vector_length);  // 生成するベクトルのサイズ

// 指定されたサイズの行列を生成する。
la_object_t la_matrix_from_splat(
    la_object_t splat,        // Splat
    la_count_t matrix_rows,   // 生成する行列の行数
    la_count_t matrix_cols);  // 生成する行列の列数

実際の演算では、いろんな演算の一方のサイズがいろいろ変わるということはあまり考えられないので、これらの関数を使ってSplatを（固定サイズの）行列やベクトルに変換してしまった方がいいと思う。

今日はここまで！