Tactic for determinants

Stream: mathlib4

[Jonathan Reich]

While working on Cartan matrix theorems for exceptional Lie algebras (#33115), we encountered a limitation with decide for computing determinants of larger matrices.
The situation:

• G₂.det = 1 (2×2): decide works 

• F₄.det = 1 (4×4): decide works

• E₆.det = 3 (6×6): decide fails (max recursion depth)

• E₇.det = 2 (7×7): decide fails

• E₈.det = 1 (8×8): decide fails

These are integer matrices with small entries (only 0, ±1, ±2, ±3), so the computation is straightforward by hand. We've added these as proof_wanted for now.
Is there interest in developing a more efficient determinant tactic? Possible approaches:

• Cofactor expansion with better evaluation strategy
• LU decomposition for integer matrices
• Exploiting sparsity (Cartan matrices are tridiagonal-ish)

The matrices are defined in Mathlib.Algebra.Lie.CartanMatrix if anyone wants to experiment. Happy to help test any proposed solutions.

[Ruben Van de Velde]

How far do simp and norm_num get? (Possibly specifying the "determinant is a sum" lemma)

[Jonathan Reich]

Thanks @Ruben Van de Velde. I tested a few approaches:
simp/norm_num: Haven't found the right lemma combination yet. Do you have a specific determinant expansion lemma in mind?
native_decide (discovered while testing):

• E₆.det = 3 works

• E₇.det = 2 works

• E₈.det = 1 fails (too many IR definitions)

The difference from what I can gather from decide is that native_decide compiles to native code, bypassing kernel recursion limits.

So we can fill in E₆ and E₇ with native_decide. E₈ still needs something smarter - maybe the determinant expansion lemmas you mentioned?

[Robin Carlier]

The native_decide tactic is not allowed in mathlib.

[Jonathan Reich]

Ah, good to know native_decide isn't allowed in mathlib. Thanks for the heads up.

So we're still looking for a proper solution. simp alone doesn't get far. Do you have a specific lemma in mind for the determinant expansion? Something like repeatedly applying Matrix.det_succ_column_zero for cofactor expansion?

For now the PR keeps these as proof_wanted. Happy to implement whatever approach is recommended.

[Li Xuanji]

That was a fun nerd snipe :) it occurred to me that you could write the L and U matrices for say E₈ and have Lean verify that L*U = E₈. The time complexity for a general matrix M would be O(n^3) (bottleneck is verifying that LU = M), which is the same as finding the decomposition, but the code to do so can live outside Lean.

So here's a proof that E₆.det = 3, E₇.det = 2.

Unfortunately L and U are necessarily matrices over ℚ, the proof for det_map_ℚ is pretty horrible, and there's something strange going on with prod_helper that screws up the proof for E₈.

[Yaël Dillies]

Li Xuanji said:

The time complexity for a general matrix M would be O(n^3)

Technically it's O(n^2.38) :slight_smile:

[Jonathan Reich]

Thanks @Li Xuanji - this is really elegant! The LU verification approach is exactly what's needed.

For the E₈ prod_helper issue - have you tried just exact Matrix.det_mul L₈ U₈? Curious if it's a timeout or something more subtle with the 8×8 case? Anyway, I'm sure you tried this.

Are you going to PR this directly to close out the proof_wanted items, or should we coordinate? Either way is good with me =)

[Li Xuanji]

Feel free to take my code and do whatever you want, including opening a PR! I have no plans to do so

[Eric Wieser]

Li Xuanji said:

the proof for det_map_ℚ is pretty horrible

This golfs to

lemma det_map_ℚ (n : ℕ) (M : Matrix (Fin n) (Fin n) ℤ) : (M.map (Int.castRingHom ℚ)).det = M.det :=
  Int.cast_det _ |>.symm

(docs#Int.cast_det)

[Oliver Nash]

It's very impressive, and great work by @Li Xuanji that we have proofs for the determinants. Thanks!

However I'm not absolutely sure we want these in Mathlib (I might feel different if / when we need to know these results for an application). What we really want is a tactic.

[Oliver Nash]

Incidentally if we were to go down the adhoc route of custom proofs for these matrices, it seems to me that E₈ could be the easiest: because its determinant is 1, the inverse will be a matrix of integers so if we supply this data in the proof and calculate that the product of E₈ with its inverse is 1 then using docs#Matrix.det_mul we immediately get that the determinant of E₈ divides 1 in the integers, so we have it up to a sign.

[Jonathan Reich]

det_lu.lean
I've spent a bit of time on this. 

The LU approach works almost perfectly. I have:

• L₈ and U₈ defined as the lower and upper triangular factors over ℚ

• E₈_lu : E₈Q = L₈ * U₈ — compiles fine

• L₈_det : L₈.det = 1 — proven via det_of_lowerTriangular + norm_num

• U₈_det : U₈.det = 1 — proven via det_of_upperTriangular + norm_num

• E₈_det_eq_one_or_neg_one — proven using the integer inverse (which also compiles)

The problem is the final connection step. Any attempt to use E₈_lu in a subsequent proof causes a stack overflow in the kernel. I tried:

• congrArg Matrix.det E₈_lu then rewriting

• Direct rw [E₈_lu]

• calc chains

• Various maxRecDepth and maxHeartbeats settings (up to 50k/2M)

All hit stack overflow, not just recursion depth limits.

I also tried the modular arithmetic route: compute det(E₈) in ZMod 3 (where 1 ≠ -1) using decide. This also fails with stack overflow / timeout as the 8! = 40,320 permutations in the Leibniz formula are too much even for a small finite field.

The proof is verified in pieces. The kernel just can't verify the equality E₈Q = L₈ * U₈ efficiently enough to use it in a subsequent theorem. This really does seem to need a dedicated determinant tactic that can exploit sparsity or use LU decomposition without requiring kernel verification of the full matrix product.

For now I've left it as proof_wanted with the LU infrastructure in place for when such a tactic exists. If others have more thoughts on how to make progress on this, please let me know.

[Li Xuanji]

Yeah, I can see why having random LU factorizations in Mathlib might not be desired.

Would Mathlib welcome a tactic to compute determinants of Matrix (Fin n) (Fin n) ℚ (0 < n) which internally computes the LU factorization? I can forsee that users not understanding why it doesn't work for finite fields, or ℚ[i], or ℚ[√2]... I'm not totally sure what the most general type over which this would work is, it seems to be some kind of field with computable field operations and equality.

Technically it's O(n^2.38) :slight_smile:

I meant that my proof that L*U = E₆ is O(n^3) if generalized to larger matrices. It would be fun to write a different O(n^2.38) proof and see if it has an actual speedup for concrete matrices :smile:

Curious if it's a timeout or something more subtle with the 8×8 case?

I get a (kernel) deep recursion detected even in term mode

theorem prod_helper : (L₈ * U₈).det = L₈.det * U₈.det := Matrix.det_mul L₈ U₈

I think the definition of det might be getting unfolded for concrete matrices.

[Oliver Nash]

The (kernel) deep recursion detected error is pretty bad. I've opened #33590 to address this.

It is likely that deeper issues exist (e.g., perhaps some day somebody working with docs#MultilinearMap.alternatization on a large concrete matrix might hit a similar issue) but making docs#Matrix.det into a def makes sense in its own right so I think we can probably merge this.

[Oliver Nash]

Li Xuanji said:

Would Mathlib welcome a tactic to compute determinants of Matrix (Fin n) (Fin n) ℚ (0 < n) which internally computes the LU factorization? I can forsee that users not understanding why it doesn't work for finite fields, or ℚ[i], or ℚ[√2]... I'm not totally sure what the most general type over which this would work is, it seems to be some kind of field with computable field operations and equality.

Yes we absolutely would love to have such a tactic. Regarding coefficients, I expect it should be similar to norm_num: general coefficients are fine as long as you don't expect it to do anything fancy for non-integral elements of your ring.

[Yaël Dillies]

Probably this shouldn't be a tactic but a simproc! See the blog

[Li Xuanji]

Regarding coefficients, I expect it should be similar to norm_num: general coefficients are fine as long as you don't expect it to do anything fancy for non-integral elements of your ring.

The issue I am thinking of is with these two matrices

def A : Matrix (Fin 2) (Fin 2) ℚ := 1
def A' : Matrix (Fin 2) (Fin 2) ℝ := 1

#eval A.det
#eval A'.det

I'm not really sure how to compute that A' is lower triangular since (for e.g.) #eval A' 0 0 fails.

In any case, I am very unlikely to have time to try to write a simproc, so anyone reading should feel free to give it a shot :slight_smile:

[Oliver Nash]

I'm not sure how strictly to interpret your use of the word "compute" but I'll just note in passing that it is trivial to "prove" that A' is lower triangular:

import Mathlib.Data.Real.Basic
import Mathlib.LinearAlgebra.Matrix.Block

open Matrix

def A' : Matrix (Fin 2) (Fin 2) ℝ := 1

-- `A'` is upper triangular
example : A'.BlockTriangular id :=
  fun i j hij ↦ by aesop (add simp [A', Matrix.one_apply])

-- `A'` is lower triangular (unsure of idiomatic spelling)
example : A'ᵀ.BlockTriangular id :=
  fun i j hij ↦ by aesop (add simp [A', Matrix.one_apply])

[Bhavik Mehta]

I have a proof of the determinant of a matrix corresponding to E8 in the sphere packing project, and a prototype implementation for much larger matrices elsewhere too. (In fact I also have a prototype for the same thing with real interval matrices, giving an interval output too)

[Bhavik Mehta]

Oliver Nash said:

The (kernel) deep recursion detected error is pretty bad. I've opened #33590 to address this.

It is likely that deeper issues exist (e.g., perhaps some day somebody working with docs#MultilinearMap.alternatization on a large concrete matrix might hit a similar issue) but making docs#Matrix.det into a def makes sense in its own right so I think we can probably merge this.

This exact issue was discussed on zulip a few months ago (and first noticed here https://leanprover.zulipchat.com/#narrow/channel/287929-mathlib4/topic/kernel.20deep.20recursion.20detected/near/453282247): it's specific to large concrete real matrices, and making it a def does indeed fix this, but when I tried to do so last time it had other unfortunate consequences. Maybe those are gone now though?

[Oliver Nash]

Thanks @Bhavik Mehta there seemed to be no negative consequences of making Matrix.det a def. It solves the issue here though not the root cause. I'll comment further in the other thread from ~18 months ago (active again less than three months ago) which I hadn't seen before.

[Johan Commelin]

I know this is not a tactic either, but I couldn't help myself:

variable {n : Type*} [Fintype n] [DecidableEq n] in
lemma E₈_det_aux {M : Matrix n n ℤ} {k : ℤ} (i j : n) (a b : ℤ)
    (h : (M.updateRow i (a • M i + b • M j)).det = a * k)
    (hij : i ≠ j := by decide) (ha : a ≠ 0 := by decide) :
    M.det = k := by
  simpa only [M.det_updateRow_add, M.det_updateRow_smul, M.det_updateRow_eq_zero hij.symm,
    M.updateRow_eq_self, add_zero, mul_zero, mul_eq_mul_left_iff, ha, or_false] using h

theorem E₈_det : E₈.det = 1 := by
  apply E₈_det_aux 2 0 (2:ℤ) (1:ℤ)
  apply E₈_det_aux 3 1 (2:ℤ) (1:ℤ)
  apply E₈_det_aux 3 2 (3:ℤ) (2:ℤ)
  apply E₈_det_aux 4 3 (5:ℤ) (1:ℤ)
  apply E₈_det_aux 5 4 (4:ℤ) (1:ℤ)
  apply E₈_det_aux 6 5 (3:ℤ) (1:ℤ)
  apply E₈_det_aux 7 6 (2:ℤ) (1:ℤ)
  rw [Matrix.det_of_upperTriangular]
  · decide +kernel
  · unfold BlockTriangular; decide +kernel

Note that the intermediate goals are completely unreadable.