Complex number support in Math library autodiff


I am yet to be convinced there is sufficient need to expose complex numbers. Maybe some exposed functions can return a tuple with real and imaginary parts (perhaps not in the transformed parameters or model blocks), but I think most useful things can be implemented internally in Stan Math.


I am in agreement - unless someone implements a full complex algebra everything will be piecemeal implementations that are more appropriately wrapped in the function for which the partial complex algebra implementation is needed. In that case we should be able to get away with functions that return real outputs (i.e. power and phase functions from an FFT, magnitude and phase functions for input complex numbers, etc) instead of a partial complex algebra implementation that is going to be extraordinarily confusing to users.


On the original question, when I asked for this, I just wanted the largest eigenvalue for certain nonsymmetric matrices. This was actually for cases where you can guarantee they are real, but it seems you cannot compute them without complex numbers.


(post withdrawn by author, will be automatically deleted in 24 hours unless flagged)


I would like to request a code review of the associated PR

@seantalts @Bob_Carpenter ??
Also, I would like to use Sean’s option mentioned in the PR about taking the tests as the spec.


:))) Travis build is ALL GREEN now (because the apt-repo was available for it to install gcc this time, yay).
Trying one more build for Jenkins’ Doxygen thing


Thanks for putting this together. This is definitely going to go beyond my analysis competence, so someone else should review it.

Spec for what?

The summary for the PR should be more specific as to what “initial attempt” means and what the scope of the PR is. The cited issue is just a bug report about a particular function. Is the PR to fix that function or is it more general? Will the behavior of other functions change?

We don’t want to start merging works in progress until they provide useful functionality.


(post withdrawn by author, will be automatically deleted in 24 hours unless flagged)


I think what people were saying about complex derivatives in the related issue 123 was complicating things more than they have to be. A multiple derivative could be rewritten in terms of a complex number and its complex conjugate, or it could just as easily (and more sanely) be rewritten in terms of the explicitly real and explicitly imaginary components of the complex number. They are two sides of the same coin. Either way, we still must end up with the derivative with respect to two different components, simply because imaginary numbers always have two components. The std::complex class tracks explicitly real and explicitly imaginary components, so that is what we allow in the stan complex code also. It really is nothing more than tracking var + imaginary I * var, where the vars are treated as real numbers (and the complex class handles the constant I multiplication that makes the 2nd var purely imaginary).

Hrm, I would say it was a lot harder to get it to a state for merging than just a state that would temporarily work against a given development branch. For me, the goal has been to make it into stan-math so that I would be able to rely on the functions delivered for years to come. I feel like I have done that, but will do whatever y’all need to give the go-ahead.

I updated the PR first comment itself to put some answers there. I have also added more below for context.

This is Sean’s comment on specs for context. He said this after I asked for a format for the specs in the PR thread.

We don’t have ay specific format, just something that lets us know what operations you want to work, what’s not supported, that kind of thing. If you’re dead set on coding before getting feedback on a spec, your tests might function as a sort of very detailed and functional spec, which I think would be great.

There are two simple tests, one is for reverse, and one is for forward.
reverse: test/unit/math/rev/core/var_test.cpp
forward: test/unit/math/fwd/scal/fun/abs_test.cpp

At first they set the real component only. This gives the opportunity to test that the null dereference is sidestepped correctly for var.

Then the imaginary component is set to another number that makes it possible to easily check that Re^2 + Im^2 = abs^2 (essentially, 3,4,5 right triangle).

After that, a derivative is propagated through the system, and the comments show the calculations for what the derivative should be.

Lastly, addition, multiplication, and division are performed in left and right fashions to arrive at a number that when 1 * i is subtracted from it, equals zero. This is to help check the relevant operator templates.

Here is the mixed mode test: test/unit/math/mix/mat/functor/autodiff_test.cpp

I chose the more complex mixed mode gradient hessian test because people have specifically complained about Eigen’s desire to return eigen-decompositions of matrices, including ones of stan’s vars and fvars, as std::complex data types, which were not previously implemented.

This allows us to illustrate that the complex code could run the gauntlet of hellacious operations and still satisfy the eigendecomposition identities that are being checked on each iteration of newton’s method.

Newton’s method is used to drive one of the returned eigenvalues to a critical point where its real and imaginary components are equal (via a merit function that compares them).

Edit: Rotation matrix info: See the last two math pictures and the last 3 paragraphs of this section


This is a Mathematica notebook that I used to help myself derive the solution for the mixed mode derivative test.