Cylindrical Algebraic Decomposition tactic

Stream: general

[Luigi Massacci]

Bolton Bailey said:

To return to the topic of the thread, a more ambitious missing tactic would be something that resolves statements in the first-order theory of the reals. I seem to recall watching a video of a talk someone gave about work formalizing CAD algorithms but I unfortunately am having trouble finding the link.

The thing with CAD is that Collins' algorithm is kinda monstrous. There were two separate attempts (by some experienced users too) to make a Coq tactic but neither came to completion. Also if you consider the theoretical complexity and the fact that it hasn't really been tried extensively in practice (so someone might implement it only to come to the conclusion it's too slow), I doubt that it will get done any time soon

[Luigi Massacci]

At least not by one or two people. Maybe as a bigger team effort? But then the problem is finding volunteers

[Bolton Bailey]

Luigi Massacci said:

Bolton Bailey said:

To return to the topic of the thread, a more ambitious missing tactic would be something that resolves statements in the first-order theory of the reals. I seem to recall watching a video of a talk someone gave about work formalizing CAD algorithms but I unfortunately am having trouble finding the link.

The thing with CAD is that Collins' algorithm is kinda monstrous. There were two separate attempts (by some experienced users too) to make a Coq tactic but neither came to completion. Also if you consider the theoretical complexity and the fact that it hasn't really been tried extensively in practice (so someone might implement it only to come to the conclusion it's too slow), I doubt that it will get done any time soon

Very good points. Although I would point out that, in theory, the problem that polyrith solves is apparently also doubly-exponential. polyrith is still useful though, I presume because:

• Heuristics make it fast in many cases
• Perhaps people are simply using the tactic on very small problem instances and still getting usefulness out of the tactic that way.

[Luigi Massacci]

Yes, that's why I mentioned the fact that (unlike Buchberger's & co) CAD hasn't really been tried extensively, so you don't have that "it sucks in theory but it's not bad in practice" defense. As far as I know of course

[Luigi Massacci]

Maybe getting a tactic for computing Grobner bases or something similar into Lean first might be a good trial for trying something for CAD later

[Mario Carneiro]

Regarding a CAD tactic, I've wanted to embark on that on multiple occasions, but the big missing piece is sufficient theory regarding the algebraic numbers. Last I checked we still don't have a proof that the algebraic numbers are a field

[Mario Carneiro]

I don't think implementation concerns are the blocker

[Johan Commelin]

@Mario Carneiro you mean with computable field operations? Otherwise I'm sure we have it

[Mario Carneiro]

We don't

[Bolton Bailey]

docs#AlgebraicClosure.commRing is the algebraic closure of a field a field in general?

[Mario Carneiro]

the ring axioms are the hard ones anyway

[Mario Carneiro]

it's not clear that this can be used to prove that the algebraic numbers, as a subset of the complex numbers, are closed under addition and multiplication

[Mario Carneiro]

AlgebraicClosure constructs a completely new set with unclear relation to the original

[Bolton Bailey]

I seem to recall that Descartes' rule of signs was a prerequisite to this.

[Mario Carneiro]

last time I looked into this we couldn't find anything in present mathlib

[Johan Commelin]

Can we take the normal closure, or the integral closure? Those are relative constructions so they give a subtype of C

[Johan Commelin]

@Mario Carneiro Are these helpful?

lemma foo (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x + y) := by
  rw [isAlgebraic_iff_isIntegral] at hx hy ⊢
  exact hx.add hy

lemma bar (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x * y) := by
  rw [isAlgebraic_iff_isIntegral] at hx hy ⊢
  exact hx.mul hy

[Mario Carneiro]

yes, what are the imports there?

[Johan Commelin]

def AlgNums := (integralClosure ℚ ℂ).copy {x : ℂ | IsAlgebraic ℚ x} <| by
  ext
  apply isAlgebraic_iff_isIntegral

example : CommRing AlgNums := Subalgebra.toCommRing _

[Johan Commelin]

Ooh, can't shake tell you the imports?

[Mario Carneiro]

I don't know, you didn't write a MWE

[Johan Commelin]

import Mathlib.Data.Complex.Basic
import Mathlib.RingTheory.Algebraic

section test

variable {x y : ℂ}

def AlgNums := (integralClosure ℚ ℂ).copy {x : ℂ | IsAlgebraic ℚ x} <| by
  ext
  apply isAlgebraic_iff_isIntegral

example : CommRing AlgNums := Subalgebra.toCommRing _

lemma foo (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x + y) := by
  rw [isAlgebraic_iff_isIntegral] at hx hy ⊢
  exact hx.add hy

lemma bar (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x * y) := by
  rw [isAlgebraic_iff_isIntegral] at hx hy ⊢
  exact hx.mul hy

end test

[Johan Commelin]

I had import Mathlib before

[Johan Commelin]

Note that AlgNums is even computable

[Mario Carneiro]

Is addition on it computable?

[Mario Carneiro]

well, I guess it is since it's complex addition

[Mario Carneiro]

We do need computable addition and multiplication for this application, so now I'm seeing whether the proof method supports this at all

[Mario Carneiro]

Having the noncomputable proof is the main thing though since you can turn it into a constructive proof fairly mechanically and use pieces from the original proof for any parts that aren't part of the computation

[Mario Carneiro]

It looks like most of the magic is in docs#LinearMap.exists_monic_and_coeff_mem_pow_and_aeval_eq_zero_of_range_le_smul

[Johan Commelin]

import Mathlib.Data.Complex.Basic
import Mathlib.FieldTheory.IntermediateField
import Mathlib.RingTheory.Algebraic

section test

variable {x y : ℂ}

def AlgNums := (integralClosure ℚ ℂ).copy {x : ℂ | IsAlgebraic ℚ x} <| by
  ext
  apply isAlgebraic_iff_isIntegral

example : CommRing AlgNums := Subalgebra.toCommRing _

lemma foo (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x + y) :=
  AlgNums.add_mem hx hy

lemma bar (hx : IsAlgebraic ℚ x) (hy : IsAlgebraic ℚ y) : IsAlgebraic ℚ (x * y) :=
  AlgNums.mul_mem hx hy

lemma baz (hx : IsAlgebraic ℚ x) : IsAlgebraic ℚ (-x) :=
  AlgNums.neg_mem hx

lemma qux (hx : IsAlgebraic ℚ x) : IsAlgebraic ℚ (x⁻¹) := by
  obtain rfl|hx₀ := eq_or_ne x 0
  · simpa only [inv_zero] using isAlgebraic_zero
  rcases hx with ⟨φ, h₀, h⟩
  refine ⟨φ.reverse, ?_, ?_⟩
  · simpa only [ne_eq, Polynomial.reverse_eq_zero] using h₀
  have : Invertible x := hx₀.isUnit.invertible
  rw [← invOf_eq_inv]
  exact (Polynomial.eval₂_reverse_eq_zero_iff _ _ _).mpr h

def AlgebraicNumber : IntermediateField ℚ ℂ :=
  { AlgNums with
    inv_mem' := fun _ hx ↦ qux hx }

noncomputable
instance : Field AlgebraicNumber := inferInstance

end test

[Johan Commelin]

Voila, there you have your field.

[Notification Bot]

33 messages were moved here from #general > "Missing Tactics" list by Mario Carneiro.

[Adam Topaz]

Is the fact that AlgNum is algebraically closed also needed?

[Adam Topaz]

How hard would that be?

[Aaron Anderson]

Cylindrical Algebraic Decomposition is close to my research interests - I’m busy right now, but if someone gets a project going, I’d definitely volunteer

[Johan Commelin]

Shouldn't be hard. If a complex number x is algebraic over AlgNum then by transitivity it is algebraic over Q. Hence it is already an element of AlgNum

[Adam Topaz]

Do we have the appropriate notion of transitivity?

[Johan Commelin]

I hope so!

[Adam Topaz]

I.e. over the field generated by the coefficients in the polynomial equation…

[Adam Topaz]

Or do you have something else in mind?

[Johan Commelin]

https://leanprover-community.github.io/mathlib4_docs/Mathlib/RingTheory/Algebraic.html#Algebra.IsAlgebraic.trans is the best I can find

[Adam Topaz]

I suppose AlgNum could(should?) be defined using docs#integralClosure ?

[Mario Carneiro]

It was, above

[Mario Carneiro]

I did a deep dive with Johan on the proof here to find out how to make it constructive, and we found that the construction essentially uses A.charpoly where A is the matrix which describes how to express $(x+y)x^iy^j$ in terms of the basis elements $x^iy^j$ for i,j less than the degree of their respective minimal polynomials. Unfortunately this seems quite a bit less efficient than https://en.wikipedia.org/wiki/Resultant which is the standard reference for the proof that algebraic numbers are closed under addition, and it's also not obviously related (one is the determinant of a $mn\times mn$ matrix, the other one is the determinant of a completely different $(m+n)\times(m+n)$ matrix)

[Mario Carneiro]

So if we want the resultant construction there is still some mathlib theory work to be done here

[Mario Carneiro]

The resultant has a few interesting properties and also comes up in the definition of the discriminant, so I think there are independent reasons to want to have its basic theory

[Mario Carneiro]

Wikipedia also mentions an n-ary version of the resultant called "Macaulay's resultant" which also might be interesting, although it's not as computationally motivated since it's a ridiculously large polynomial

[Mario Carneiro]

Another thing that is missing for turning this into a practical algorithm is that det uses the O(n!) definition and we don't have a proof of the gauss definition of the determinant yet

[Mario Carneiro]

Do we have LU decomposition, QR decomposition or Cholesky decomposition? Apparently they can all be used for computing determinants

[Kevin Buzzard]

Fields Medallist Alan Baker used to teach the explicit resultant proof that the sum of two algebraic numbers was algebraic in his Cambridge Part III course and we all used to laugh at him for not using the "subspace of a finite-dimensional vector space is finite-dimensional" proof -- he claimed that he didn't like to use "modern concepts" when dealing with the material (this was in the 1980s).

[Geoffrey Irving]

I think the mental model here may be that everything needs to be rewritten to be not just computational but actually fast, so the existence of current theory gives you only so much. You'll also likely want approximate arithmetic to speed things up, as in practice working with exact algebraic numbers is very slow (you need them, but you also need shortcuts).

[Geoffrey Irving]

The interval arithmetic I currently have is one component, though for CAD you almost certainly need the arbitrary precision version. If it seems like that's useful I can do the conversion: the code actually gets easier going from fixed precision to infinite precision due to side conditions disappearing.

[Antoine Chambert-Loir]

Mario Carneiro said:

So if we want the resultant construction there is still some mathlib theory work to be done here

Do resultants allow to prove that the product of two algebraic elements is algebraic ? — I'm not sure…

[Antoine Chambert-Loir]

@Johan Commelin , your code indicates that mathlib is missing IsAlgebraic.add, etc.

[Patrick Massot]

Yes resultants also do products.

[Eric Rodriguez]

The resultant is in flt-reg, iirc

[Eric Rodriguez]

There's been a big need to port it for a while but it's still not been done

[Bolton Bailey]

I finally remembered the talk I saw about CAD, it was James Davenport's talk at the IPAM conference.

[Tomaz Mascarenhas]

Mario Carneiro said:

I did a deep dive with Johan on the proof here to find out how to make it constructive, and we found that the construction essentially uses A.charpoly where A is the matrix which describes how to express $(x+y)x^iy^j$ in terms of the basis elements $x^iy^j$ for i,j less than the degree of their respective minimal polynomials. Unfortunately this seems quite a bit less efficient than https://en.wikipedia.org/wiki/Resultant which is the standard reference for the proof that algebraic numbers are closed under addition, and it's also not obviously related (one is the determinant of a $mn\times mn$ matrix, the other one is the determinant of a completely different $(m+n)\times(m+n)$ matrix)

I may have missed something, but why is it necessary to prove these facts constructively? Is it for users that want to use the tactic to produce constructive proofs?

[Mario Carneiro]

it's not, but having a constructive proof means that it is really easy to make some light modifications and recombine the lemmas to obtain a verified algorithm instead

[Mario Carneiro]

For example, if I wanted to write an algebraic numbers representation, I'd want to use minimal polynomials and add them using resultants, but then I'd get stuck either defining the operation or proving correctness without having a proof that the resultant has the required properties. If there was a constructive proof of closure of algebraic numbers which used resultants, the necessary lemma would be right there and I could reuse it

[Mario Carneiro]

the existing proof is "constructive enough", but the actual algorithm that it implicates is kind of bad