Matrix eigendecomposition (#95)

Added Two function that calculate the eigen-values and right eigen vectors for a given Matrix.
Jounce · Feb 9, 2019 · a6dd5f6 · a6dd5f6
1 parent 3a99c1d
commit a6dd5f6
Show file tree

Hide file tree

Showing 2 changed files with 251 additions and 0 deletions.
diff --git a/Sources/Surge/Matrix.swift b/Sources/Surge/Matrix.swift
@@ -117,6 +117,42 @@ public struct Matrix<Scalar> where Scalar: FloatingPoint, Scalar: ExpressibleByF
     }
 }
 
+/// Holds the result of eigendecomposition. The (Scalar, Scalar) used
+/// in the property types represents a complex number with (real, imaginary) parts.
+public struct MatrixEigenDecompositionResult<Scalar> where Scalar: FloatingPoint, Scalar: ExpressibleByFloatLiteral {
+    let eigenValues: [(Scalar, Scalar)]
+    let leftEigenVectors: [[(Scalar, Scalar)]]
+    let rightEigenVectors: [[(Scalar, Scalar)]]
+
+    public init(eigenValues: [(Scalar, Scalar)], leftEigenVectors: [[(Scalar, Scalar)]], rightEigenVectors: [[(Scalar, Scalar)]]) {
+        self.eigenValues = eigenValues
+        self.leftEigenVectors = leftEigenVectors
+        self.rightEigenVectors = rightEigenVectors
+    }
+
+    public init(rowCount: Int, eigenValueRealParts: [Scalar], eigenValueImaginaryParts: [Scalar], leftEigenVectorWork: [Scalar], rightEigenVectorWork: [Scalar]) {
+        // The eigenvalues are an array of (real, imaginary) results from dgeev
+        self.eigenValues = Array(zip(eigenValueRealParts, eigenValueImaginaryParts))
+
+        // Build the left and right eigenvectors
+        let emptyVector = [(Scalar, Scalar)](repeating: (0, 0), count: rowCount)
+        var leftEigenVectors = [[(Scalar, Scalar)]](repeating: emptyVector, count: rowCount)
+        buildEigenVector(eigenValueImaginaryParts: eigenValueImaginaryParts, eigenVectorWork: leftEigenVectorWork, result: &leftEigenVectors)
+
+        var rightEigenVectors = [[(Scalar, Scalar)]](repeating: emptyVector, count: rowCount)
+        buildEigenVector(eigenValueImaginaryParts: eigenValueImaginaryParts, eigenVectorWork: rightEigenVectorWork, result: &rightEigenVectors)
+
+        self.leftEigenVectors = leftEigenVectors
+        self.rightEigenVectors = rightEigenVectors
+    }
+}
+
+/// Errors thrown when a matrix cannot be decomposed.
+public enum EigenDecompositionError: Error {
+    case matrixNotSquare
+    case matrixNotDecomposable
+}
+
 // MARK: - Printable
 
 extension Matrix: CustomStringConvertible {
@@ -469,6 +505,97 @@ public func det(_ x: Matrix<Double>) -> Double? {
     return det
 }
 
+// Convert the result of dgeev into an array of complex numbers
+// See Intel's documentation on column-major results for sample code that this
+// is based on:
+// https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev.htm
+private func buildEigenVector<Scalar: FloatingPoint & ExpressibleByFloatLiteral>(eigenValueImaginaryParts: [Scalar], eigenVectorWork: [Scalar], result: inout [[(Scalar, Scalar)]]) {
+    // row and col count are the same because result must be square.
+    let rowColCount = result.count
+
+    for row in 0..<rowColCount {
+        var col = 0
+        while col < rowColCount {
+            if eigenValueImaginaryParts[col] == 0.0 {
+                // v is column-major
+                result[row][col] = (eigenVectorWork[row+rowColCount*col], 0.0)
+                col += 1
+            } else {
+                // v is column-major
+                result[row][col] = (eigenVectorWork[row+col*rowColCount], eigenVectorWork[row+rowColCount*(col+1)])
+                result[row][col+1] = (eigenVectorWork[row+col*rowColCount], -eigenVectorWork[row+rowColCount*(col+1)])
+                col += 2
+            }
+        }
+    }
+}
+
+/// Decomposes a square matrix into its eigenvalues and left and right eigenvectors.
+/// The decomposition may result in complex numbers, represented by (Float, Float), which
+///   are the (real, imaginary) parts of the complex number.
+/// - Parameters:
+///   - x: a square matrix
+/// - Returns: a struct with the eigen values and left and right eigen vectors using (Float, Float)
+///   to represent a complex number.
+public func eigenDecompose(_ x: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> {
+    var input = Matrix<Double>(rows: x.rows, columns: x.columns, repeatedValue: 0.0)
+    input.grid = x.grid.map { Double($0) }
+    let decomposition = try eigenDecompose(input)
+
+    return MatrixEigenDecompositionResult<Float>(
+        eigenValues: decomposition.eigenValues.map { (Float($0.0), Float($0.1)) },
+        leftEigenVectors: decomposition.leftEigenVectors.map { $0.map { (Float($0.0), Float($0.1)) } },
+        rightEigenVectors: decomposition.rightEigenVectors.map { $0.map { (Float($0.0), Float($0.1)) } }
+    )
+}
+
+/// Decomposes a square matrix into its eigenvalues and left and right eigenvectors.
+/// The decomposition may result in complex numbers, represented by (Double, Double), which
+///   are the (real, imaginary) parts of the complex number.
+/// - Parameters:
+///   - x: a square matrix
+/// - Returns: a struct with the eigen values and left and right eigen vectors using (Double, Double)
+///   to represent a complex number.
+public func eigenDecompose(_ x: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> {
+    guard x.rows == x.columns else {
+        throw EigenDecompositionError.matrixNotSquare
+    }
+
+    // dgeev_ needs column-major matrices, so transpose x.
+    var matrixGrid: [__CLPK_doublereal] = transpose(x).grid
+    var matrixRowCount = __CLPK_integer(x.rows)
+    let matrixColCount = matrixRowCount
+    var eigenValueCount = matrixRowCount
+    var leftEigenVectorCount = matrixRowCount
+    var rightEigenVectorCount = matrixRowCount
+
+    var workspaceQuery: Double = 0.0
+    var workspaceSize = __CLPK_integer(-1)
+    var error: __CLPK_integer = 0
+
+    var eigenValueRealParts = [Double](repeating: 0, count: Int(eigenValueCount))
+    var eigenValueImaginaryParts = [Double](repeating: 0, count: Int(eigenValueCount))
+    var leftEigenVectorWork = [Double](repeating: 0, count: Int(leftEigenVectorCount * matrixColCount))
+    var rightEigenVectorWork = [Double](repeating: 0, count: Int(rightEigenVectorCount * matrixColCount))
+
+    let decompositionJobV = UnsafeMutablePointer<Int8>(mutating: ("V" as NSString).utf8String)
+    // Call dgeev to find out how much workspace to allocate
+    dgeev_(decompositionJobV, decompositionJobV, &matrixRowCount, &matrixGrid, &eigenValueCount, &eigenValueRealParts, &eigenValueImaginaryParts, &leftEigenVectorWork, &leftEigenVectorCount, &rightEigenVectorWork, &rightEigenVectorCount, &workspaceQuery, &workspaceSize, &error)
+    if error != 0 {
+        throw EigenDecompositionError.matrixNotDecomposable
+    }
+
+    // Allocate the workspace and call dgeev again to do the actual decomposition
+    var workspace = [Double](repeating: 0.0, count: Int(workspaceQuery))
+    workspaceSize = __CLPK_integer(workspaceQuery)
+    dgeev_(decompositionJobV, decompositionJobV, &matrixRowCount, &matrixGrid, &eigenValueCount, &eigenValueRealParts, &eigenValueImaginaryParts, &leftEigenVectorWork, &leftEigenVectorCount, &rightEigenVectorWork, &rightEigenVectorCount, &workspace, &workspaceSize, &error)
+    if error != 0 {
+        throw EigenDecompositionError.matrixNotDecomposable
+    }
+
+    return MatrixEigenDecompositionResult<Double>(rowCount: x.rows, eigenValueRealParts: eigenValueRealParts, eigenValueImaginaryParts: eigenValueImaginaryParts, leftEigenVectorWork: leftEigenVectorWork, rightEigenVectorWork: rightEigenVectorWork)
+}
+
 // MARK: - Operators
 
 public func + (lhs: Matrix<Float>, rhs: Matrix<Float>) -> Matrix<Float> {

diff --git a/Tests/SurgeTests/MatrixTests.swift b/Tests/SurgeTests/MatrixTests.swift
@@ -93,4 +93,128 @@ class MatrixTests: XCTestCase {
         XCTAssertEqual(result.rows, 1)
         XCTAssertEqual(result.columns, 0)
     }
+
+    func complexVectorMatches<T: FloatingPoint>(_ a: [(T, T)], _ b: [(T, T)], accuracy: T) -> Bool {
+        guard a.count == b.count else {
+            return false
+        }
+        return (zip(a, b).first { a, e -> Bool in
+            !(abs(a.0 - e.0) < accuracy && abs(a.1 - e.1) < accuracy)
+        }) == nil
+    }
+
+    func testEigenDecompositionTrivialGeneric<T: FloatingPoint & ExpressibleByFloatLiteral>(defaultAccuracy: T, eigendecompostionFn: ((Matrix<T>) throws -> MatrixEigenDecompositionResult<T>)) throws {
+        let matrix = Matrix<T>([
+            [1, 0, 0],
+            [0, 2, 0],
+            [0, 0, 3],
+        ])
+        let ed = try eigendecompostionFn(matrix)
+
+        let expectedEigenValues: [(T, T)] = [(1, 0), (2, 0), (3, 0)]
+        XCTAssertTrue(complexVectorMatches(ed.eigenValues, expectedEigenValues, accuracy: defaultAccuracy))
+
+        let expectedEigenVector: [(T, T)] = [(1, 0), (0, 0), (0, 0), (0, 0), (1, 0), (0, 0), (0, 0), (0, 0), (1, 0)]
+        let flatLeft = ed.leftEigenVectors.flatMap { $0 }
+        XCTAssertTrue(complexVectorMatches(flatLeft, expectedEigenVector, accuracy: defaultAccuracy))
+
+        let flatRight = ed.rightEigenVectors.flatMap { $0 }
+        XCTAssertTrue(complexVectorMatches(flatRight, expectedEigenVector, accuracy: defaultAccuracy))
+    }
+
+    func testEigenDecompositionTrivial() throws {
+        try testEigenDecompositionTrivialGeneric(defaultAccuracy: doubleAccuracy) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
+            try eigenDecompose(m)
+        }
+        try testEigenDecompositionTrivialGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
+            try eigenDecompose(m)
+        }
+    }
+
+    // Example from https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html
+    func testEigenDecompositionComplexResultsNumpyExampleGeneric<T: FloatingPoint & ExpressibleByFloatLiteral>(defaultAccuracy: T, eigendecompostionFn: ((Matrix<T>) throws -> MatrixEigenDecompositionResult<T>)) throws {
+        let matrix = Matrix<T>([
+            [1, -1],
+            [1, 1],
+        ])
+        let ed = try eigendecompostionFn(matrix)
+
+        let expectedEigenValues: [(T, T)] = [(1, 1), (1, -1)]
+        XCTAssertTrue(complexVectorMatches(ed.eigenValues, expectedEigenValues, accuracy: defaultAccuracy))
+
+        let expectedEigenVector: [(T, T)] = [(0.707_106_78, 0), (0.707_106_78, 0), (0, -0.707_106_78), (0, 0.707_106_78)]
+
+        // Numpy only gives the right eigenvector
+        let flatRight = ed.rightEigenVectors.flatMap { $0 }
+        XCTAssertTrue(complexVectorMatches(flatRight, expectedEigenVector, accuracy: 0.000_001))
+    }
+
+    func testEigenDecompositionComplexResultsNumpyExample() throws {
+        try testEigenDecompositionComplexResultsNumpyExampleGeneric(defaultAccuracy: doubleAccuracy) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
+            try eigenDecompose(m)
+        }
+        try testEigenDecompositionComplexResultsNumpyExampleGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
+            try eigenDecompose(m)
+        }
+    }
+
+    // Example from Intel's DGEEV documentation
+    // https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev.htm
+    func testEigenDecompositionComplexResultsDGEEVExampleGeneric<T: FloatingPoint & ExpressibleByFloatLiteral>(defaultAccuracy: T, eigendecompostionFn: ((Matrix<T>) throws -> MatrixEigenDecompositionResult<T>)) throws {
+        let matrix = Matrix<T>(rows: 5, columns: 5, grid: [
+            -1.01, 0.86, -4.60, 3.31, -4.81,
+            3.98, 0.53, -7.04, 5.29, 3.55,
+            3.30, 8.26, -3.89, 8.20, -1.51,
+            4.43, 4.96, -7.66, -7.33, 6.18,
+            7.31, -6.43, -6.16, 2.47, 5.58
+        ])
+        let ed = try eigendecompostionFn(matrix)
+
+        let expectedEigenValues: [(T, T)] = [(2.86, 10.76), (2.86, -10.76), (-0.69, 4.70), (-0.69, -4.70), (-10.46, 0)]
+        XCTAssertTrue(complexVectorMatches(ed.eigenValues, expectedEigenValues, accuracy: 0.02))
+
+        let expectedLeftEigenVectors: [[(T, T)]] = [
+            [(0.04, 0.29), (0.04, -0.29), (-0.13, -0.33), (-0.13, 0.33), (0.04, 0)],
+            [(0.62, 0.00), (0.62, 0.00), (0.69, 0.00), (0.69, 0.00), (0.56, 0)],
+            [(-0.04, -0.58), (-0.04, 0.58), (-0.39, -0.07), (-0.39, 0.07), (-0.13, 0)],
+            [(0.28, 0.01), (0.28, -0.01), (-0.02, -0.19), (-0.02, 0.19), (-0.80, 0)],
+            [(-0.04, 0.34), (-0.04, -0.34), (-0.40, 0.22), (-0.40, -0.22), (0.18, 0)],
+        ]
+        XCTAssertEqual(ed.leftEigenVectors.count, expectedLeftEigenVectors.count)
+        for i in 0..<ed.leftEigenVectors.count {
+            XCTAssertTrue(complexVectorMatches(ed.leftEigenVectors[i], expectedLeftEigenVectors[i], accuracy: 0.02))
+        }
+
+        let expectedRightEigenVectors: [[(T, T)]] = [
+            [(0.11, 0.17), (0.11, -0.17), (0.73, 0.00), (0.73, 0.00), (0.46, 0.0)],
+            [(0.41, -0.26), (0.41, 0.26), (-0.03, -0.02), (-0.03, 0.02), (0.34, 0.0)],
+            [(0.10, -0.51), (0.10, 0.51), (0.19, -0.29), (0.19, 0.29), (0.31, 0.0)],
+            [(0.40, -0.09), (0.40, 0.09), (-0.08, -0.08), (-0.08, 0.08), (-0.74, 0.0)],
+            [(0.54, 0.00), (0.54, 0.00), (-0.29, -0.49), (-0.29, 0.49), (0.16, 0.0)],
+        ]
+        XCTAssertEqual(ed.rightEigenVectors.count, expectedRightEigenVectors.count)
+        for i in 0..<ed.rightEigenVectors.count {
+            XCTAssertTrue(complexVectorMatches(ed.rightEigenVectors[i], expectedRightEigenVectors[i], accuracy: 0.02))
+        }
+    }
+
+    func testEigenDecompositionComplexResultsDGEEVExample() throws {
+        try testEigenDecompositionComplexResultsDGEEVExampleGeneric(defaultAccuracy: doubleAccuracy) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
+            try eigenDecompose(m)
+        }
+        try testEigenDecompositionComplexResultsDGEEVExampleGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
+            try eigenDecompose(m)
+        }
+    }
+
+    func testEigenDecomposeThrowsWhenNotSquare() throws {
+        let matrix = Matrix<Double>([
+            [1, 0, 0],
+            [0, 2, 0],
+        ])
+
+        XCTAssertThrowsError(try eigenDecompose(matrix), "not square") { e in
+            XCTAssertEqual(e as? EigenDecompositionError, EigenDecompositionError.matrixNotSquare)
+        }
+    }
 }