Engine prototype #13

Vectornaut · 2024-08-20T09:09:14Z

Vectornaut commented

2024-08-20 09:09:14 +00:00

The branch to be merged tracks our work on prototyping a geometric constraint-solving engine. This work is in the engine-proto folder. Most of the code is in Julia (version 1.10.0).

Approaches

We tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising.

Algebraic (In the alg-test subfolder). Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques.
- Gröbner bases (Engine.Algebraic.jl). Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract.
- Homotopy continuation (Engine.Numerical.jl). Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice.
A few notes about our experiences can be found on the engine prototype wiki page.
Gram matrix (in the gram-test subfolder). A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques.
- LDL decomposition (gram-test.sage, gram-test.jl, overlap-test.jl). Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in overlap-test.jl.
- Low-rank factorization (source files listed in findings section). Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (low-rank-test.jl) and an visualization of the loss function landscape near global minima (basin-shapes.jl).
The Gram matrix parameterization wiki page contains detailed notes on this approach.

Findings

With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world.

With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (Engine.jl) and used it to solve five example problems:

overlapping-pyramids.jl
circles-in-triangle.jl
sphere-in-tetrahedron.jl
tetrahedron-radius-ratio.jl
irisawa-hexlet.jl

We plan to use low-rank factorization of the Gram matrix in our first app prototype.

Visualizations

We used the visualizer in the ganja-test folder to visually check our low-rank factorization results. The visualizer runs Ganja.js in an Electron app, made with Blink. Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production.

It seems to have precision problems with low-curvature spheres.
We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details.

The branch to be merged tracks our work on prototyping a geometric constraint-solving engine. This work is in the `engine-proto` folder. Most of the code is in Julia (version 1.10.0). ### Approaches We tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising. - **Algebraic** *(In the `alg-test` subfolder).* Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques. - **Gröbner bases** *(`Engine.Algebraic.jl`).* Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract. - **Homotopy continuation** *(`Engine.Numerical.jl`).* Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice. A few notes about our experiences can be found on the [engine prototype](wiki/Engine-prototype) wiki page. - **Gram matrix** *(in the `gram-test` subfolder).* A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques. - **LDL decomposition** *(`gram-test.sage`, `gram-test.jl`, `overlap-test.jl`).* Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in `overlap-test.jl`. - **Low-rank factorization** *(source files listed in findings section).* Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (`low-rank-test.jl`) and an visualization of the loss function landscape near global minima (`basin-shapes.jl`). The [Gram matrix parameterization](wiki/Gram-matrix-parameterization) wiki page contains detailed notes on this approach. ### Findings With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world. With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (`Engine.jl`) and used it to solve five example problems: - `overlapping-pyramids.jl` - `circles-in-triangle.jl` - `sphere-in-tetrahedron.jl` - `tetrahedron-radius-ratio.jl` - `irisawa-hexlet.jl` We plan to use low-rank factorization of the Gram matrix in our first app prototype. ### Visualizations We used the visualizer in the `ganja-test` folder to visually check our low-rank factorization results. The visualizer runs [Ganja.js](https://enkimute.github.io/ganja.js/) in an Electron app, made with [Blink](https://github.com/JuliaGizmos/Blink.jl). Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production. - It seems to have precision problems with low-curvature spheres. - We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details.

Vectornaut added 113 commits 2024-08-20 09:09:16 +00:00

Start drafting engine prototype b864cf7866

Realize geometric elements as symbolic vectors 4d5aa3b327

Realize relations as equations 463a3b21e1

Order variables by coordinate and then element 86dbd9ea45

In other words, order coordinates like
  (rₛ₁, rₛ₂, sₛ₁, sₛ₂, xₛ₁, xₛ₂, xₚ₃, yₛ₁, yₛ₂, yₚ₃, zₛ₁, zₛ₂, zₚ₃)
instead of like
  (rₛ₁, sₛ₁, xₛ₁, yₛ₁, zₛ₁, rₛ₂, sₛ₂, xₛ₂, yₛ₂, zₛ₂, xₚ₃, yₚ₃, zₚ₃).

In the test cases, this really cuts down the size of the Gröbner basis.

Write a simple solver for the hitting set problem c29000d912

I think we need this to find the dimension of the solution variety.

Correct Minkowski product; build chain of three spheres 59a527af43

Correct relation equations 0731c7aac1

Extend AbstractAlgebra ideals to rational coefficients 6349f298ae

The extension should also let us work over finite fields of prime order,
although we don't need to do that.

Find dimension of solution variety 4e02ee16fc

Use module names as filenames 65d23fb667

You're right: this naming convention seems to be standard for Julia
modules now.

Order spheres before points within each coordinate block a3f3f6a31b

In the cases I've tried so far, this leads to substantially smaller
Gröbner bases.

Switch element abbreviation from "elem" to "elt" 21f09c4a4d

Add relations to center and orient the construction 43cbf8a3a0

Seek sample solutions by cutting with a hyperplane 45aaaafc8f

The example hyperplane yields a single solution, with multiplicity six. You can
find it analytically by hand, and homotopy continuation finds it numerically.

Try a cut that goes through the trivial solution f97090c997

The previous cut was supposed to do this, but I was missing some parentheses.

Show explicitly that all coefficients are 1 in first cut equation 95c0ff14b2

Find witnesses on random rational hyperplanes 34358a8728

Choose hyperplanes that go through the trivial solution.

Try switching to compiled system becefe0c47

Say how many sample solutions we found 06872a04af

Systematically try out different cut planes 8e33987f59

Test a scale constraint af1d31f6e6

In all but a few cases (for example, a single point on a plane), we
should be able to us the radius-coradius boost symmetry to make the
average co-radius—representing the "overall scale"—roughly one.

Separate the algebraic and numerical parts of the engine b3b7c2026d

Try uniformly distributed hyperplane orientations 621c4c5776

Unit normals are uniformly distributed over the sphere.

Test lots of uniformly distributed hyperplanes 6f18d4efcc

Move random cut routine into engine 1f173708eb

Evaluate and display elements 6cf07dc6a1

Try displaying a chain of spheres a450f701fb

For three mutually tangent spheres, I couldn't find real solutions.

Play with two points on two spheres 31d5e7e864

Guess conditions that make the scaling constraint impossible to satisfy.

Explore the performance wall e41bcc7e13

Three points on two spheres is too much.

Study mutually tangent spheres with two fixed 291d5c8ff6

Store elements in arrays to keep order stable 8d8bc9162c

This seems to restore reproducibility.

Make results reproducible ae5db0f9ea

Find real solutions for three mutually tangent spheres ba365174d3

I'm not sure why the solver wasn't working before. It might've been just
an unlucky random number draw.

Test different sign patterns for cosines f2000e5731

It seems like there are real solutions if and only if the product of the
cosines is positive.

Clean up example of three mutually tangent spheres 3170a933e4

Try out the Gram matrix approach 16826cf07c

Extend Gram matrix automatically 717e5a6200

The signature of the Minkowski form on the subspace spanned by the Gram
matrix should tell us what the big Gram matrix has to look like

Correct signature ef33b8ee10

Try numerical low-rank factorization 58a5c38e62

The best technique I've found so far is the homemade gradient descent
routine in `descent-test.jl`.

Get a Ganja.js visualization running in Blink d1ce91d2aa

Get familiar with Ganja.js inline syntax 3c34481519

Clean up examples 3b10c95d5f

Declare JavaScript variables. Revise Julia comments to match new code.

Add live updates to Ganja.js visualization 69a9baa8ee

Play with reflections 06a9dda5bb

Try configuration of five tangent spheres.

Load elements from Julia into Ganja.js b7b5b9386b

Generate palette automatically 182b5bb9f6

Correct indentation of CSS 665cb30ce0

Build construction viewer module a3b1f4920c

Update namespace of AbstractAlgebra.Rationals 4a28a47520

Merge branch 'ganja' into gram 5aadfecf6c

Visualize low-rank factorization results.

Avoid naming conflict with identity transformation 2b6c4f4720

Switch to Ganja.js basis ordering c933e07312

Size the viewer window automatically 7aaf134a36

Add element visibility controls 3eb4fc6c91

Streamline visibility controls 5ea32ac53c

Merge branch 'ganja' into gram a113f33635

Get visibility controls.

Let visibility controls scroll 05a824834d

Merge branch 'ganja' into gram 8eb1ebb8d2

Get Ganja.js to display planes 242d630cc6

Do gradient descent entirely in BigFloat e7dde5800c

The previos version accidentally returned steps in Float64.

Encapsulate gradient descent code 133519cacb

The completed gram matrix from this commit matches the one from commit
e7dde58 to six decimal places.

Name gradient descent test more specifically 17fefff61e

Sketch random vector generator abc53b4705

This needs to be rewritten: it can fail at generating spacelike vectors.

Improve random vector generator 7e94fef19e

Host Ganja.js locally d39244d308

Randomize guess in gradient descent test ea354b6c2b

Randomly perturb the pre-solved part of the guess, and randomly choose
the unsolved part.

Test gradient descent for sphere in tetrahedron 736ac50b07

Add sphere and plane utilities to engine 828498b3de

Test gradient descent for circles in triangle 9efa99e8be

Add required package to "sphere in tetrahedron" example 93dd05c317

Track slope in gradient descent history 610fc451f0

Start "circles in triangle" from a very close guess 023759a267

Change loss function to match gradient 77bc124170

Visualize neighborhoods of global minima f84d475580

Require triangle sides to be planar 5652719642

Record gradient and last line search in history 4d5ea062a3

Correct improvement threshold by using unit step d538cbf716

Our formula for the improvement theshold works when the step size is
an absolute distance. However, in commit `4d5ea06`, the step size was
measured relative to the current gradient instead. This commit scales
the base step to unit length, so now the step size really is an absolute
distance.

Use Newton's method for polishing 3910b9f740

Sketch backtracking Newton's method 25b09ebf92

This code is a mess, but I'm committing it to record a working state
before I start trying to clean up.

Clean up backtracking gradient descent code 7b3efbc385

Drop experimental singularity handling strategies. Reduce the default
tolerance to within 64-bit floating point precision. Report success.

Preserve explicit zeros in Gram matrix conversion 53d8c38047

In previous commits, the `circles-in-triangle` example converged much
more slowly in BigFloat precision than in Float64 precision. This
turned out to be a sign of a bug in the Float64 computation: converting
the Gram matrix using `Float64.()` dropped the explicit zeros, removing
many constraints and making the problem much easier to solve. This
commit corrects the Gram matrix conversion. The Float64 search now
solves the same problem as the BigFloat search, with comparable
performance.

Switch back to BigFloat precision in examples 94e0d321d5

Switch to backtracking Newton's method in Optim b185fd4b83

This performs much better than the trust region Newton's method for the
actual `circles-in-triangle` problem. (The trust region method performs
better for the simplified problem produced by the conversion bug.)

Implement frozen variables 1ce609836b

Don't constrain self-product of frozen vector 7c77481f5e

Make tetrahedron faces planar e6cf08a9b3

Switch engine to light cone basis bde42ebac0

Write examples directly in light cone basis 2038103d80

Give spheres positive radii in examples 4728959ae0

This changes the meaning of `indep_val` in the overlapping pyramids
example, so we adjust `indep_val` to get a nice-looking construction.

Start tetrahedron radius ratio example ea640f4861

Add the vertices of the tetrahedron to the `sphere-in-tetrahedron`
example.

Revert "Give spheres positive radii in examples" 5abd4ca6e1

This reverts commit 4728959ae0, which
actually gave the spheres negative radii! I got confused by the sign
convention differences between the notes and the engine.

Tetrahedron radius ratio: correct signs 6d233b5ee9

Correct point utility d51d43f481

Tetrahedron radius ratio: correct vertex guesses 6e719f9943

Tetrahedron radius ratio: add circumscribed sphere a02b76544a

Tetrahedron radius ratio: tweak guess 96ffc59642

Jiggle the vertex guesses. Put the circumscribed sphere guess on-shell.

Tetrahedron radius ratio: find radius ratio 01f44324c1

Use notes' sign convention for light cone basis 69a704d414

Correct point utility again d0340c0b65

The balance between the light cone basis vectors was wrong, throwing the
point's coordinates off by a factor of two.

Correct sign of normal in plane utility 74c7f64b0c

Clarify the relevant notes too.

Clarify notes on tangency 24dae6807b

Clean up formatting of notes 3764fde2f6

Circles in triangle: correct frozen variables a7f9545a37

Since the self-product of the point at infinity is left unspecified, the
first three components can vary without violating any constraints. To
keep the point at infinity where it's supposed to be, we freeze all of
its components.

Circles in triangle: jiggle the guess 9007c8bc7c

Drop old code from examples b040bbb7fe

Report success correctly when step limit is reached b24dcc9af8

Correct scope of guess constants 33c09917d0

Overlapping pyramids: drop outdated comment 71c10adbdd

Clean up example formatting 19a4d49497

Add Irisawa hexlet example a26f1e3927

Hat tip Romy, who sent me the article on sangaku that led me to this
problem.

Irisawa hexlet: drop unviable approach 8a77cd7484

The approach in the deleted file can't work, because the "sun" and
"moon" spheres can't be placed arbitrarily.

Irisawa hexlet: use Abe's terminology in comments 9d69a900e2

Abe uses the names "sun" and "moon" for what Wikipedia calls the nucleus
spheres.

Stow algebraic engine prototype d7dbee4c05

We're using the Gram matrix engine for the next stage of development,
so the algebraic engine shouldn't be at the top level anymore.

glen added 1 commit 2024-09-19 02:50:49 +00:00

doc: Slight mods to the inversive notes 8084fdeab0

glen added 1 commit 2024-09-19 02:52:16 +00:00

doc: typo in inversive notes bd3e3506e5

glen added 1 commit 2024-09-19 03:12:29 +00:00

doc: More elaboration of plane coordinates in inversive notes a182b66301

glen added 1 commit 2024-09-19 03:27:21 +00:00

doc: Another note about inversive coordinates 23ecca3963

glen commented

2024-09-19 03:39:35 +00:00

OK, by way of starting the review of this PR, I looked over the first six rows of notes/inversive.md, and had a few questions:

I made a few changes, please double-check me, thanks.
The right column of the third row had become garbled. Is there in fact a scaling choice in the point coordinates? I think so, I think they are like projective coordinates where we could represent (1,4,8) by (9,1,1,4,8) as described, but just as well by (18,2,2,4,8), etc. We simply pick the one with second coordinate 1 for canonicalness. Please fill in the ... as you feel appropriate.
In the sixth row, you give a nice example about orientations matching. It seems to beg the question though about whether our sphere coordinates are similarly oriented, or whether in these coordinates they always have outward-pointing normals. Is there another representation of the unit sphere centered at (0,0,1) but with inward-pointing normals? Is that represented by negative radius, i.e. (0, -1, 0, 0, 1)? That vector has Q-norm -1, but sadly it still has Q-product -1 with the plane (0,0,0,0,1), which would seem to refute that now their normals are opposite at the point of tangency. So does (0, -1, 0, 0, 1) represent anything in our inversive coordinates, and if so what? Please put appropriate notes in the first few rows of these notes to clarify these points, thanks.

I will finish the review when you've had a chance to do these things, thanks.

OK, by way of starting the review of this PR, I looked over the first six rows of notes/inversive.md, and had a few questions: 1) I made a few changes, please double-check me, thanks. 2) The right column of the third row had become garbled. Is there in fact a scaling choice in the point coordinates? I think so, I think they are like projective coordinates where we could represent (1,4,8) by (9,1,1,4,8) as described, but just as well by (18,2,2,4,8), etc. We simply pick the one with second coordinate 1 for canonicalness. Please fill in the ... as you feel appropriate. 3) In the sixth row, you give a nice example about orientations matching. It seems to beg the question though about whether our sphere coordinates are similarly oriented, or whether in these coordinates they always have outward-pointing normals. Is there another representation of the unit sphere centered at (0,0,1) but with inward-pointing normals? Is that represented by negative radius, i.e. (0, -1, 0, 0, 1)? That vector has Q-norm -1, but sadly it still has Q-product -1 with the plane (0,0,0,0,1), which would seem to refute that now their normals are opposite at the point of tangency. So does (0, -1, 0, 0, 1) represent anything in our inversive coordinates, and if so what? Please put appropriate notes in the first few rows of these notes to clarify these points, thanks. I will finish the review when you've had a chance to do these things, thanks.

glen commented

2024-09-20 02:48:06 +00:00

P.S. I have now understood from today's discussion that the coordinates in the current prototype are not quite exactly these. I think we should probably leave these notes alone, but just clearly document wherever the coordinate entities are declared exactly what they are are and how they relate to these coordinates. Thanks!

This pull request can be merged automatically.

You are not authorized to merge this pull request.

View command line instructions.