Teach Lean that (Matrix) types are equivalent?

Stream: new members

[Jake Jongejans]

I am working with a lot of matrices, and one problem I am continuously running into is that I need to reindex the matrices I have to give them the correct types. This was fine with a single reindex that I defined, but now, I need to keep track of 3 types of reindexes, and it is growing difficult to manage.

I was wondering if it is possible to teach Lean that these types are identical so reindex is not necessary anymore, or if there is a better way to apply reindex, so I only have one single custom reindex command. Here's some code:

import Mathlib.Data.Matrix.Basic
import Mathlib.Data.Matrix.Notation
import Mathlib.Data.Matrix.Kronecker
import Mathlib.Data.Complex.Basic

open Matrix
open Kronecker

variable { m n : ℕ }

-- Num of qubits
abbrev QCount (n : ℕ) := Fin (2 ^ n)
abbrev nMatrix (n : ℕ) := Matrix (QCount n) (QCount n) ℂ
abbrev nRMatrix (n : ℕ) := Matrix (QCount n) (QCount n) (Finset ℂ)
abbrev mnMatrix (m n : ℕ) := Matrix (QCount m × QCount n) (QCount m × QCount n) ℂ

def QCount_mul_QCount { m n : ℕ } : (QCount m × QCount n) ≃ QCount (m + n) := by
  simp [QCount]
  rw [Nat.pow_add]
  exact @finProdFinEquiv (2 ^ m) (2 ^ n)

def QCount_mul_one { m : ℕ } : (QCount (m * 1)) ≃ QCount m := by
  simp [QCount]
  rfl

/-- Reindex a circuit matrix to Fin 2 ^ n × Fin 2 ^ n dimensions -/
def reindex (A : mnMatrix m n) : nMatrix (m + n) :=
  Matrix.reindex QCount_mul_QCount QCount_mul_QCount A

/-- Deal with * 1 in the matrix dimensions -/
def reindex' (A : nMatrix (m * 1)) : nMatrix m :=
  Matrix.reindex QCount_mul_one QCount_mul_one A


def Fin_mul_Fin_one { m : ℕ } : Fin 1 × Fin m ≃ Fin m := by
  exact Equiv.uniqueProd (Fin m) (Fin 1)

/-- This one is from trying to prove that (1 : Nat) ⊗ₖ M = M... -/
def reindex'' (M : Matrix (Fin 1 × Fin m) (Fin 1 × Fin n) ℂ) : Matrix (Fin m) (Fin n) ℂ :=
  Matrix.reindex Fin_mul_Fin_one Fin_mul_Fin_one M

As can be seen, I already have three different reindex definitions that I need to keep track of, and proving they work together with one another becomes a terrible mess as it grows exponentially.

[Kevin Buzzard]

Can you make Lean's more flexible index system work for you rather than just following the books who will insist on matrices being labelled by {1,...,n}?

[Jake Jongejans]

I came to my abbreviations (nMatrix, Count, etc.) based on different things I tried making. Working with the kronecker product gives back different types of matrices based on its input, and the current way I am handling it seemed the easiest to me.
Kronecker product example: given A B : Matrix (Fin n) (Fin n) \C, kroneckerMap A B has the type Matrix ((Fin n) x (Fin n)) ((Fin n) x (Fin n)) \C) instead of having (Fin (n * n)) or (Fin (n ^ 2)). When I am working with a tensor power, this became impossible to manage, because for each n, the type was different. Hence the reindex, but that only handles a single equivalence in Fin types.
I am not sure what I could do different to utilise the index system from Lean in a better way. Do I need to change out my tensor power to not be a recursive kroneckerMap? Am I missing something else?

[Yaël Dillies]

Jake Jongejans said:

given A B : Matrix (Fin n) (Fin n) ℂ kroneckerMap A B has the type Matrix ((Fin n) x (Fin n)) ((Fin n) x (Fin n)) ℂ instead of having (Fin (n * n)) or (Fin (n ^ 2)).

Can you simply not try to index by Fin (n * n)?

[Yaël Dillies]

Note that docs#Matrix.kronecker takes in more general things than A B : Matrix (Fin n) (Fin n) ℂ

[Jake Jongejans]

I could try, but I still think the kronecker product will mess with the dimensions of the matrix type, because it always handles it like so: given (A : Matrix (Fin a1) (Fin a2) \C) (B : Matrix (Fin b1) (Fin b2) \C)we get kroneckerMap A B as type Matrix ((Fin a1) x (Fin b1)) ((Fin a2) x (Fin b2)) \C, and I'd still need to use a reindex to get it to Matrix (Fin (a1 * b1)) (Fin (a2 * b2)) \C.

One place I was struggling is here:

import Mathlib.Data.Matrix.Basic
import Mathlib.Data.Matrix.Notation
import Mathlib.Data.Matrix.Kronecker
import Mathlib.Data.Complex.Basic

open Matrix
open Kronecker
open Complex

theorem kronecker_natOne {A B : Type} [DecidableEq A] [DecidableEq B] (M : Matrix A B ℂ) : (1 : Matrix (Fin 1) (Fin 1) ℂ) ⊗ₖ M = M := by
-- type mismatch
--    M
-- has type
--    Matrix A B ℂ : Type
-- but is expected to have type
--    Matrix (1 × A) (1 × B) ℂ : Type

I am trying to prove that kronecker with a 1 that is just a natural number (i.e. Matrix (Fin 1) (Fin 1)) with any matrix is equal to that matrix. but I cannot because of the type change.

[Yaël Dillies]

Jake Jongejans said:

it always handles it like so: given (A : Matrix (Fin a1) (Fin a2) ℂ (B : Matrix (Fin b1) (Fin b2) ℂ get kroneckerMap A B as type Matrix ((Fin a1) x (Fin b1)) ((Fin a2) x (Fin b2)) ℂ

No it doesn't? Can you give a concrete example where this happens?

[Jake Jongejans]

The last code I sent shows the problem I am having—I updated it to be ever more clear. I'm getting Matrix (1 x A) (1 x B) \C instead of Matrix A B \C.

[Eric Wieser]

(A misreading of this thread inspired me to make #12703)

[Eric Wieser]

Here's the working version of your statement:

import Mathlib.Data.Matrix.Basic
import Mathlib.Data.Matrix.Notation
import Mathlib.Data.Matrix.Kronecker
import Mathlib.Data.Complex.Basic

open Matrix
open Kronecker
open Complex

theorem kronecker_natOne {A B : Type} [DecidableEq A] [DecidableEq B] (M : Matrix A B ℂ) :
    (1 : Matrix Unit Unit ℂ) ⊗ₖ M = reindex (Equiv.punitProd _).symm (Equiv.punitProd _).symm M := by
  rw [one_kronecker]
  ext ⟨⟨⟩, i⟩ ⟨⟨⟩, j⟩
  simp [blockDiagonal]

[Kyle Miller]

Or leaving it with Fin 1:

theorem kronecker_natOne {A B : Type} [DecidableEq A] [DecidableEq B] (M : Matrix A B ℂ) :
    (1 : Matrix (Fin 1) (Fin 1) ℂ) ⊗ₖ M =
      Matrix.reindex (Equiv.uniqueProd A (Fin 1)).symm (Equiv.uniqueProd B (Fin 1)).symm M := by
  sorry

[Kyle Miller]

Or generalized to both:

theorem kronecker_natOne {U A B : Type} [Unique U] [DecidableEq A] [DecidableEq B] (M : Matrix A B ℂ) :
    (1 : Matrix U U ℂ) ⊗ₖ M =
      Matrix.reindex (Equiv.uniqueProd _ _).symm (Equiv.uniqueProd _ _).symm M := by
  sorry

[Kevin Buzzard]

I think that even if there is a better indexing set for these matrices than (bad-looking choices such as) Fin (a * b) (and there may or may not be), you are probably still going to run into trouble identifying (X x Y) x Z with X x (Y x Z) and Fin 1 x X with X (the "monoid laws" for the monoidal category of finite types)

[Jake Jongejans]

I agree, Kevin. I am trying to figure out how to avoid these troubles. Any suggestions?

[Eric Wieser]

I think the actual missing lemma here is about blockDiagonal:

theorem blockDiagonal_unique {α U A B : Type} [Unique U] [DecidableEq A] [DecidableEq B]
  (d : U → Matrix A B α) [Zero α]:
  blockDiagonal d =
    Matrix.reindex (Equiv.prodUnique _ _).symm (Equiv.prodUnique _ _).symm (d default) := by
  ext ⟨a, ua⟩ ⟨b, ub⟩
  obtain rfl := Subsingleton.elim ua default
  obtain rfl := Subsingleton.elim ub default
  rfl

theorem kronecker_natOne {U A B : Type} [Unique U] [DecidableEq A] [DecidableEq B] (M : Matrix A B ℂ) :
    (1 : Matrix U U ℂ) ⊗ₖ M =
      Matrix.reindex (Equiv.uniqueProd _ _).symm (Equiv.uniqueProd _ _).symm M := by
  rw [one_kronecker, blockDiagonal_unique]
  rfl

[Eric Wieser]

I think one answer here might be a prod_simp% helper that automatically glues uniqueProd, prodComm, prodAssoc etc together into the reindex args, but this is still painful when working with rewrites

[Kyle Miller]

@Jake Jongejans This is one of those fundamental issues with dependent types, and unfortunately we don't have any way to make association come for free.

However, we do have some automation here and there to help out. For example, there's prod_assoc%, which automatically constructs equivalences for any reassociation:

example (A B C D : Type) : (A × B) × C × D ≃ A × (B × C) × D := prod_assoc%

[Kyle Miller]

(@Adam Topaz Should we make prod_assoc% be aware of Unit? Or [Unique _] in general? Or do you think that would be out of scope?)

[Jake Jongejans]

I see, so there's no way of solving this outright. I will have to look into these helpers, see if they can make my life less of a pain. Any place where I can find out more on prod_simp% and prod_assoc%? I guess I can throw those terms in Loogle.

[Kevin Buzzard]

Jake Jongejans said:

I agree, Kevin. I am trying to figure out how to avoid these troubles. Any suggestions?

These troubles really are there, and to make proofs rigorous you might have to deal with these issues. Maybe you need explicit reindex maps corresponding to one_mul, mul_one, mul_comm and mul_assoc, but these bijections will be there already.

[Kyle Miller]

prod_simp% doesn't exist, Eric was just saying "wouldn't it be nice if it did"

[Kyle Miller]

For prod_assoc%, it's probably better looking at Mathlib/Tactic/ProdAssoc.lean and reading the documentation, and looking at test/ProdAssoc.lean to see examples

[Kyle Miller]

And welcome to the world of non-strict symmetric monoidal categories! It's known that there's a coherence theorem here, which is essentially that there's just one obvious map when reassociating and canceling "monoidal units" (types with one element) when dealing with cartesian products.

[Kyle Miller]

There are some techniques for dealing with associativity, but I don't know how well they'd work in practice for matrices. For example, rather than types, you work with functors. Then you can get defeq associativity:

def F := Type → Type

def F.lift (m : Type) : F := fun r => m × r

def F.prod (m n : F) : F := m ∘ n

def F.toType (m : F) : Type := m Unit

instance : Coe Type F := ⟨F.lift⟩

infixr:35 " ×f " => F.prod

example : ((Nat ×f Nat) ×f Nat).toType = (Nat ×f (Nat ×f Nat)).toType := rfl

Maybe this would be a workable definition of a matrix? I don't know:

def Matrix' (m n : F) (α : Type*) := ∀ a b, m a → n b → α

In any case, this would be going in a completely new direction from what's already there in the library.

[Jake Jongejans]

Ah, non-strict symmetric monoidal categories. I will do some reading on that, and see how I can perhaps implement it and utilise it. Stepping away from mathlib will be a big loss, though. Thank you for the suggested code!

[Eric Wieser]

Could you share what your original goal is here, and perhaps an equation from it that brought you to this pain?

[Jake Jongejans]

Sure, though it's quite some code in separate files. Would you like me to condense it and remove unnecessary code, or would it be better to set my repository to public and share the link, and then pinpoint where I am struggling?

[Eric Wieser]

I mean maybe just share a line or two of LaTeX (either in $$ $$s or in a screenshot) that captures the sort of equation where you run into trouble

[Jake Jongejans]

Ah, currently I am trying to create a set of theorems that help with formalising quantum algorithms (basic ones for now, like Bernstein-Vazirani), so it's less of a single goal and more of a family of theorems and definitions I am trying to create, so I am not quite capable of giving a single "goal".
I am incorporating definitions such as pow_kronecker (a tensor power), and qapply that allows for applying a gate (or circuit, at some point) to certain qubits. I haven't got the logic of the last part quite right, yet, though.
The idea is to build up a library of gates and circuits that can be applied, to define the working of algorithms, and to formally prove them through Lean. In the grand scheme of things, Quantum Hoare Logic would be implemented, but I have insufficient time to implement that, in my thesis.

Anyhow, one issue right now is, where the kronecker_natOne does not recognise the 1 I get from pow_kronecker' 0 M. But in general, the reindexing and different forms matrices have been difficult to work with.

def pow_kronecker' { m : ℕ } (n : ℕ) (M : nMatrix m) : nMatrix (m * n) :=
  match n with
    | 0 => 1
    | (n + 1) => reindex (pow_kronecker' n M ⊗ₖ M)

theorem pow_kronecker_zero' { M : nMatrix m } : pow_kronecker' 0 M = 1 := by
  simp [pow_kronecker']

theorem pow_kronecker_one' { i : nMatrix m } : Circuits.reindex' (pow_kronecker' 1 i) = i := by
  rw [pow_kronecker', pow_kronecker_zero']
  -- This next rewrite fails, because the type of my 1 and its 1 is different.
  rw [kronecker_natOne i, one_kronecker_one, smul_reindex, reindex_one]
  norm_num

Does this answer your request? I am happy to offer more information.

[Eric Wieser]

For that particular example it would be tempting to use

variable {m n : Type*} {R : Type*} [Semiring R]

def powKronecker (k : ℕ) (A : Matrix m n R) :
    Matrix (Fin k → m) (Fin k → n) R :=
  fun fm fn => ((List.finRange k).map fun i => A (fm i) (fn i)).prod

[Eric Wieser]

But I think ultimately you hit the same issue

[Jake Jongejans]

At some point, I also think I will come across a wall again. The mention of the monoidal categories did make me wonder whether my approach to the problem is wrong, but I couldn't say right now.

[Jake Jongejans]

Hey, it's me again. I hope reviving this thread is fine.

I am trying to understand a bit how to utilise these functors properly, using Kyle's code, and I am trying to implement that into basically a shadow version of the Mathlib Matrix Basic file. Because the matrix now isn't just a function anymore with two variables and one output, I am quite lost as I don't know how to retrieve any value from this new matrix shape anymore.
I thought I could do, given M : Matrix' m n a, M Unit Unit to get the matrix at its current size, then M Unit Unit i j for the i and j index of the matrix, but it keeps complaining about types. Maybe a small code example (note I renamed to Matrix without '):

def map (M : Matrix m n α) (f : α → β) : Matrix m n β :=
  of fun i j => f (M Unit Unit i j)

This gives me the following error:

application type mismatch
  M Unit Unit I
argument
  I
has type
  Type : Type 1
but is expected to have type
  m Unit : Type

Maybe I am missing something trivial but I have tried quite a few things already. Perhaps there is another part of Mathlib that utilises functors more? That would be great to study. Any help is greatly appreciated!

[Kyle Miller]

An #mwe would be helpful.

Here's a guess:

def map (M : Matrix m n α) (f : α → β) : Matrix m n β :=
  of fun x y i j => f (M x y i j)

[Kyle Miller]

Or in other words,

def map (M : Matrix m n α) (f : α → β) : Matrix m n β :=
  of (f ∘ M)

[Kyle Miller]

I hope this functor idea doesn't lead you too far down a dead end!