Matrix eigendecomposition #95
Conversation
…ctors for a given Matrix. Further changes may be appropriate
Refactoring and reformatting
- dgeev wants matrices in column major form - dgeev returns 3 things, so add all to return - decomposition results in complex results, which we represent as a 2-tuple
Sources/Surge/Matrix.swift
Outdated
| // https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgeev.htm | ||
| private func buildEigenVector(wi: [Double], v: [Double], result: inout [[(Double, Double)]]) { | ||
| let n = wi.count | ||
| for r in 0..<n { |
There was a problem hiding this comment.
I'm not a fan of one-letter variable names, let use better names.
Sources/Surge/Matrix.swift
Outdated
| // .0: The eigenvalues represented as an array of complex numbers | ||
| // .1: The left eigenvectors represented as a 2-dimensional array of complex numbers | ||
| // .2: The right eigenvectors represented as a 2-dimensional array of complex numbers | ||
| public func eigendecompostion(_ x: Matrix<Float>) -> ([(Float, Float)], [[(Float, Float)]], [[(Float, Float)]]) { |
There was a problem hiding this comment.
Use three slashes for doc comments /// also use swift doc comment style. See https://nshipster.com/swift-documentation/
There was a problem hiding this comment.
Let's define a struct for the return type.
Sources/Surge/Matrix.swift
Outdated
| // .0: The eigenvalues represented as an array of complex numbers | ||
| // .1: The left eigenvectors represented as a 2-dimensional array of complex numbers | ||
| // .2: The right eigenvectors represented as a 2-dimensional array of complex numbers | ||
| public func eigendecompostion(_ x: Matrix<Double>) -> ([(Double, Double)], [[(Double, Double)]], [[(Double, Double)]]) { |
There was a problem hiding this comment.
Same comments as above.
Sources/Surge/Matrix.swift
Outdated
| // dgeev_ needs column-major matrices | ||
| var mat: [__CLPK_doublereal] = transpose(x).grid | ||
|
|
||
| var N1 = __CLPK_integer(x.rows) |
There was a problem hiding this comment.
variable names should be lowercase, also let's use better names.
|
The build failure is because of strict lint issues with the 3-tuples I used I'd like some guidance on the lint warning about using a tuple of 3 items. I could
|
- better variable names - use a struct for the decomposition result - use standard Swift doc comments
|
@alejandro-isaza I implemented your suggestions. |
Sources/Surge/Matrix.swift
Outdated
| /// - 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 eigendecompostion(_ x: Matrix<Double>) -> MatrixEigenDecompositionResult<Double> { |
There was a problem hiding this comment.
This is a lot better. Let's make this function a verb. What about decompose? or eigenDecompose? or decomposeSpectrally? Also missing triple slash in the first three lines.
Sources/Surge/Matrix.swift
Outdated
| 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) | ||
| assert(error == 0, "Matrix not eigen decomposable") |
There was a problem hiding this comment.
Asserts are bad here because they will be ignored in release build which will probably make the below code crash. Let's make a enum EigenDecompositionError and throw that.
There was a problem hiding this comment.
I will fix -- I was trying to follow the conventions in the rest of the file (and Surge seems generally not to throw), but I agree that it's better.
- rename func to eigenDecompose - fix some doc comments - throw instead of asserting - added a test for non-square matrices
|
@alejandro-isaza ok -- throw is in (with a test for non-square). As far as I can tell, there is no real matrix that isn't decomposable to complex results, but since DGEEV could return an error, I throw -- there is no test I could write to test it, though. |
|
Thanks again! |
This PR tries to bring the almost-done eigendecomposition branch into Surge
I expect there to be some comments (it is not lint warning free right now)
But,
(Double, Double)for complex numbers. I don't see a convention in Surge for this, but if there is one, let me know[[(Double, Double)]]