New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Matrix eigendecomposition #95

Merged
merged 10 commits into from Feb 9, 2019

Conversation

Projects
None yet
5 participants
@loufranco
Copy link
Contributor

loufranco commented Feb 8, 2019

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)

  1. I got the older branch updated and fixed conflicts
  2. I changed the interface to support complex vector and matrix returns
  3. I added tests based on Intel LAPACK and numpy docs

But,

  1. My return is a 3-tuple, which generates a lint warning
  2. I am using (Double, Double) for complex numbers. I don't see a convention in Surge for this, but if there is one, let me know
  3. I can't use Matrix with that, so I use [[(Double, Double)]]

Ben Robbins and others added some commits Apr 21, 2015

Ben Robbins Ben Robbins
Added Two function that calculate the eigen-values and right eigen ve…
…ctors for a given Matrix. Further changes may be appropriate
[Issue #14] Getting PR to compile
Refactoring and reformatting
Fixed bugs in eigendecomposition and added tests
- 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
@alejandro-isaza
Copy link
Collaborator

alejandro-isaza left a comment

Thanks! Just a few style comments.

// 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 {

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

I'm not a fan of one-letter variable names, let use better names.

// .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)]]) {

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

Use three slashes for doc comments /// also use swift doc comment style. See https://nshipster.com/swift-documentation/

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

Let's define a struct for the return type.

// .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)]]) {

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

Same comments as above.

// dgeev_ needs column-major matrices
var mat: [__CLPK_doublereal] = transpose(x).grid

var N1 = __CLPK_integer(x.rows)

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

variable names should be lowercase, also let's use better names.

@loufranco

This comment has been minimized.

Copy link
Contributor Author

loufranco commented Feb 9, 2019

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

  1. Suppress it
  2. Make a struct to hold the return
  3. Only return the right eigenvector (this is what numpy does), and use a 2-tuple -- this would make the function have a very analogous signature to numpy. I don't know how important having the left eigenvector is.
Implement PR Feedback
- better variable names
- use a struct for the decomposition result
- use standard Swift doc comments
@loufranco

This comment has been minimized.

Copy link
Contributor Author

loufranco commented Feb 9, 2019

@alejandro-isaza I implemented your suggestions.

@alejandro-isaza
Copy link
Collaborator

alejandro-isaza left a comment

Just a couple more things.

/// - 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> {

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

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.

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")

This comment has been minimized.

@alejandro-isaza

alejandro-isaza Feb 9, 2019

Collaborator

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.

This comment has been minimized.

@loufranco

loufranco Feb 9, 2019

Author Contributor

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.

Respond to PR Feedback
- rename func to eigenDecompose
- fix some doc comments
- throw instead of asserting
- added a test for non-square matrices
@loufranco

This comment has been minimized.

Copy link
Contributor Author

loufranco commented Feb 9, 2019

@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.

@alejandro-isaza alejandro-isaza merged commit a6dd5f6 into mattt:master Feb 9, 2019

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@alejandro-isaza

This comment has been minimized.

Copy link
Collaborator

alejandro-isaza commented Feb 9, 2019

Thanks again!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment