Announcing `PolyChaos.jl`

Package Announcements
1

I am happy to announce PolyChaos.jl, a collection of functions for orthogonal polynomials, quadrature rules, and polynomial chaos expansions.

What’s the idea?

Let’s say you have some non-negative function w: \Omega \longrightarrow \mathbb{R} and you would like to find the polynomals that are orthogonal relative to w. That’s the core functionality of PolyChaos.jl.

Why is that useful? Once you know the orthogonal polynomials, you can compute quadrature rules such aus Gauss, Gauss-Radau, Gauss-Lobatto. These quadrature rules allow you to solve integrals

\int_\Omega f(x) w(x) \mathrm{d} x \approx \sum_{i=1}^N \omega_i f(\tau_i).

The numbers (\omega_i, \tau_i) for i = 1, \dots, N are the quadrature rule that PolyChaos.jl computes.

Finally, orthgonal polynomials are intricately linked to random variables. For example, Hermite polynomials (more specifically, probabilists’ Hermite polynomials) happen to be the orthogonal polynomials relative to the probability distribution of a standard Gaussian random variable, the well-known density w(x) = \exp(-x^2/2)/\sqrt{2 \pi}. The mathematical method that investigates the relation between orthogonal polynomials and random variables is called polynomial chaos expansion, dating back to the infamous Norbert Wiener. Polynomial chaos is to random variables what Fourier series is to periodic signals – it allows you to represent some complicated mathematical object by a few deterministic numbers. Too abstract? Think of a Gaussian random variable – everything you need to know about it are its mean and variance. PolyChaos.jl let’s you compute polynomial chaos expansions of random variables.

I am fully aware that there are several great packages for quadrature out there such as QuadGK, FastGaussQuadrature.jl or ApproxFun.jl. As I see it, PolyChaos.jl is a complement, and by no means a replacement for those excellent packages. If you are interested in how these quadrature rules come about, then perhaps PolyChaos.jl can provide insight(s).

The idea for the package was born out of the need to have a Julia package for polynomial chaos expansions (there was also a poster at JuliaCon 2018). There is still a lot that can be added (and I will create issues for that) so any form of contribution is welcome!

This being my first software project I am sure there will be unexpected stuff to happen. Regardless, it has been a great experience so far and I am actively using the package for my research. Let’s see what happens!

23 Likes
2

Nice! This’ll be a good starting place for SPDE expansions!

3

Yes!

A similar to-do is this: auto-generate Galerkin projections of ODEs. In the documentation of PolyChaos I did it by hand for a scalar linear ODE.

4

Awesome stuff and great tutorials too. I’ve never used these types of approaches but I’ve looked into them in the past. Your documentation gives one of the most approachable introductions to the topic that I’ve seen.

5

Good news: PolyChaos.jl v0.1.1 is out. The improvements as summarized in NEWS.md:

  • improved performance of recurrence coefficients, quadrature rules, evaluation of orthogonal polynomials
  • added considerable number of tests for recurrence coefficients and quadrature rules
  • added code coverage
  • added GNU General Public License v3.0
  • removed dependency on IterTools
6

Is there any particular reason you chose GPL? In the Julia community it is more usual to use MIT.

7

The GaussQuadrature.jl package can generate all of the classical Gauss rules in multiple precision (as BigFloat arrays), and also rules for a log weight function.

2 Likes
8

It’s mostly because my bosses said so (great reason, I know…). I am aware of the discussions about the licenses and I am doing by best to persuade them. :slight_smile:

9

This is great, thanks! I have to admit that our focus was not so much on arbitrary precision of quadrature rules, but having a workflow from weight → recursion coefficients → quadrature rule → polynomial chaos. It just so happened that quadrature rules became an important part. PolyChaos also has support for Lobatto and Radau quadrature.

There is a difference though: I saw that you implemented the Chebyshev and modified Chebyshev algorithms. In PolyChaos I chose to go with Stieltjes and Lanczos procedures. The reason being that the Chebyshev algorithm can lead to ill-conditioned problems. There is also an implementation of Gautschi’s multiple discretization procedure which allows you do to this for example.
Having said all that, we should perhaps get in touch and chat! There is interesting overlap.

10

(offtopic) Another discussion about licenses: Bitcoin.jl released as GPL, later switched to MIT.

11

My GaussQuadrature.jl package is mostly a re-write of a Fortran 95 code descended from the original implementation gaussq.f of the Golub-Welsch algorithm.

12

It would be interesting the replace the current golubwelsch() with your special_eigenproblem!().

13

I am pleased to announce that PolyChaos.jl is now available in version 0.2.0. The major difference to the 0.1.x versions is a cleaner type hierarchy, building on top of AbstractMeasure, AbstractOrthoPoly, AbstractQuad, and AbstractTensor. Check out the docs for details.

There are also new examples: orthogonal bases for Gaussian mixture models, and polynomial chaos for optimal power flow problems. Other power flow problems are implemented here.

2 Likes