The key issue is that std::complex<T> being instantiated is unspecified behavior. So what I’m providing is a complete specialization so nothing relies on unspecified behavior.
For a primer on the evils of function template specializations:
Overloads of missing std:: functions like isinf in namespace stan::math for argument-dependent lookup; these are not related to std::complex but necessary for some of Eigen’s algorithms; also specializes iterator_traits.
Specializations of std::complex<var> and std::complex<fvar<U>>. The constructors initialize real and imaginary values to zero where necessary (this is roughly what’s suggested in the last bullet point of @syclik’s summary on the issue).
Code duplication in the std::complex<T> class overloads was eliminated by factoring the class overloads into a CRTP-based base class. Code dupe in functions was eliminated by defining fully templated function templates for all complex operations in stan::math to which all other specializations and overloads delegate. I simplified signatures as much as possible to remove ambiguities while covering the combinatorics. There are notes on how that’s done as code comments, since it won’t be obvious from just reading the code (8 overloads should be required for ADL, but the specializations need to be redundantly overloaded to work in both g++ and clang++)
Specializations of std::complex functions where available. Most functions are not class member functions. This is necessary to get around libstdc++'s (used in g++) poor coding of std::complex for argument-dependent lookup—they hardcode std:: qualifiers on many complex functions. clang++ works without any specializations.
Overloads of all functions defined in std:: in namespace stan::math for all possible complex instantiations with autodiff variables and for all combinations of pairs of operands (int, double, std::complex<double>, Tandstd::complex, where Tis one of our autodiff types). These will be picked up by argument-dependent lookup. They will prevent overpromotion tocomplex` or to an autodiff type.
Some supporting metaprograms to deal with return types in complex template functions; these could potentially be merged into general things like return_type_t but I didn’t want to extend this mega-PR beyond what was necessary for std::complex (my plan’s to enlist @stevebronder’s help on integrating the metaprograms into the general Stan metaprogramming framework—for now they’re just coded to work simply when at least one complex type is involved)
Complete unit tests for constructors, standalone functions, and overloads. This uses the autodiff test framework to test derivatives at all levels of autodiff: var, fvar<double>, fvar<fvar<double>>, fvar<var>, and fvar<fvar<var>>. That makes sure everything gets instantiated with both complex<T> for autodiff type T and mixed with complex<double> and all other types like double, int, and T itself. Where autodiff’s not involved, operaitons like constructors and setters/getters are tested directly for reference equality.
The final tests are borrowed from the previous PR and instantiate three solvers at all levels of autodiff to make sure they get the right result
Eigendecomposition for asymmetric matrices
Pseudoeigendecomposition
Complex Schur decomposition
Staging a PR
I plan to break this monolithic branch down in order to create a PR in three steps.
namspace stan::math: basic Stan infrastructure that doesn’t depend on complex overloads
namespace std: specializations of std::complex class for var and fvar
namespace stan::math: template functions operating on std::complex
namespace std: specializations of std::complex functions plus namespace stan::math overloads of std::complex functions for var and fvar
I wrote a bunch more tests on one of the branches somewhere that tested whether the derivatives were accurate. Basically, for each (unary for simplicity) function, foo, I would evaluate it with a varx and for h = 1e-8 also evaluate it with a std::complex<var, var>(x, h). And then test whether the sensitivity to x for the former was near the imaginary part of the latter when divided by h. I will try to find those tests again, but it would be good if that could be meta-programmed somehow so that it would happen automatically.
complex<var>. It helps for me to think of a complex<var> as a shorthand for the function z\left(x,y\right) = x + iy with either (or both) of x and y being a (real) var and i = \sqrt{-1} a constant. So, then if you also had a (real) var, w, and wanted to autodiff some function of w and z\left(x,y\right) like w \times z\left(x,y\right) = wx + iwy, you should get back a complex<var>.
Had some free time tonight so I started goofing around with adding the stan math meta programming stuff. Not all the way done but I got pretty far with std_complex.hpp. I figure a complex<var> is a pretty heavy type so I’d be nice to have the perfect forwarding stuff here so we can avoid copies where possible
Done! As long as I test {x, 0} (that is, a complex number with imaginary component zero), I get exactly those tests scripted for x in the test framework for first, second, and third order derivatives through our functionals. This exercises types int, double, var, fvar<double>, fvar<fvar<double>>, fvar<var>, fvar<fvar<var>>. Each test case makes sure that if a function throws, it throws at each autodiff level; if it suceeds, it has to have the same value and derivatives that match finite differences for each of the above types. The base double versions are implemented using std::complex<double>, so we make sure they match the library implementations on each platform.
Of course! Didn’t you have an application that was driving you to do this in the first place? If you could test that, it’d be super helpful. As is, I only have versions of your tests for the asymmetric eigendecomposition, pseudoeigendecomposition, and complex Schur decomposition. So anything else you can suggest would be helpful. I just solved the type algebra problem for the classes—I don’t really understand complex numbers at much deeper a level than i = \sqrt{-1} and having a sense of what a dual-number algebra should do (too much pure discrete math in my background and not enough continuous applied math for my current job).
Then the next thing to think about is how to expose complex-based functions like the test cases in Stan. I have no idea what people want there, so I stopped at the test cases. One other function I know people want from this is fast fourier transforms. We use those to compute autocorrelations in Stan already using the Eigen unsupported implementation (which seems to work and be supported just fine).
The main thing to do for this code is code review, either formally or informally. So if you have comments on how to make the code better or simpler, I’d love to hear about it.
The code should look pretty familiar to you as you’d gotten a lot of the way there implementing the class specializations and overloads. I followed your implementation in terms of what would need to be implemented (like the iterator traits and std::isinf, etc.). I just went through and brute-forced every single specialization and overload for g++ and clang++ compatibility using the testing framework; see my reply to @bgoodri for more details on what’s being tested.
I had to bend over backward to fight with clang++'s missing templates and g++'s ADL-unfriendly implementation. The testing framework was really key in finagling the edge cases—don’t know how I could’ve done this writing test cases one by one at five levels of autodiff. I was very confused about why things weren’t working until I read the libstdc++ code from g++. The end result is that g++ takes fifteen minutes to compile the tests, whereas clang++ takes less than one. I haven’t tested run time performance of the generated code.
A complex<T> is two instances of T, and a var is a single pointer, so `complex is two pointers or 16 bytes.
We’ve been coding from the get-go as if the pass-by-value efficiency cut-off is 8 bytes. I just picked this notion up through folk wisdom about passing primitives by value and everything else by constant reference.
That’s how I learned it in high school and when we had to do examples for basic abstract algebra in college. It doesn’t help (me at least), when you have to figure out how pow or atanh are supposed to work. Euler, I aint. Luckily, the spec and the Wikipedia have all the necessay formulas, which are as simple to code as they are difficult to comprehend. For example, it was news to me that complex cosine is defined by
\textrm{cos}(z) = \textrm{cosh}(i \cdot z),
but it’s trivial to code (with a helper function i_times that’s more efficient than literally multiplying by i).
that is complicated. Maybe we could pass those off to Boost and do the derivatives analytically, but that could be a subsequent PR since I don’t think those functions are used much in statistics.
Thanks for the link. There’s definitely room for imporvement on the numerics. I just used the “textbook” implementations everywhere that were described in the spec as “function behaves as if implemented by …”.
It looks like the boost functions are templated appropriately that we could autodiff through them. You could probably just replace my own implementations in complex_asin, etc., which have exactly the same signatures as the Boost implementations. I should’ve looked there first!
It’ll definitely be easier to wait until everything’s working before trying to optimize, though the optimizations should be dead easy. In fact, now that I think about it a bit more, it might be easier to just plug their implementations in because then I wouldn’t have to test 25 of my own infrastructure functions and can stick to testing the interface functions that call them.
Hmm. They only supply the arc trig functions, but I can use those and a lot of the other ones depend on these. Let me give it a try.
Yes sir. I have an application at work that sits on top of my complex implementation. When I get back in to the office on Monday, I can try to swap your code underneath the app and tell you what blows up. It currently uses complex math in two areas: the calculation and optimization of response amplitude operators, and the calculation/optimization of some approximations to tapered beam & tapered spaceframe properties, where the approximations use complex numbers.
You found something missing from mine?
BTW, I read a little bit of your new complex<T> code this morning. I can see what you mean by needing the new autodiff testing framework to batch process the test cases. Since you actually implement the complex<T> functions rather than delegating back to the original functions in libc++/libstdc++ (like mine), the new code is actually switching function implementations when using Stan types vs double. I can see how this would be difficult without being able to rely on something that would compare the var output with the double output automatically. However, it does mean we would be more exposed to tricky parts of the function domains. (It amplifies the question people always have about test coverage: “Is this test passing because of where we are in the function domain, or is it passing because the function is actually correct?”)
Is the testing framework univariate in the independent variable? If so, maybe the framework could be proxied / multiplexed for test coverage into the imaginary part?
We could, although the analytical derivatives are typically very easy and better behaved numerically than the function values. You can get the derivative of asin, for example, from
The analytic derivatives should all be pretty easy because the functions are all pretty simple. I’d rather wait on doing the optimizations until the basic implementations are working.
Most of the function specializations and overloads were missing. It looks like there are implementations inside of g++ and clang++ that provide them, but it’s hard to tell where that is and isn’t going to work, so I just implemented everything. We’ll eventually want to drop in custom autodiff anyway for efficiency.
The main difference is that I never had to take extreme measures to set things to zero. I just defined the constructors in std::complex<var> and std::complex<fvar<T>> to do that.
I agree that’s always a problem with using only finitely many tests.
It can handle arbitrary types of up to two variables (three when the latest PR gets merged). But those can be any type. So I treat a univariate complex function as a function working on a vector of two elements and returning a std::vector of two elements. Then it all fits in the framework. I do that with the to_array and from_array functions. So all the gradients for both real and imaginary parts get tested.
We’d want to be testing all this anyway because we never know if an external implementation of a function will be autodiffable.
Specifically, the default implementations run into trouble when you combine double and complex<var>, or var and complex<double>, or complex<var> and complex<double>.
It looks to me like what was going on in @ChrisChiasson’s version is that only things where the core g++ and clang++ libs didn’t work were patched.
It would be good if we had a bunch of unit tests that integrated a function from \mathbb{R} to \mathbb{C} over the whole real line or at least the whole positive real line. Tests like that in
were really good for ensuring that our PDFs integrated to 1. Unless the inaccuracy cancels, that suggests the accuracy is good over the whole domain.
Yes, they’re provided by the C++ libraries in mine. I had thought you meant like something was verified broken in mine. Understood now.
Ok, that sounds great. At first I was getting concerned when I read about the zero imaginary part above.
It’s true; they will blow up if they’re not pre-empted carefully. We fixed everything we could find, though. The ctors and the auxillary functions (copysign, isinf, val, etc), were set up slightly differently than the current branch in order to enable that to happen. (?? Let me know if there is some specific spot that is wrong in mine-- The code is being used for offshore structures that get made in real life that mortal people use… Safety first…)
Even though it isn’t fool-proof, one thing that might be worth checking is that if we’re only going to test one point in a domain (not that we couldn’t test more), we might choose a spot that has non-zero 1st and 2nd derivatives in both directions at that location, to make sure we capture changes.