From a6dd5f622ed5f36c74ead2dc54329a2d8f39f855 Mon Sep 17 00:00:00 2001 From: Lou Franco Date: Sat, 9 Feb 2019 10:18:07 -0500 Subject: [PATCH] Matrix eigendecomposition (#95) Added Two function that calculate the eigen-values and right eigen vectors for a given Matrix. --- Sources/Surge/Matrix.swift | 127 +++++++++++++++++++++++++++++ Tests/SurgeTests/MatrixTests.swift | 124 ++++++++++++++++++++++++++++ 2 files changed, 251 insertions(+) diff --git a/Sources/Surge/Matrix.swift b/Sources/Surge/Matrix.swift index d7fcbb16..1c13beea 100644 --- a/Sources/Surge/Matrix.swift +++ b/Sources/Surge/Matrix.swift @@ -117,6 +117,42 @@ public struct Matrix 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 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? { 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(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..) throws -> MatrixEigenDecompositionResult { + var input = Matrix(rows: x.rows, columns: x.columns, repeatedValue: 0.0) + input.grid = x.grid.map { Double($0) } + let decomposition = try eigenDecompose(input) + + return MatrixEigenDecompositionResult( + 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) throws -> MatrixEigenDecompositionResult { + 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(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(rowCount: x.rows, eigenValueRealParts: eigenValueRealParts, eigenValueImaginaryParts: eigenValueImaginaryParts, leftEigenVectorWork: leftEigenVectorWork, rightEigenVectorWork: rightEigenVectorWork) +} + // MARK: - Operators public func + (lhs: Matrix, rhs: Matrix) -> Matrix { diff --git a/Tests/SurgeTests/MatrixTests.swift b/Tests/SurgeTests/MatrixTests.swift index 1b357b6b..48e009b5 100644 --- 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(_ 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(defaultAccuracy: T, eigendecompostionFn: ((Matrix) throws -> MatrixEigenDecompositionResult)) throws { + let matrix = Matrix([ + [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) throws -> MatrixEigenDecompositionResult in + try eigenDecompose(m) + } + try testEigenDecompositionTrivialGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix) throws -> MatrixEigenDecompositionResult in + try eigenDecompose(m) + } + } + + // Example from https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html + func testEigenDecompositionComplexResultsNumpyExampleGeneric(defaultAccuracy: T, eigendecompostionFn: ((Matrix) throws -> MatrixEigenDecompositionResult)) throws { + let matrix = Matrix([ + [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) throws -> MatrixEigenDecompositionResult in + try eigenDecompose(m) + } + try testEigenDecompositionComplexResultsNumpyExampleGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix) throws -> MatrixEigenDecompositionResult 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(defaultAccuracy: T, eigendecompostionFn: ((Matrix) throws -> MatrixEigenDecompositionResult)) throws { + let matrix = Matrix(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..) throws -> MatrixEigenDecompositionResult in + try eigenDecompose(m) + } + try testEigenDecompositionComplexResultsDGEEVExampleGeneric(defaultAccuracy: floatAccuracy) { (m: Matrix) throws -> MatrixEigenDecompositionResult in + try eigenDecompose(m) + } + } + + func testEigenDecomposeThrowsWhenNotSquare() throws { + let matrix = Matrix([ + [1, 0, 0], + [0, 2, 0], + ]) + + XCTAssertThrowsError(try eigenDecompose(matrix), "not square") { e in + XCTAssertEqual(e as? EigenDecompositionError, EigenDecompositionError.matrixNotSquare) + } + } }