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

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

今日はここまで！