computing submatrix for ring_theory.matrix

Stream: new members

[Tobias Grosser]

On another lean saturday afternoon, I continued looked forther into formalizing another basic lean maths. This time I want to compute sub-matrixes for ring_theory.matrix.

This is what I came up with:

def rsubmx {m n_left n_right: nat} :
Π (m n_left n_right),
matrix (fin m) (fin (n_right + n_left)) α → matrix (fin m) (fin (n_right)) α
| x y_l y_r A :=
  λ i j,
  let (offset : fin (y_r + y_l)) := (fin.add_nat j y_l) in
  A i offset

It compiles without flaws, but looks very messy. Two things I feel should still be changed are:
1) It would be nice if the function would not require m, n_right, n_left to be passed explicitly, but would instead automatically derive the relevant offsets.
2) I keep getting type errors due to a fin(x + y) being not equal to fin(y + x) when using it. It would be great if one could automatically convert between the two.

[Tobias Grosser]

Feedback very much appreciated.

[Kevin Buzzard]

I think that I would have been tempted to move away from just adding an offset and would instead consider any increasing maps fin m' -> fin m and fin n' -> fin n. Isn't that what people usually look at -- all the m' x n' minors or something?

[Tobias Grosser]

I am not sure I understand what interface you have in mind.

[Bryan Gin-ge Chen]

For example, suppose I have a 5x5 matrix and want to extract the 3x3 matrix consisting of rows 1,2,4 and columns 1,3,4.

[Tobias Grosser]

Right.

[Tobias Grosser]

I mean, it makes sense to have a more generic interface that computes a minor.

[Tobias Grosser]

Such a functionality could even be the basis for this function, I assume.

[Tobias Grosser]

Still, I feel that it would not completely replace this function as the interface for the special case I propose could be simpler, as it could just derive the sizes from the types.

[Tobias Grosser]

More from a practical perspective, I am not sure how a general minorification interface should look like. Would I ask for two more parameters fin m -> Prop and fin n -> Prop, which evaluate to true if the specific column should be included?

[Tobias Grosser]

Computing the necessary offsets seems a little involved. Not sure how this will effect performance.

[Bryan Gin-ge Chen]

I think the interface Kevin suggests is the following: to specify a general $m' \times n'$ submatrix of an $m \times n$ matrix, it is sufficient to provide two increasing maps, one of type fin m' → fin m and one of type fin n' → fin n. No need to compute any offsets, so far as I can tell; simply pull out the specified matrix entries using the two maps to grab the row and column indices.

[Kevin Buzzard]

And then the function you've already implemented could be a special case of this.

[Tobias Grosser]

Right.

[Tobias Grosser]

Here what I came up with:

import ring_theory.matrix

variables {α : Type}
variables {m n m' n': Type} [fintype m] [fintype n] [fintype m'] [fintype n']

def minormx : matrix m n α → (m' → m) → (n' → n) → matrix m' n' α :=
λ A trans_col trans_row i j, A (trans_col i) (trans_row j)

def rsubmx {m n_left n_right: nat} :
  matrix (fin m) (fin (n_right + n_left)) α → matrix (fin m) (fin (n_right)) α
| A := minormx A (λ i, i) (λ j, fin.add_nat j _)

[Tobias Grosser]

(deleted)

[Tobias Grosser]

(deleted)

[Tobias Grosser]

I created a bunch of tests:

def rsubmx_test_join_fail:
matrix (fin 5) (fin 5) α → matrix (fin 5) (fin 3) α :=
λ A, rsubmx A

def rsubmx_test_split_const_ok (x : nat) :
matrix (fin 5) (fin (3 + 2)) α → matrix (fin 5) (fin 3) α :=
λ A, rsubmx A

def rsubmx_test_split_const_fail (x : nat) :
matrix (fin 5) (fin (2 + 3)) α → matrix (fin 5) (fin 3) α :=
λ A, rsubmx A

def rsubmx_test_split_param_ok (x y : nat) :
matrix (fin 5) (fin (x + y)) α → matrix (fin 5) (fin (x)) α :=
λ A, rsubmx A

def rsubmx_test_split_param_fail (x y : nat) :
matrix (fin 5) (fin (y + x)) α → matrix (fin 5) (fin (x)) α :=
λ A, rsubmx A

Only if the fine type is split into a clear "+" with the correct order of operands in the plus (the LHS operand must appear in the result type), this compiles. Automatic splitting or exploiting commutativity does not work.

Here the associated error messages:

test.lean:16:5: error

type mismatch at application
  rsubmx A
term
  A
has type
  matrix (fin 5) (fin 5) α
but is expected to have type
  matrix (fin 5) (fin (3 + ?m_1)) α

type mismatch at application
  rsubmx A
term
  A
has type
  matrix (fin 5) (fin (2 + 3)) α
but is expected to have type
  matrix (fin 5) (fin (3 + ?m_1)) α
test.lean:24:5: error

type mismatch at application
  rsubmx A
term
  A
has type
  matrix (fin 5) (fin (y + x)) α
but is expected to have type
  matrix (fin 5) (fin (x + ?m_1)) α

[Tobias Grosser]

@Mario Carneiro , any idea how to make my tests independent of the order of fin? This seems to require lean foo I am not yet capable of.

[Mario Carneiro]

Use Kevin's suggestion

[Mario Carneiro]

You should avoid algebra inside types like the plague

[Tobias Grosser]

With

def minormx : matrix m n α → (m' → m) → (n' → n) → matrix m' n' α :=
λ A trans_col trans_row i j, A (trans_col i) (trans_row j)

I tried to implement Kevin's suggestion.

[Tobias Grosser]

I now try to implement ssreflects

   [l|r]submx A == the left/right submatrices of a row block matrix A.
                   Note that the type of A, 'M_(m, n1 + n2) indicates how A
                   should be decomposed.

on top of it.

[Tobias Grosser]

It seems a nice feature that they can split just based on the type sizes.

[Tobias Grosser]

Are you saying the interface they use in coq is not a good interface for lean?

[Mario Carneiro]

That will only work when n1 and n2 are explicitly given, and they also have to be numerals not variables for it to be useful

[Mario Carneiro]

Do you have any examples of use of this notation in Coq?

[Mario Carneiro]

Note that the type of A, 'M_(m, n1 + n2) indicates how A should be decomposed.

This sounds to me like lean's behavior is consistent with coq's: You have to give an input of a matrix matrix 5 (3 + 2) rather than matrix 5 5 to get the function to figure out n1 and n2

[Tobias Grosser]

Sure, my guaussien elimination code:

Fixpoint Gaussian_elimination {m n} : 'M_(m, n) → 'M_m × 'M_n × nat :=
  match m, n with
  | _.+1, _.+1 ⇒ fun A : 'M_(1 + _, 1 + _) ⇒
    if [pick ij | A ij.1 ij.2 != 0] is Some (i, j) then
      let a := A i j in let A1 := xrow i 0 (xcol j 0 A) in
      let u := ursubmx A1 in let v := a^-1 *: dlsubmx A1 in
      let: (L, U, r) := Gaussian_elimination (drsubmx A1 - v ×m u) in
      (xrow i 0 (block_mx 1 0 v L), xcol j 0 (block_mx a%:M u 0 U), r.+1)
    else (1%:M, 1%:M, 0%N)
  | _, _ ⇒ fun _ ⇒ (1%:M, 1%:M, 0%N)
  end.

it uses dlsubmx which again is defined as:

Definition ulsubmx := lsubmx (usubmx A).
Definition ursubmx := rsubmx (usubmx A).
Definition dlsubmx := lsubmx (dsubmx A).
Definition drsubmx := rsubmx (dsubmx A).

which again uses

Fact lsubmx_key : unit.
Definition lsubmx (A : 'M[R]_(m, n1 + n2)) :=
  \matrix[lsubmx_key]_(i, j) A i (lshift n2 j).

Fact rsubmx_key : unit.
Definition rsubmx (A : 'M[R]_(m, n1 + n2)) :=
  \matrix[rsubmx_key]_(i, j) A i (rshift n1 j).

Fact usubmx_key : unit.
Definition usubmx (A : 'M[R]_(m1 + m2, n)) :=
  \matrix[usubmx_key]_(i, j) A (lshift m2 i) j.

Fact dsubmx_key : unit.
Definition dsubmx (A : 'M[R]_(m1 + m2, n)) :=
  \matrix[dsubmx_key]_(i, j) A (rshift m1 i) j.

Lemma lshift_subproof m n (i : 'I_m) : i < m + n.

Lemma rshift_subproof m n (i : 'I_n) : m + i < m + n.

Definition lshift m n (i : 'I_m) := Ordinal (lshift_subproof n i).
Definition rshift m n (i : 'I_n) := Ordinal (rshift_subproof m i).

[Tobias Grosser]

Note that the type of A, 'M_(m, n1 + n2) indicates how A should be decomposed.

This sounds to me like lean's behavior is consistent with coq's: You have to give an input of a matrix matrix 5 (3 + 2) rather than matrix 5 5 to get the function to figure out n1 and n2

If that's as good as it gets that's fine. I just had hoped lean would be able to figure out that 5 can be split into 3 + 2. If this not the case, I need to do this manually. Is there a function that does this type splitting for me?

[Mario Carneiro]

what is the type of A1?

[Mario Carneiro]

or xrow and xcol

[Mario Carneiro]

Lean can figure out that 5 is 3 + 2, but it is ambiguous what to split it into

[Tobias Grosser]

I see.

[Tobias Grosser]

In my implementation in lean it is A1 : matrix (fin (x + 1)) (fin (y + 1)) α

[Tobias Grosser]

So what I have today might already be enough.

[Tobias Grosser]

I learn that generalizing this further is not a good idea.

[Tobias Grosser]

Thanks! I leave it at this generality.

[Tobias Grosser]

Now I am only left with the problem of casting. I need the matrix to be of type matrix (fin (1 + x)) (fin (1 + y)))

[Tobias Grosser]

I tried to type-cast using:

    let (A1': matrix (fin (1+x)) (fin (1+y)) α) := A1 in

but this just gives

type mismatch at application
  _let_match A1'
term
  A1'
has type
  matrix (fin (1 + x)) (fin (1 + y)) α
but is expected to have type
  matrix (fin (x + 1)) (fin (y + 1)) α

[Tobias Grosser]

You said lean should be able to understand this.

[Tobias Grosser]

This is no ambiguity

[Tobias Grosser]

I just need to show that these two are equivalent types.

[Mario Carneiro]

this is not ambiguous, but since x is not a numeral these aren't defeq

[Mario Carneiro]

How do the xrow i 0 (xcol j 0 A) functions work?

[Mario Carneiro]

are you building a new row on A or taking a minor?

[Tobias Grosser]

They just interchange rows. AFAIU they don;t thange the type.

[Mario Carneiro]

The fact that in coq 1+x = succ x is definitional makes this translation more complicated, unless you reverse the + arguments in the minors function, or work from the bottom rather than the top of the matrix

[Mario Carneiro]

In lean x+1 = succ x is definitional

[Tobias Grosser]

My lean implementation is:

def xrow [decidable_eq m] (row1: m) (row2: m) (A: matrix m n α) : matrix m n α :=
λ x y, if x = row1
         then
           A row2 y
         else
           if x = row2
             then
               A row1 y
             else
               A x y

[Mario Carneiro]

right, I see, that shouldn't be a problem

[Mario Carneiro]

you are still using the code you posted above for rsubmx?

[Mario Carneiro]

with n_right + n_left in the type I mean

[Tobias Grosser]

Yes.

[Tobias Grosser]

I could move the zero columns/rows to position x+1 / y +1, then the types would work out.

[Tobias Grosser]

But then the result matrices are layed out in a non-standard way.

[Mario Carneiro]

so since A has type (n+1) (n+1) you can naturally strip off the left/top end

[Tobias Grosser]

Right.

[Tobias Grosser]

That's easy. Just that most textbook implementations of gaussian elimination build upper triangular matrices.

[Tobias Grosser]

So the zeros should be at the lower left corner.

[Mario Carneiro]

I think you can still build upper triangular if you also transpose the rows and columns

[Tobias Grosser]

Sure, I could probably hack my way though this.

[Tobias Grosser]

Just don't feel this is nice.

[Mario Carneiro]

But I think that having reversed arguments to the minors function is also fine

[Tobias Grosser]

Why again does this work in coq?

[Mario Carneiro]

because they defined addition by recursion on the first argument rather than the second

[Mario Carneiro]

Personally I prefer lean's definition (i.e. x + succ y = succ (x + y))

[Tobias Grosser]

They do?

[Tobias Grosser]

I feel they also have this cast in their code.

[Tobias Grosser]

 | _.+1, _.+1 ⇒ fun A : 'M_(1 + _, 1 + _) ⇒

[Mario Carneiro]

That's why in the induction you notice the type is 1+_ instead of _ +1

[Mario Carneiro]

The .+ is suspicious

[Mario Carneiro]

do you know what it means?

[Tobias Grosser]

Notation "n .+1" := (succn n) (at level 2, left associativity,
  format "n .+1") : nat_scope.

[Mario Carneiro]

heh

[Tobias Grosser]

So they define notation to allow RHS addition of +1?

[Mario Carneiro]

so 1+n = n .+1

[Mario Carneiro]

looks like it

[Tobias Grosser]

This stuff is all crazy!

[Tobias Grosser]

Now based on the way nat is defined I cannot write my algorithm the way I want it.

[Tobias Grosser]

:(

[Mario Carneiro]

You can also write fin.add_nat with an appropriate type

[Mario Carneiro]

You can also cast to get A1' the way you were attempting, but this will cause more problems later on

[Tobias Grosser]

OK, I will try to write add_nat after breakfast.

There is already

def add_nat {n} (i : fin n) (k) : fin (n + k) :=
⟨i.1 + k, nat.add_lt_add_right i.2 _⟩

[Tobias Grosser]

Seems trivial to write (assuming one knows what to write)

[Tobias Grosser]

Now, I am not sure what I actually want here.

[Tobias Grosser]

In some way I just want to cast from n + 1 to 1 + n.

[Tobias Grosser]

Or pattern match based on 1 + n.

[Tobias Grosser]

Would add_nat help me with the pattern matching?

[Mario Carneiro]

just swap the arguments, like this

def nat_add {n} (k) (i : fin n) : fin (k + n) :=
⟨k + i.1, nat.add_lt_add_left i.2 _⟩

[Mario Carneiro]

then use it in your minors function

[Tobias Grosser]

Are you suggeting to change the signature of my functions.

def ursubmx {m_bottom m_top n_left n_right: nat} : Π (m_bottom m_top n_left n_right),
  matrix (fin (m_bottom + m_top)) (fin (n_left + n_right)) α → matrix (fin m_top) (fin (n_right)) α
:= sorry

def ursubmx2 {m_bottom m_top n_left n_right: nat} : Π (m_bottom m_top n_left n_right),
  matrix (fin (m_top + m_bottom)) (fin (n_right + n_left)) α → matrix (fin m_top) (fin (n_right)) α
:= sorry

[Mario Carneiro]

yes... although I think this will make a cast required

[Tobias Grosser]

I feel I have enough ideas to make things work. I can then later try to make the code nice.

[Tobias Grosser]

This just screems ugly all over the place.

[Tobias Grosser]

Need to get used to this more to see a beautiful solution.

[Tobias Grosser]

Need to have breakfast now. Thanks a lot for the helpful discussion. I learned a lot!

[Mario Carneiro]

Here's a trick: You can use Kevin's minors function as a cast

[Mario Carneiro]

the identity function fin (m+n) -> fin (n+m) is monotonic :)

[Kevin Buzzard]

Is there an argument for not even using fin n at all and having a general class indexed over a pair of finite types?

Here's a result I'd like to see in Lean one day: if f(X) is the char poly of the n x n matrix M then the coefficient of X^(n-i) in f(X) is the sum of the (appropriately signed?) determinants of the i x i principal minors.

[Mario Carneiro]

We already have that; Tobias is explicitly avoiding that because he wants to do induction on fin n

[Kevin Buzzard]

Oh -- I'm behind the times :-)

[Mario Carneiro]

If this is to be computationally efficient, it should probably use list representation though (a.k.a fast_matrix), and it might be possible to just ignore the dependent types and work with list (list A)

[Tobias Grosser]

I see.

[Tobias Grosser]

I would certainly like this to be computationally efficient, but will take this a step at a time. I feel I am close to sth that could work. I will finish this up and polish it.

[Tobias Grosser]

Then we have a first baseline implementation. I am then happy to take feedback on how to improve interface / performance / ...

[Mario Carneiro]

I think that's a good idea

[Tobias Grosser]

OK, I think I got a first implementation: https://gist.github.com/tobig/376e9a394c674474b8c1f6ecf9555478
It compiles and gives expected results. Still a couple of sorry, but seems to work.

[Tobias Grosser]

Need to clean up everything at some point.

[Tobias Grosser]

But no today any more.

[Tobias Grosser]

Also I learned it is indeed _very_ slow.

[Tobias Grosser]

I need to wait 10s of seconds for results on a 3x3 matrix.

[Tobias Grosser]

Oh dear, seems there is optimization potential.

[Kevin Buzzard]

I think for your test matrix I would have been tempted to send x and y to C[3*x+y] with C a vector of the nine entries (def C := [3,3,3,3,2,3,2,1,1])...hmm...one might then have to prove 3*x+y<=8...

[Andrew Ashworth]

hmm, are you aiming for speed? I believe vector has special support in the VM, so you might try defining matrix using that instead of an arbitrary fintype

[Andrew Ashworth]

iirc, and you would be more familiar than me on this topic, but doesn't SSReflect also define matrix as a vector of vectors?

[Tobias Grosser]

I am not aiming for speed atm. I just want to get things working / defined and improve my lean-knowledge.

[Tobias Grosser]

Eventually, I would like this to be fast and maybe even support exporting code to C++.

[Tobias Grosser]

However, this is a medium term only.

[Tobias Grosser]

My first pull request: https://github.com/leanprover/mathlib/pull/387

Still very small, but gets things rolling.

[Tobias Grosser]

I got an amazing review from @Simon Hudon and @Sean Leather. Thanks a lot!

[Simon Hudon]

The PR has my stamp of approval and, I think, that of @Sean Leather, now we just need @Mario Carneiro or @Johannes Hölzl to merge it.

[Sean Leather]

Agreed.

[Mario Carneiro]

yeah okay

[Simon Hudon]

This must be some sort of record :D

[Tobias Grosser]

Impressive!

[Tobias Grosser]

I pushed anther pull request for fin in (https://github.com/leanprover/mathlib/pull/388) as well as the submatrix definitions this thread was originally about in: https://github.com/leanprover/mathlib/pull/389. The last PR uses @Kevin Buzzard's idea of using a general function minormx to implement this functionality.

[Tobias Grosser]

Especially the last one is more a RFC, as I am not sure if limiting these functionality to (fin n) makes sense.

[Tobias Grosser]

My lean time is over for today. Will look at potential feedback later tonight. Thanks everyone again!