Fully functional std::complex specializations and overloads

Developers
21

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.

22

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.

23

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.

24

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

25

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.

26

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.

27

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.

28

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…

29

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>
30

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.

31

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?

32

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.

33

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,
                         ^
34

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.

35

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.

36

A finite set of random tests is still measure zero. How about rather than having unit tests, we run a million trials of something random under development? That way, we don’t have random inexplicable failures of unit tests unrelated to the branch being developed.

37

I called it CRTP-like in the doc. It would be an instance of CRTP if the template type were std::complex rather than T; but it’s used just like the CRTP in terms of casting to the subclass.

I’ll check out what’s going on. I may not have been testing that case.

I can try to fire up your PR and see if it passes all of my tests. I thought this was going to be a lot easier than it is.

38

Also, as of now, everything’s broken while I try to sort out the ambiguities in the header tests.

39

Indeed, pow(var, int) is also ambiguous for me with clang++:

  auto fx = pow(x, 2);
            ^~~
./stan/math/rev/core/std_complex.hpp:2599:26: note: candidate function
inline std::complex<var> pow(const var& lhs, const std::complex<double>& rhs) {
                         ^
./stan/math/rev/core/std_complex.hpp:2616:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, int rhs) {
                         ^
./stan/math/rev/scal/fun/pow.hpp:124:12: note: candidate function
inline var pow(const var& base, double exponent) {
           ^
./stan/math/rev/core/std_complex.hpp:2620:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, double rhs) {
                         ^
./stan/math/rev/scal/fun/pow.hpp:108:12: note: candidate function
inline var pow(const var& base, const var& exponent) {
           ^
./stan/math/rev/core/std_complex.hpp:2624:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs, const var& rhs) {
                         ^
./stan/math/rev/core/std_complex.hpp:2628:26: note: candidate function
inline std::complex<var> pow(const std::complex<var>& lhs,

Looks like I have some template wrangling to do to fix the calls. This one can be fixed by defining pow(var, int).

40

I’ve pushed fixes so that the header tests now pass and pow(int, var) and pow(var, int) work.

The unit tests are still failing to compile in various places due to ambiguities.