complex numbers with vars

syclik commented 9 years ago

From @bgoodri on July 18, 2015 15:10

The following C++ program compiles but segfaults (with clang++-3.6 and g++-5)

#include <stan/math.hpp>
#include <iostream>
#include <complex>

int main() {
  stan::math::var x = 1;
  std::complex<stan::math::var> z = x;
  std::cout << "imaginary part is " << z.imag().val() << std::endl;
  return 0;
}

I think we should try to prevent that from compiling, or at least dying gracefully, but failing that I guess we should initialize the imaginary part to zero?

Copied from original issue: stan-dev/stan#1554

There were a lot of things learned by @ChrisChiasson's effort to implement std::complex for Stan's autodiff variables in PR #789. Some of the important takeaways:

The scope is for both the real and imaginary parts of std::complex<T> to be of type T. We can either autodiff the real component, the imaginary component, or the result of a function that takes a std::complex<T> and returns a T. This means we don't have to extend autodiff to deal with adding a .grad() function to std::complex<T> and we can focus on making sure std::complex<T> compiles and behaves appropriately.
std::complex<T> where T is our autodiff types (e.g. stan::math::var, stan::math::fvar<U>) are explicitly unspecified (as of C++11, 08/2018). "The effect of instantiating the template complex for any other type is unspecified." https://en.cppreference.com/w/cpp/numeric/complex
Given that std::complex<T> is unspecified, it is a requirement for any implementation to have tests for the basic operations of std::complex<T>. The list of functions can be found here: https://en.cppreference.com/w/cpp/numeric/complex. That means that we should be testing for construction, the member functions, the non-member functions, etc. It's ok to do this in small chunks, but since we're relying on unspecified behavior, we need to be able to tell when our implementation fails in the basic tests before trying to push this through Eigen functions.
The core reason why std::complex<T> can't just be instantiated with one of our autodiff types is because in the current implementation, there are calls to T() to construct real or imaginary components. In our current autodiff design, stan::math::var() results in an uninitialized variable. This is intentional so we don't create a vari on the stack. In order for this to work, we would need to replace T() with T(0) or some other way of instantiating an autodiff variable on the stack when instantiating a new std::complex<T>. For specific places where this happens (especially in the default constructor), search for _Tp() in the implementation file: https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/std/complex
@ChrisChiasson's PR #789 cleverly worked around the above problem. std::complex<stan::math::var> was implemented by relying on stan::math::complex<stan::math::var> which then relied on std::complex<stan::zeroing<stan::math::var>>. stan::zeroing<stan::math::var> was a type that had a default constructor that actually made it to the stack. There was a lot of template magic to make this work properly. There is a prototype of this in the closed PR. In order to finish out this design, the full design would need to be written up (it's tricky templating and it extends our meta traits) and fully tested.
An alternative design might be to write explicit template instantiations for our types for places where it constructs uninitialized autodiff variables for different components. See this comment for an example of how to template specialize the constructor. If this works, we may only need to specialize a few functions within std::complex<T> and not have to use tricky templating. (This is speculation and should be verified with tests.)
We may need to implement functions in the std namespace like isnan(std::complex<T>) and isinf(std::complex<T>) because Eigen relies on them being in the std namespace. @ChrisChiasson was able to do this using ADL by adding these functions to the stan::math namespace, but we've already gone ahead and implemented std::isnan(T) for both our autodiff types.

syclik commented 9 years ago

@bgoodri, do you know if it's a default constructor issue? Can you point a link to the doc for complex, specifically what constructor it uses and the contract for classes?

bgoodri commented 9 years ago

According to

http://en.cppreference.com/w/cpp/numeric/complex

The specializations std::complex, std::complex, and std:: complex are LiteralTypes http://en.cppreference.com/w/cpp/concept/LiteralType for representing and manipulating complex numbers http://en.wikipedia.com/wiki/Complex_number.

The effect of instantiating the template complex for any other type is unspecified.

But even if it were defined, we don't have any implementations for anything with complex numbers, so who knows what would result if someone tried to utilize them (which Eigen does occasionally).

On Sat, Jul 18, 2015 at 11:31 AM, Daniel Lee notifications@github.com wrote:

@bgoodri https://github.com/bgoodri, do you know if it's a default constructor issue? Can you point a link to the doc for complex, specifically what constructor it uses and the contract for classes?

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122554346.

bob-carpenter commented 9 years ago

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complex<var> constructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

bgoodri commented 9 years ago

When Eigen just extracts the real and imagainary parts to do something with them, we would probably be okay. But I would worry about our autodiff doing the right thing with complex.

On Sat, Jul 18, 2015 at 12:00 PM, Bob Carpenter notifications@github.com wrote:

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complexconstructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122558670.

bob-carpenter commented 9 years ago

Do complex numbers have different rules for derivatives than just autodiffing the real and imaginary components? I got lost in the theory straightway, but this looks promising:

https://en.wikipedia.org/wiki/Differentiation_rules#Elementary_rules_of_differentiation

If we can just apply the usual chain rule, we'd be OK, right? I know a lot of people reading this know complex analysis!

On Jul 18, 2015, at 9:21 AM, bgoodri notifications@github.com wrote:

When Eigen just extracts the real and imagainary parts to do something with them, we would probably be okay. But I would worry about our autodiff doing the right thing with complex.

On Sat, Jul 18, 2015 at 12:00 PM, Bob Carpenter notifications@github.com wrote:

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complexconstructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122558670.

— Reply to this email directly or view it on GitHub.

bgoodri commented 9 years ago

It might work out-of-the-box for some things but not others. For example, log() is defined for any non-negative real number but any complex number. So, we would have to do template specializations for a lot of functions and worry about things like branch cuts.

On Sat, Jul 18, 2015 at 2:04 PM, Bob Carpenter notifications@github.com wrote:

Do complex numbers have different rules for derivatives than just autodiffing the real and imaginary components? I got lost in the theory straightway, but this looks promising:

https://en.wikipedia.org/wiki/Differentiation_rules#Elementary_rules_of_differentiation

If we can just apply the usual chain rule, we'd be OK, right? I know a lot of people reading this know complex analysis!

On Jul 18, 2015, at 9:21 AM, bgoodri notifications@github.com wrote:

When Eigen just extracts the real and imagainary parts to do something with them, we would probably be okay. But I would worry about our autodiff doing the right thing with complex.

On Sat, Jul 18, 2015 at 12:00 PM, Bob Carpenter < notifications@github.com> wrote:

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complexconstructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122558670.

— Reply to this email directly or view it on GitHub.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122572067.

betanalpha commented 9 years ago

GIANT RABBIT HOLE WARNING.

Derivatives in complex spaces are actually more constrained that derivatives in real spaces, which makes them both easier and harder to use. As Ben alluded to, once you have to worry about branch cuts and the like your life gets equally more difficult and more fun.

But if we simplify things and just consider a complex number as a tuple of two real functions then we can naively apply autodiff and be fine (mathematically this corresponds to differentiating with respect to real variables only).

On Jul 18, 2015, at 7:36 PM, bgoodri notifications@github.com wrote:

It might work out-of-the-box for some things but not others. For example, log() is defined for any non-negative real number but any complex number. So, we would have to do template specializations for a lot of functions and worry about things like branch cuts.

On Sat, Jul 18, 2015 at 2:04 PM, Bob Carpenter notifications@github.com wrote:

Do complex numbers have different rules for derivatives than just autodiffing the real and imaginary components? I got lost in the theory straightway, but this looks promising:

https://en.wikipedia.org/wiki/Differentiation_rules#Elementary_rules_of_differentiation

If we can just apply the usual chain rule, we'd be OK, right? I know a lot of people reading this know complex analysis!

On Jul 18, 2015, at 9:21 AM, bgoodri notifications@github.com wrote:

When Eigen just extracts the real and imagainary parts to do something with them, we would probably be okay. But I would worry about our autodiff doing the right thing with complex.

On Sat, Jul 18, 2015 at 12:00 PM, Bob Carpenter < notifications@github.com> wrote:

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complexconstructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122558670.

— Reply to this email directly or view it on GitHub.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122572067.

— Reply to this email directly or view it on GitHub.

bob-carpenter commented 9 years ago

I wasn't even considering differentiating w.r.t. the complex part --- no wonder I couldn't understand anything I was seeing :-)

Bob

On Jul 18, 2015, at 8:03 PM, Michael Betancourt notifications@github.com wrote:

GIANT RABBIT HOLE WARNING.

Derivatives in complex spaces are actually more constrained that derivatives in real spaces, which makes them both easier and harder to use. As Ben alluded to, once you have to worry about branch cuts and the like your life gets equally more difficult and more fun.

But if we simplify things and just consider a complex number as a tuple of two real functions then we can naively apply autodiff and be fine (mathematically this corresponds to differentiating with respect to real variables only).

On Jul 18, 2015, at 7:36 PM, bgoodri notifications@github.com wrote:

It might work out-of-the-box for some things but not others. For example, log() is defined for any non-negative real number but any complex number. So, we would have to do template specializations for a lot of functions and worry about things like branch cuts.

On Sat, Jul 18, 2015 at 2:04 PM, Bob Carpenter notifications@github.com wrote:

Do complex numbers have different rules for derivatives than just autodiffing the real and imaginary components? I got lost in the theory straightway, but this looks promising:

https://en.wikipedia.org/wiki/Differentiation_rules#Elementary_rules_of_differentiation

If we can just apply the usual chain rule, we'd be OK, right? I know a lot of people reading this know complex analysis!

On Jul 18, 2015, at 9:21 AM, bgoodri notifications@github.com wrote:

When Eigen just extracts the real and imagainary parts to do something with them, we would probably be okay. But I would worry about our autodiff doing the right thing with complex.

On Sat, Jul 18, 2015 at 12:00 PM, Bob Carpenter < notifications@github.com> wrote:

This would make a great project for someone to explore over a summer. The first thing to do is to overload the std::complexconstructor to do the right thing with real and imaginary parts (yes, init imaginary to 0 with one input) and then see what happens in Eigen.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122558670.

— Reply to this email directly or view it on GitHub.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-122572067.

— Reply to this email directly or view it on GitHub.

— Reply to this email directly or view it on GitHub.

bgoodri commented 8 years ago

I think we should put off the idea of auto-diffing complex vars but fix the original issue that Stan Math can be compiled a std::complex<stan::math::var>.

betanalpha commented 8 years ago

Exposing a complex type which is just a pair would be handy as we could then expose FFTs.

On Mar 4, 2016, at 10:38 PM, bgoodri notifications@github.com wrote:

I think we should put off the idea of auto-diffing complex vars but fix the original issue that Stan Math can be compiled a std::complexstan::math::var.

— Reply to this email directly or view it on GitHub.

davidmanheim commented 7 years ago

This should allow using eigen for finding arbitrary eigenvalues, (i.e. for nonsymmetric matrices,) right?

So if/when this project is done, it would be nice to provide that functionality to STAN as well.

bob-carpenter commented 7 years ago

@davidmanheim We could already implement something that returned the complex and real components. The complexity comes in from trying to extend derivatives for functions of complex variables.

davidmanheim commented 7 years ago

If you can already make a change to the current implementation to allow finding eigenvalues to nonsymmetric matrices, that would be helpful for me at least...

bob-carpenter commented 7 years ago

@davidmanheim --- We could probably do that. Would you mind opening a separate more limited issue for just eigenvalues of non-symmetric matrices? What do you expect to get as a return (that is, what's the signature as a Stan function in terms of argument and return types) and what do you expect the derivatives to be w.r.t. the inputs.

bgoodri commented 7 years ago

We cannot do that because the eigenvalues of non-symmetric matrices are generally complex and we cannot autodiff that. Request the pseudo-eigendecomposition, which is (or could be) done with only real numbers.

On Mon, May 22, 2017 at 12:39 PM, Bob Carpenter notifications@github.com wrote:

@davidmanheim https://github.com/davidmanheim --- We could probably do that. Would you mind opening a separate more limited issue for just eigenvalues of non-symmetric matrices? What do you expect to get as a return (that is, what's the signature as a Stan function in terms of argument and return types) and what do you expect the derivatives to be w.r.t. the inputs.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-303154126, or mute the thread https://github.com/notifications/unsubscribe-auth/ADOrqiAOumLzW8G4Cff7PcdkDqwGf5urks5r8bpKgaJpZM4FbP69 .

bob-carpenter commented 7 years ago

What's the use case here? It might help to start with that and work down (I have a blog post on that scheduled for tomorrow!).

I assume the problem with returning complex numbers, even as two separate arrays of components, is that we can't just operate on the imaginary component without bringing in proper complex derivatives, which we don't have.

Are there situations where users are only going to care about the real component and can we autodiff through that and just ignore the complex compnents?

davidmanheim commented 7 years ago

My use case is actually fairly simple. In infectious disease modeling, I'm hoping to generate the largest eigenvalue of the non-symmetric "next generation matrix" - which, under reasonable assumptions, is always real - even if other components are not.

So yes, I only care about the real component. I'm unsure if that's a really unique use case, or if others would find this useful as well.

sakrejda commented 7 years ago

even in those population models the imaginary components of sub-dominant eigenvalues are related to transient oscillations so I don't know that there's much of an argument for ignoring them in general. It's easier to find the dominant eigenvalue than anything else so I wonder if a good compromise for David's use-case might be a function that just extracts the dominant eigenvalue (and throws if it's not real).

On Mon, May 22, 2017 at 2:55 PM David Manheim notifications@github.com wrote:

My use case is actually fairly simple. In infectious disease modeling, I'm hoping to generate the largest eigenvalue of the non-symmetric "next generation matrix" - which, under reasonable assumptions, is always real - even if other components are not.

So yes, I only care about the real component. I'm unsure if that's a really unique use case, or if others would find this useful as well.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-303189174, or mute the thread https://github.com/notifications/unsubscribe-auth/AAfA6UXGUs_yAbyb4_042b4CRNX49EtGks5r8do7gaJpZM4FbP69 .

bgoodri commented 7 years ago

It is difficult to ignore the imaginary component because it causes a segfault. I think we should static_assert on the var constructor for complex double once we get to C++11. The pseduo-eigendecomposition suffices for the "I am only interested in the real eigenvalue(s)" use case and the "All the eigenvalues are real, I promise" use case.

On Mon, May 22, 2017 at 2:48 PM, Bob Carpenter notifications@github.com wrote:

What's the use case here? It might help to start with that and work down (I have a blog post on that scheduled for tomorrow!).

I assume the problem with returning complex numbers, even as two separate arrays of components, is that we can't just operate on the imaginary component without bringing in proper complex derivatives, which we don't have.

Are there situations where users are only going to care about the real component and can we autodiff through that and just ignore the complex compnents?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-303187155, or mute the thread https://github.com/notifications/unsubscribe-auth/ADOrqnNEo6VK36DZQUXxIBPEE2YiPmAZks5r8dhigaJpZM4FbP69 .

bbbales2 commented 7 years ago

How much of a C++ fan are you David? There's some stuff with cmdstan that you can hook in some C++ code that interfaces with the autodiff. I've used it a couple times now and it works pretty slick.

If the parameter here is real, then the derivative of the largest eigenvalue with respect to this parameter should be possible to get at analytically (eq. 14 here http://www.win.tue.nl/analysis/reports/rana06-33.pdf).

You could make sure it's real on the fly like Krzysztof said. What's the matrix size?

Ben

On Mon, May 22, 2017 at 1:14 PM, bgoodri notifications@github.com wrote:

It is difficult to ignore the imaginary component because it causes a segfault. I think we should static_assert on the var constructor for complex double once we get to C++11. The pseduo-eigendecomposition suffices for the "I am only interested in the real eigenvalue(s)" use case and the "All the eigenvalues are real, I promise" use case.

On Mon, May 22, 2017 at 2:48 PM, Bob Carpenter notifications@github.com wrote:

What's the use case here? It might help to start with that and work down (I have a blog post on that scheduled for tomorrow!).

I assume the problem with returning complex numbers, even as two separate arrays of components, is that we can't just operate on the imaginary component without bringing in proper complex derivatives, which we don't have.

Are there situations where users are only going to care about the real component and can we autodiff through that and just ignore the complex compnents?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-303187155, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADOrqnNEo6VK36DZQUXxIBPEE2YiPmAZks5r8dhigaJpZM4FbP69

.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-303207817, or mute the thread https://github.com/notifications/unsubscribe-auth/AEhdGMNyX7mo-8t4TakekLdRv95Rwz93ks5r8eywgaJpZM4FbP69 .

bob-carpenter commented 6 years ago

I didn't think you could just plug an autodiff type into std::complex<T> as that first paper suggests (it uses ADOL-C's adouble type). See, e.g.,the Wikipedia definition of the complex chain rule or this math overflow post.

bob-carpenter commented 6 years ago

I think given the discussion above from @betanalpha that this will be OK if you only need derivatives w.r.t. the real parts.

I'm not clear on what a bunch of that C++ is doing, such as the zeroing template and its specializations---given that they only have typedefs, I'm assuming it's some kind of traits metaprogram, but I don't see where the typedefs get used anywhere.

I also didn't understand the purpose of zvar, which just seems to mark var for some reason (there's no other interface or implementation beyond its superclass's).

bob-carpenter commented 6 years ago

Our documentation's still way short of where it should be for onboarding new devs. Let us know if you need help. And we'd love to hear back if you get it working on an interesting problem you can share.

There's an arXiv paper that goes over the whole structure of our autodiff if you haven't found that yet.

seantalts commented 6 years ago

Yes please! Thanks! On Wed, Mar 7, 2018 at 13:03 ChrisChiasson notifications@github.com wrote:

When I try to push my feature branch:

chris.chiasson@ MINGW64 ~/workspace/stan-dev/math (feature/issue-123-complex-numbers-with-vars) $ git push origin feature/issue-123-complex-numbers-with-vars ERROR: Permission to stan-dev/math.git denied to ChrisChiasson. fatal: Could not read from remote repository.

Please make sure you have the correct access rights and the repository exists.

I am sure I'm doing something wrong because I don't usually contribute to other projects... Should I fork your repo?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-371227903, or mute the thread https://github.com/notifications/unsubscribe-auth/AAxJ7MpagqdWcrOdhm6_QYuPH8s7sy5Qks5tcCDegaJpZM4FbP69 .

seantalts commented 6 years ago

I think Jenkins should auto-format that stuff soon, so no worries.

On Wed, Mar 7, 2018 at 3:32 PM, ChrisChiasson notifications@github.com wrote:

Sorry about what appear to be tabs in the diff... In my editor, tabs are shown as single spaces, and I guess I forgot to tell it to always save files with spaces only.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/123#issuecomment-371274728, or mute the thread https://github.com/notifications/unsubscribe-auth/AAxJ7MFct1vSJZd7s6hNtTpskiJvhbmvks5tcEPNgaJpZM4FbP69 .

syclik commented 6 years ago

I just updated the top-level comment. @bgoodri and @ChrisChiasson, could you check to see if that's appropriate? Were there things I missed? @seantalts, mind taking a look to see if that's clear enough to keep the issue open? I just reopened it.

bgoodri commented 6 years ago

I think the original issue 123 is basically moot. It segfaults, we once tried to prevent it from compiling with a std::complex<stan::math::var>, we failed to prevent it from compiling, but nothing you can do in the Stan language results in a std::complex<T>.

If we wanted to reopen #789 and support complex operations in Stan Math, that would be a feature rather than a bugfix.

bbbales2 commented 6 years ago

The list of functions can be found here

I had a look at the functions that need tested. I can do some work on these. I presume it's just simply checking values + gradients for each one?

What branch should I work from? #789 is showing "unknown repository"

bbbales2 commented 6 years ago

Also I'll say it again here just so it's written down somewhere, but I like the idea of only implementing std::complex<stan::math::var>. I don't think I'd need std::complex<stan::math::fvar> or the more complicated stuff for what I want to do.

bob-carpenter commented 6 years ago

What do you want to do?

We're going to have a bunch of things already that we're going to have to disallow for samplers or optimizers using higher-order autodiff.

On Aug 10, 2018, at 1:59 PM, Ben Bales notifications@github.com wrote:

Also I'll say it again here just so it's written down somewhere, but I like the idea of only implementing std::complex. I don't think I'd need std::complex or the more complicated stuff for what I want to do.

bob-carpenter commented 2 years ago

This has been resolved.

stan-dev / math

complex numbers with vars #123