Repository navigation
complex numbers with vars #123
Description
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 typeT. We can either autodiff the real component, the imaginary component, or the result of a function that takes astd::complex<T>and returns aT. This means we don't have to extend autodiff to deal with adding a.grad()function tostd::complex<T>and we can focus on making surestd::complex<T>compiles and behaves appropriately. -
std::complex<T>whereTis 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 ofstd::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 toT()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 avarion the stack. In order for this to work, we would need to replaceT()withT(0)or some other way of instantiating an autodiff variable on the stack when instantiating a newstd::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 onstan::math::complex<stan::math::var>which then relied onstd::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
stdnamespace likeisnan(std::complex<T>)andisinf(std::complex<T>)because Eigen relies on them being in thestdnamespace. @ChrisChiasson was able to do this using ADL by adding these functions to thestan::mathnamespace, but we've already gone ahead and implementedstd::isnan(T)for both our autodiff types.
Activity
syclik commentedon Jul 18, 2015
@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 commentedon Jul 18, 2015
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 [email protected]
wrote:
bob-carpenter commentedon Jul 18, 2015
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 commentedon Jul 18, 2015
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 [email protected]
wrote:
bob-carpenter commentedon Jul 18, 2015
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!
bgoodri commentedon Jul 18, 2015
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 [email protected]
wrote:
betanalpha commentedon Jul 19, 2015
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 [email protected] wrote:
bob-carpenter commentedon Jul 19, 2015
I wasn't even considering differentiating w.r.t. the
complex part --- no wonder I couldn't understand anything
I was seeing :-)
52 remaining items
Fixed operator* and operator/ implementation and tests (Issue #123)
Added specializations of operator!= to complex<var> to avoid segfault…
Added specializations for std::complex<stan::math::var> operator*= an…
Moved all Bgoodrich's tests over. Added sqrt and pow (Issue #123)
Made mingw implementation of pow inline (issue #123)
bob-carpenter commentedon Mar 18, 2022
This has been resolved.