Fully functional std::complex specializations and overloads

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.

BTW, I should probably clarify this just in case: the branch ChrisChiasson-feature/issue-123-complex-numbers-with-vars was created by someone else, not me, even though it has my name on it.
The branch I made that appears in PR 789’s diffs can be checked out out with the following commands:

git fetch origin pull/789/head:sometempbranchname
git checkout sometempbranchname

This part of its tests shows int->double->complex promotion & interaction with var. There is also a similar one for fvar.

I think the only matrix decomposition test I wrote was optimization related, because I wanted to be sure that the derivative information would be carried properly all the way through Eigen’s eigenvalue decomposition. It 2nd order minimizes–via Hessian inverse dot negative gradient–the difference between the real and imaginary parts of the first eigenvalue of a rotation matrix. It starts at a rotation angle of 3pi/8 and tests that the optimization stops at pi/4 (where the real and imaginary parts of the eigenvalue are equal). During every iteration, it also checks the entries of the matrix decomposition by the identity for an eigensystem.

For a long while, these were the only tests in my branch. IMO, these three tests would be worth including in yours, because they test the interaction between the different types of variables you mentioned, and because they test that the derivative information is perserved even during the matrix manipulations that Eigen performs.

I just wanted to say that I am glad it looks like we will have a version everybody is happy with, and I am glad that @ChrisChiasson stayed with us for the whole journey and is still contributing constructively to the effort. Also thank everybody else who has worked hard to figure out the best way to involve complex numbers in Stan over the time.

Thanks for the correction—I updated the original post.

I figured someone would insist on testing the derivatives of those complex matrix decomposition functions :-). So far, I’ve just been checking values and assuming autodiff will work through everything as there aren’t any custom derivatives anywhere.

I’ll change that. With the autodiff testing framework and C++11, it’s easy to write thorough unit-level autodiff tests at all autodiff levels. For example, here’s the test I just added for real components of eigenvectors:

template <typename T>
Eigen::EigenSolver<T> eigen_solver(const T& a) {
  Eigen::EigenSolver<T> s(a);
  return s;
}
...
TEST(...)
  auto f1 = [](const auto& a) {
    return eigen_solver(a).eigenvectors().real().eval();
  };
  Eigen::MatrixXd a(2, 2);
  a << 1, 2, 3, 0.7;

  stan::test::ad_tolerances tols;
  tols.hessian_hessian_ = 1e-2;
  tols.hessian_fvar_hessian_ = 1e-2;
  stan::test::expect_ad(tols, f1, a);
}

This’ll get similar tests for eigenvalues and imaginary components.

The lowered tolerances are for Hessian finite difference relative error comparisons. I’m thinking the default may need to be 1e-2 or 5e-3 as there are a lot of functions where I have to lower the tolerances like this. Or maybe I could replace the double finite diffs with finite diff on the gradient—all the tests would still be well founded and there should be much more accuracy.

I needed the factory function eigen_solver to do type inference for the class template parameter T for Eigen::EigenSolver<T>. Does anyone know if Is there an easier way to just pull the type off an auto variable? I used the old-school factory pattern because instantiating EigenSolver<decltype(a)> for an const auto& a function argument was giving me incomplete type errors.

I want to test all the binary functions for all combinations of int, double, complex<double>, T, and complex<T> inputs for all autodiff types T up to fvar<fvar<var>>. For unary functions, we only need to test complex<T> as the other tests are elsewhere. Each test is against double, so those all get exercised to make sure they can be instantiated and provide the same values and finite diff derivatives.

The expect_identity function should also be checking that the Jacobian is a zero matrix.

I’m inclined to just get rid of these algebra-driven tests and stick with the direct function tests through the autodiff test framework. They’re nice shortcuts to get off the ground, but with brute-force autodiff testing, I don’t think they’re adding anything.

Well, I intended all those tests to be evaluated across a bunch of random matrices in order to find out whether there are places where the accuracy gets unacceptable. But people hated non-deterministic unit tests.

lol, sorry :) The reason I tested them in mine is that I didn’t know exactly what Eigen was doing, specifically if there were any numbers initialized through identities they exploited in code without calculation, possibly breaking the adjoint flow. Might not be possible anyway; didn’t know / wanted to program defensively.

Hrm, if a is a reference, maybe the result of decltype(a) could be post-processed via decay or remove_reference?

Are the tolerances treated as absolute residuals (abs(left-right)<tol), as relative ones (abs(left-right)<tol*max(abs(left),abs(right))), or combined like (abs(left-right)<tol*max(1.0,max(abs(left),abs(right))))? I usually use a combined type of check to handle cases with large and small magnitudes, with tol around half the precision of a double (about 1.0e-8, or more closely pow(2.,-53./2.)). If it still isn’t loose enough, I provide a multiplier on top of tol. A combined approach should allow the tol to be a fairly small number unless the finite difference is in a really bad spot.

Just think how far we’ve come from when people were debating about whether or not it was even possible to have stan do this at all…

I’m reading your code more closely today before swapping out my implementation for yours in my application at work.

Right now, the inheritance structure I’m reading is:
std::complex<stan::math::var> inherits from stan::math::complex_base<stan::math::var>. The operations of stan::math::complex_base statically downcast their return references back to std::complex<stan::math::var>.

For reference, my inheritance structure is:
std::complex<stan::math::var> inherits from stan::math::complex<stan::math::var>, which in turn inherits from std::complex<stan::math::zeroing<stan::math::var>>>. There is also an operator std::complex() defined on stan::math::complex, which implicitly statically downcasts to std::complex for both base classes. (The extra level of inheritance together with the implicit downcasting delegates all operations back to the standard library as if they were operating on a zeroing version of stan’s types. This is why there is so little code, excluding the Jenkins autoformat misconfiguration that had touched a lot of files.)

While trying to test your code:
I was trying to build the libraries. I ended up partially following a page I wrote in the math wiki a while back, along with an example Jenkins test run. make test-headers fails with the error below. ./runTests.py fails because due to a missing Google Test library. I’m on commit 2c746c2cc9d19feb036dae474444b86d940be210 (HEAD -> feature/0123-complex-var, origin/feature/0123-complex-var). Edit: Forgot to add that I’m on Windows with MSYS2 x86_64 using GCC 9.1.

g++ -std=c++1y -m64 -D_REENTRANT -Wall -Wno-unused-function -Wno-uninitialized -Wno-unused-but-set-variable -Wno-unused-variable -Wno-sign-compare -Wno-unused-local-typedefs      -I lib/tbb_2019_U8/include -O3  -I . -I lib/eigen_3.3.3 -I lib/boost_1.69.0 -I lib/sundials_4.1.0/include  -D_USE_MATH_DEFINES  -DBOOST_DISABLE_ASSERTS      -c -std=c++1y -m64 -D_REENTRANT -Wall -Wno-unused-function -Wno-uninitialized -Wno-unused-but-set-variable -Wno-unused-variable -Wno-sign-compare -Wno-unused-local-typedefs      -I lib/tbb_2019_U8/include -O3  -I . -I lib/eigen_3.3.3 -I lib/boost_1.69.0 -I lib/sundials_4.1.0/include -O0 -include stan/math/mix/mat/fun/typedefs.hpp test/dummy.cpp -o nul
In file included from ./stan/math/rev/core.hpp:49,
                 from ./stan/math/mix/mat/fun/typedefs.hpp:6,
                 from <command-line>:
./stan/math/rev/core/std_complex.hpp: In instantiation of 'std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&) [with T = stan::math::var]':
./stan/math/rev/core/std_complex.hpp:693:48:   required from 'std::complex<_Tp> stan::math::complex_asinh(const std::complex<_Tp>&) [with V = stan::math::var]'
./stan/math/rev/core/std_complex.hpp:1537:25:   required from here
./stan/math/rev/core/std_complex.hpp:229:23: error: no matching function for call to 'value_of_rec(stan::math::complex_base<stan::math::var>::value_type)'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |           ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note: candidate: 'template<class T> std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&)'
  228 | inline std::complex<double> value_of_rec(const std::complex<T>& z) {
      |                             ^~~~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note:   template argument deduction/substitution failed:
./stan/math/rev/core/std_complex.hpp:229:23: note:   'stan::math::complex_base<stan::math::var>::value_type' {aka 'stan::math::var'} is not derived from 'const std::complex<_Tp>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |           ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:229:47: error: no matching function for call to 'value_of_rec(stan::math::complex_base<stan::math::var>::value_type)'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                   ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note: candidate: 'template<class T> std::complex<double> stan::math::value_of_rec(const std::complex<_Tp>&)'
  228 | inline std::complex<double> value_of_rec(const std::complex<T>& z) {
      |                             ^~~~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:228:29: note:   template argument deduction/substitution failed:
./stan/math/rev/core/std_complex.hpp:229:47: note:   'stan::math::complex_base<stan::math::var>::value_type' {aka 'stan::math::var'} is not derived from 'const std::complex<_Tp>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                   ~~~~~~~~~~~~^~~~~~~~~~
./stan/math/rev/core/std_complex.hpp:229:57: error: could not convert '{<expression error>, <expression error>}' from '<brace-enclosed initializer list>' to 'std::complex<double>'
  229 |   return {value_of_rec(z.real()), value_of_rec(z.imag())};
      |                                                         ^
      |                                                         |
      |                                                         <brace-enclosed initializer list>

Thanks. I hadn’t tested the headers—those are always easy to fix with more includes, so I’ll fix them and check in a cleaner version. Aside from failing the header tests, were there other problems?

The reason I was reluctant to delegate back to the standard library is that instantation of std::complex<T> for anything other than a primitive type is undefined behavior.

I just used a simple CRTP base class to avoid shared code between forward and reverse std::complex implementations, figuring that cutting-and-pasting a whole interface wouldn’t get by code review.
And I wanted to support mixed operations of things like std::complex<double> and std::complex<var> without promotion. Eventually, if we care about performance, we’ll need to reimplement reverse mode with analytic derivatives, not just autodiff.

I didn’t find the need to zero anything anywhere other than in the std::complex constructors.

I’ll give it a shot. I just took the easy way out because I was wrestling with so many type errors on the way to getting everyting to compile.

The errors are relative unless one of the values is very small at which point it kicks over to absolute errors. The tolerance is set around 1e-8 for gradients, but by the time you get to Hessians with finite differences (and they’re based on three finite diffs, not just one in a band, which is partly why they’re so slow), things are much less stable arithmetically. The tolerance there is 1e-3 and often has to go down to 5e-3 or 1e-2. Longer term, I think the right thing to do is instead do finite differences of gradients, since the gradients are being independently tested. Same for third order derivatives, but those are so often zero that it’s not nearly as much of a problem for testing as Hessians.

We could also use more of the algebraic tests if people are concerned about precision. I just figured if Hessians were right up to 1e-2 they were probably actually right. That’s obviously not a watertight line of reasoning.

I don’t like that they can fail randomly and I don’t like that the coverage is random. I’d rather have someone think out the edge cases and test those.

I’m one of those people. The reasons are that They’re unruly in that they can fail randomly, especially with finite diffs. Isn’t it also hard to actually get edge tests like highly ill-conditioned matrices randomly?

I don’t like that the test coverage is a measure zero subset of the parameter space it is going to get used for with HMC. In general, our unit tests primarily test whether stuff compiles and whether the analytical derivatives are logically correct, but there is very little testing of whether we are getting the reasonable numerical values considering that numerical accuracy is critical to HMC. Checking that you get reasonable numerical values for a small matrix often with integer values does not tell you whether you are going to get reasonable numerical values when plugged into an actual Stan program. We have some experience as to what sort of things tend to go badly in Stan programs, but no accumulated experience as to what is going to happen when you extend out into the complex plane.

To test “edge” cases, you usually have to construct things manually, but it is not as if you need an edge case like a singular matrix for things to break down numerically. Often it can be just due to having big values in some cells and small values in other cells, particularly if the matrix is large. To find things like that out, it is often good to use random matrices, which is what Eigen largely does in its unit tests with basic scalar types.

What about using randoms, but if it hits an error, making it dump the input to enable devs to add a new test with exactly those inputs that gets run every time? Over the years, it might accumulate some really good test cases just by chance.

Hrm, sounds reasonable to me.

+1

Hrm, although I would note in the first instance that the goal of both implementations is already in conflict with the ‘rule’: “Specializing the template std::complex for any type other than float, double, and long double is unspecified.”

Of course, I never follow the rules anyway. Like using auto all the time with Eigen… Eh, to each his or her own. I’m just happy someone in a position to do something about it is adding stuff I need to stan.

I got further more recently by checking out master from develop, building a core rev test to generate the libs as a side effect, and then switching back to the complex branch. To integrate with my app, I had to figure out MinGW with both static and shared libs (tbb)… (side note, rpath does not work on Windows, and TBB likes to spawn 3 threads just for a single threaded unit test… ::arches eyebrow::).

I’ll check when I get in to work tomorrow, but it was telling me that pow(var,int) (I think that’s what it said) was ambiguous, an error that didn’t show up when using my implementation. I don’t have clang++ fully set up at work because the link step fails in MinGW, but I did copy down and email to myself the list of candidate functions from it (g++ doesn’t dump that info):


../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2568:26: note: candidate function
inline std::complex<var> pow(const var& lhs, const std::complex<double>& rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2585:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, int rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/scal/fun/pow.hpp:124:12: note: candidate function
inline var pow(const var& base, double exponent) {
           ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2589:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, double rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/scal/fun/pow.hpp:108:12: note: candidate function
inline var pow(const var& base, const var& exponent) {
           ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2593:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, const var& rhs) {
                         ^
../new-stan-dev/math\stan/math/rev/core/std_complex.hpp:2597:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs,
                         ^

I dont want to derail this topic but not sure what you mean by “not working”? In order to do dynamic linking on Windows TBB needs to be in the users or systems PATH variable.

For Stan Math we let the TBB threadpool handle the number of spawned threads with their default policy. You can override that by calling stan::math::init_threadpool_tbb();. If STAN_NUM_THREADS is set, that call will spawn the number of threads set in that variable, otherwise it will spawn 1 threads. Cmdstan calls that at the start of each model.

rpath, if I understood it correctly, is different from the PATH environment variable. I believe you’re correct about PATH (in that it is one of the ways to solve the problem). I had only mentioned rpath because it is in the instructions for how to compile a stan math program in the readme.md.

It was pow(var,int). The beginning of the error snippet was as follows:

src/Nautical/MotionTables/MotionTables.hpp:55:16: error: call of overloaded 'pow(stan::math::var, int)' is ambiguous
   55 |        ess+=pow(ERSAft-check,2);
      |             ~~~^~~~~~~~~~~~~~~~
compilation terminated due to -Wfatal-errors.
In file included from src/main.cpp:9:
./src/Nautical/MotionTables/MotionTables.hpp:55:13: fatal error: call to 'pow' is ambiguous
       ess+=pow(ERSAft-check,2);

I have been wondering, why are you calling it CRTP? The pattern I see is

template <>
class complex<stan::math::var>
    : public stan::math::complex_base<stan::math::var> {

I agree with that and consider it worthwhile. I would only note the approaches are not mutually exclusive. In my implementation, sometimes we selectively avoided delegating back to the standard library, especially with g++ (as you mentioned), by writing declarations with more specific types to give them higher priority in template selection and overload deduction.