I’m wrapping up support for complex number types in the Stan math library. Now I’m wondering what we can do with complex numbers in statistical models.

**Functions operating in the complex domain**

The initial plan is to add some matrix functions that use complex numbers internally:

- fast fourier transforms
- asymmetric eigendecomposition
- Schur decomposition

The eigendecomposition and Schur decomposition are already working. We need to settle on an interface for FFTs before adding those. Are there other functions like these we should add?

**Complex random variables?**

But what if we add a complex number type to the language? That’d essentially give us complex random variables (and constants) as first-class citizens. We’d have something like the following.

complex x; // declaration x = complex(-1, 2.9); // constructor x.real = 1.3; // setters x.imag = -2.5; complex y = sqrt(x); // complex sqrt real yr = y.real; // getters real yi = y.imag;

Technically, we’re treating a complex number just like a pair (the elements being the real and imaginary components). We will be able to differentiate with respect to real and/or imaginary components, and thus differentiate through functions operating in the complex domain, but we’re not taking derivatives with respect to an entire complex number.

It’s easy in the math library to create a matrix of complex numbers, but we’ll need new types in the language: `complex_matrix`, `complex_vector`, `complex_row_vector`. Arrays will just follow from general principles.

**Applications, anyone?**

So what can we do with complex random variables if we have them?

A quick search online found a discussion of complex variables on stats.stackexchange which mirrors some of our internal discussions.

I have no idea where we’d use something like characteristic functions in the Stan language or complex polynomial roots, but am hoping someone else can clue me in.

This falls under the rubric of eigendecomposition, but the general tests for stationarity in time series models is to arrange the coefficients in a square matrix in a particular way, get the eigenvalues, and then check that all of them have an absolute value below 1. If you try to do this in the generated quantities block in Stan currently, my recollection is that you can run into problems when the eigenvalues are complex. It’s safer to just do it outside of Stan. So having a complex type in Stan is an obvious use case here.

It’s a much harder problem than FFTs, but other similar transforms, in particular the spherical harmonic transform (scalar and tensor), would be incredible useful in astronomy (and other fields, I imagine). But various issues of precision alongside the related specific ordering in which it is most efficient to actually compute the harmonics, make this non-trivial. See, for example, https://sourceforge.net/projects/libsharp and https://healpix.sourceforge.io/documentation.php.

(The hard part of this stuff isn’t really the complex number aspect…)

We can do quantum probability! (relevant for recent discussion)

Thanks, I was thinking I totally misunderstood that discussion.

There was some earlier work on fitting macroeconomic models in Stan https://pdfs.semanticscholar.org/0a9b/35f9f11ce0c489a911622d700d2a4e64f385.pdf which faced a challenge because the likelihood for these models typically involves a *generalized* Schur decomposition, which jointly decomposes a matrix pair (A,B) into an upper triangular pair in a way similar to how the Schur decomposition generates an upper triangular matrix from a single matrix. There are a few research groups I am aware of (including myself, for what it’s worth) for whom the absence of this functionality is the primary reason for not using Stan. We are currently coding up some alternatives in Julia and running AdvancedHMC.jl for HMC sampling, but full Stan functionality would save a lot of trouble. There is a real Schur decomposition as well which avoids complex numbers by producing a block triangular representation and suffices for these purposes, so you might not strictly need complex numbers to do this, but if you’re taking requests on functions to add…

Interestingly, the reason this is needed is precisely for the reason previous commenter John Hall mentions: the models impose time series stationarity in a way which requires checking the eigenvalues. But because this is done inside of the model, since you can’t even define the likelihood uniquely without pinning down that you want the stationary solution, it can’t be swept into post-processing outside of Stan.

To clarify, I have not heard of anyone asking for the Schur decomposition of a complex matrix in Eigen (although it is there):

http://eigen.tuxfamily.org/dox/classEigen_1_1ComplexSchur.html

Rather, the Schur decomposition of a real matrix in Eigen utilizes some std::complex arithmetic internally that would have caused segfaults had we tried to use it before. With a real Schur decomposition, we can do more of the Dynamic Stochastic General Equilibrium (DSGE) models.

Also, forgot, although this is just an application of FFT, we can do the Whittle [approximation](https://watermark.silverchair.com/asy071.pdf?token=AQECAHi208BE49Ooan9kkhW_Ercy7Dm3ZL_9Cf3qfKAc485ysgAAAmEwggJdBgkqhkiG9w0BBwagggJOMIICSgIBADCCAkMGCSqGSIb3DQEHATAeBglghkgBZQMEAS4wEQQMcx8C_b3PES5m_gYLAgEQgIICFLDZIYJ8j2PkfEVkOh69VM8KCPUuZOBq8VbzJiaozSjbABsfGS8d-QOFwaG_zHIBu12t7LMmvU4Wve_zyHXWEOQeHGlXyizNDLwhqebXArHgJRz0A_p8loPqRsmpD7uG19UsxicGgGBlQxtpqYboOjRyBJdQL4SX7sEb6yqkfadEWxl2x62LzfjE4GczMqQjMAVio87KsrIKo9YNrEWthCpVAWHkTNkIBZj61Ko_u1c6Hb5VpOy-OPKAayVEiW1MCsKnweJeDkS63xogc0d2QEOPR5tBoeUwp0E2e7efb-UlvEypodlbdFlzzUEWAcKrS30jhrWefIq5N_qbA1-tBxrFPAnxa_hZqdRQeAjyIAn4sbV5VTafH0uGPmCkzIbgAholz8N3jkR9AYH5tMyA0AjLGHg4OL43GppdjtV-8Tw7G4Luvgz9LuTDCYueOYFXKcWYBA88v71Gjo4J0IPBvMPh37_f8Nv-HDlmywq3sxJEBk2gSj9isPz3MV4ULdZuoGhLJw9H_j98KMj1QMvAJPCDNYkwrulDgsLIrDgpmaQ5NFQEnj97WXhKvkvObRVfu-cr9q6H3PXtwePD_iVvnH8lmLROovxFJrWBvzMc1Xt9hFZrQcxfZUrhxWLezx1mgqO1kLFHzn70Wam3MRNkd1FVOblXtROGeJgxrqnN3rIDY0bhSy22-BtUyizSOxZOz-W5Nt8) to big multivariate normal distributions

Dead link. I get the message “Your session has timed out. Please go back to the article page and click the PDF link again.”

Try https://discovery.ucl.ac.uk/id/eprint/10069159/8/Olhede_asy071.pdf or Google for “The debiased Whittle likelihood”.

I’ve used complex numbers for computing the CDF or density of a convolution, as in this paper (https://www.jstor.org/stable/2334555) which looks at linear combinations of chi-squareds, the distribution of a quadratic form in Normals. I think I’ve seen the same approach used for the CDF of the sum of bernoulli variables with different means.

The probabilistic formulation of non-negative matrix factorisation with the IS divergence (IS-NMF) uses a complex normal distribution:

https://www.irit.fr/~Cedric.Fevotte/publications/proceedings/eusipco09a.pdf

Although I don’t see this as strong motivation for putting complex numbers in Stan.

Small point: be careful choosing a branch for the complex square root, and document the choice.