Complex Support

[spinicist]

Hello,

I want an Eigen compatible way to calculate the first N singular values/vectors of a complex float matrix. This project looks great, but I think at the moment it only supports real-valued matrices? What would be needed to support complex?

PS - Any chance you can add a license file to the repo?

[LTLA]

What would be needed to support complex?

Some immediate thoughts:

• Would need to template all occurrences of the Eigen Xd classes so that they could be Xcd instead.
• Would have to overload fill_with_random_normals to properly sample in a complex space.
• There's a few places where I assume that an operation (e.g., norm) produces a double. This may or may not be true for complex matrices/vectors, depending on the operation. Would need to check.

Assuming that the Eigen complex multiplications behave as expected, the rest of it should be fine as it's just Eigen-in, Eigen-out; I don't think I have many double pointers flying around.

That said, I don't know whether the IRLBA theory is trivially extensible to complex numbers. Someone with more expertise would have to comment here, maybe @bwlewis might have some thoughts.

PS - Any chance you can add a license file to the repo?

Done.

[bwlewis]

The method works fine with complex values (indeed, the original R method explicitly supports this)--just be careful to use conjugate transpose everywhere instead of just the real-valued transpose.

Things should work in lower precision, but I don't have any experience with that, in particular the convergence of the Lanczos process used. There is some literature on that, so I suspect it to work it's just that convergence might take a bit longer perhaps?

[LTLA]

Thanks @bwlewis - happily enough, I was already using Eigen's adjoint() instead of transpose(), so that should automatically do a conjugate transpose if the real matrix is swapped with its complex counterpart.

@spinicist It seems feasible, and possibly even straightforward, to add complex support to this library. However, to be honest with you, this won't be a high priority for me. You seem to have some experience with Eigen, so perhaps you might think about a PR? If you can make a minimal working example, I can add the tests and docs and such later.

[spinicist]

Great - I will add a PR to my to-do list. I am pretty sure .norm() is fine for complex, it's essentially defined as sqrt(dot(x, x)) and the Eigen dot product is already adjoint in the second argument. The only thing I don't immediately know how to do is the random numbers.

Thanks for adding the license and picking a liberal one.

[spinicist]

Hello,

I've finished a first pass at this here: #5

It turned out to be a more invasive change than I expected, as every helper class needed to be templated as well. I am not sure I've done this in the best way so I'm very open to suggestions for improvements.

I added a simple test which passes. I haven't had time to actually try this in my project yet.

[LTLA]

It turned out to be a more invasive change than I expected, as every helper class needed to be templated as well.

No worries, this was to be expected.

I am not sure I've done this in the best way so I'm very open to suggestions for improvements.

I made some edits to automate most of the templating based on the input matrix type. It has the bonus of being compatible with existing calls, so downstream users won't have to change any code.

I will say, though, that the extra templating is pretty rough on the compiler. tests/src/irlba.cpp takes a minute to compile on my laptop, using up to 4 GB of RAM. Don't remember it taking that long before, but oh well.

@bwlewis Is there something special we should be doing with the random Normal sampling in the complex case? From looking at irlba::irlba, the initial V vector is generated with rnorm and is fully real, but should it have a random imaginary component as well? Seems kinda odd that our initial vector would have all-zero imaginary when the final result is complex.