complex numbers with vars

[syclik]

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 Feature/issue 123 complex numbers with vars #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]

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

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
#123 (comment).

[bob-carpenter]

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]

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::complex 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.

—
Reply to this email directly or view it on GitHub
#123 (comment).

[bob-carpenter]

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::complex 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.

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

[bgoodri]

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::complex 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.

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

—
Reply to this email directly or view it on GitHub
#123 (comment).

[betanalpha]

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::complex 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.

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

[bob-carpenter]

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::complex 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.

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

—
Reply to this email directly or view it on GitHub
#123 (comment).

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

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

[bgoodri]

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]

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.

[bob-carpenter]

This has been resolved.