Add Vioreanu-Rokhlin quadrature finding

[inducer]

To do:

• Get it to work
• Add a real test for the node finder (Lebesgue constant maybe?)
• Improve basis orthogonalization (@inducer thinks this can't be improved while still representing functions as linear combinations of monomials. Using PKDOs on a simplex avoids this. Unverified conjecture: an orthonormal basis on a mismatched domain will still make this better, e.g. tensor product Chebyshev on a simplex.)
• Add Gauss-Newton finder
• Add proper documentation

Closes #106.

cc @ShawnL00 @a-alveyblanc

[inducer]

It made its first nodes:

🎉

cc @alexfikl

[inducer]

One possible issue right now is the basis orthogonalization. V and R propose an SVD-based procedure that (IIUC?) seems to produce $\ell^2$ (rather than $L^2$)-orthogonal functions. I went with a Cholesky-based approach instead (somewhat analogous to CholeskyQR2, maybe), but it doesn't maintain orthogonality very well for larger (?) bases. For the order-7 example above, it only gets less than 12 digits.

[alexfikl]

I went with a Cholesky-based approach instead (somewhat analogous to CholeskyQR2, maybe), but it doesn't maintain orthogonality very well for larger (?) bases. For the order-7 example above, it only gets less than 12 digits.

Would a pivoted version work better? Although those aren't exactly available easily https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lapack.dpstrf.html#scipy.linalg.lapack.dpstrf

[inducer]

Would a pivoted version work better?

TIL that's even a thing! In CS450, I say things like "Cholesky doesn't need pivoting". And it's not just me:

Now, suppose that the Cholesky decomposition is applicable. As mentioned above, the algorithm will be twice as fast. Furthermore, no pivoting is necessary[...]

For now, I would bet that perhaps the bigger problem is that the implementation (just the test really) currently tries to go where it needs to by linearly combining monomials.

Another indication that it's perhaps not the linear algebra that is to blame here is that, at least by back-of-the-envelope analysis, the main "accuracy loser" in CholeskyQR2 is the formation of $A^TA$, which is not part of the method.

Yet another indication is that the coefficients reach a magnitude of about 3000, which might help explain the loss of about three digits.

[inducer]

This is ready for a look. @ShawnL00, while I can't formally request a review from you, I'd like you to review this code as well.

One thing that's missing is better support for (inefficiently) computing reference integrals of "obstinate" integrands, probably using adaptive quadrature. Some useful pieces exist for that, but I figured I would still leave it for future work.

[ShawnL00]

This is ready for a look. @ShawnL00, while I can't formally request a review from you, I'd like you to review this code as well.

One thing that's missing is better support for (inefficiently) computing reference integrals of "obstinate" integrands, probably using adaptive quadrature. Some useful pieces exist for that, but I figured I would still leave it for future work.

That is cool! will do!

[alexfikl]

Left a few nitpicks, but it all looks pretty straightforward. Very cool! 🎉

Some additional questions:

• How do the nodes from guess_nodes_vr compare to the stored VR nodes? Do they also miss some digits?
• Does orthogonalize_basis work to iteratively construct a basis too? I recall there was a section in Vioreanu2011 doing something like "start with f_i = 1, multiply it by x, orthogonalize, repeat".
• How slow is this to find higher order rules?

[inducer]

• How do the nodes from guess_nodes_vr compare to the stored VR nodes? Do they also miss some digits?

Not a meaningful comparison IMO. guess_nodes_vr gives interpolatory nodes. Off the bat, the nodes out of guess_nodes_vr only integrate $P^N$, nothing more, at least in the test code. The stored VR nodes are the result of additional "tightening" via Gauss-Newton. And our Gauss-Newton widget isn't yet full-featured enough to compare. For instance, if I saw this correctly, the nodes lose symmetry pretty much right away. It does manage to turn nodes for $P^7$ (36 quadrature nodes) into ones that integrate $P^{11}$ (78 functions $\approx 36\cdot2.166$. By way of conjecture, 108=36*3 integrands should be attainable.).

• Does orthogonalize_basis work to iteratively construct a basis too? I recall there was a section in Vioreanu2011 doing something like "start with f_i = 1, multiply it by x, orthogonalize, repeat".

TBH, I did not fully understand that bit of the paper. As it stands, it seems that this process would provide functions that are orthogonal in terms of the $\ell^2$ norm applied to the nodal values, not the $L^2$ sense that seems mathematically needed. It also seems require representing the integrands as a bunch of nodal values, which doesn't leave an obvious way to find basis derivatives that are needed for the Gauss-Newton process. So, in all likelihood, I'm still misunderstanding something.

• How slow is this to find higher order rules?

The tests run within a few seconds, but they only exercise up to $P^{11}$ on the triangle. So it could probably be sped up, but I'll wait with that until it's actually too slow in a specific use case. I haven't profiled, but I suspect the evaluation of the PKDOs is the bottleneck.

[inducer]

Pushed the changes addressing the concerns.

[alexfikl]

Left a few nitpicks in the docs to make them sound more sentence-y, but otherwise this looks good to me! ❤️ (As usual, not sure if these follow line lengths and stuff :(

Not a meaningful comparison IMO. guess_nodes_vr gives interpolatory nodes.

Yeah, forgot that does a lot of extra work to get the nodes. Would be cool to be able to reproduce those one day 😁

TBH, I did not fully understand that bit of the paper.

Yeah, it seemed more like they just added it to the paper to fill in the hole "what if you have a domain where you don't know the basis??".

The tests run within a few seconds, but they only exercise up to P 11 on the triangle.

That's not bad already. I assume that we'll precompute and cache the rules to some order (like the current ones).

[inducer]

@alexfikl Thanks for taking a look!

I assume that we'll precompute and cache the rules to some order (like the current ones).

I would be super excited if it were plausible to run rule generation at code startup and stuff the resulting data somewhere in a cache (instead of plastering a bunch of floating point numbers into source code files). This would of course delay cold-cache runs, but on the upside it would make sure to preserve the recipe for making the rule (as opposed to a bunch of opaque data).

[inducer]

Forgot the most important part: Advertise the existence of this in the README! 😁 3dbe7d3