Stream: math-comp users

Topic: matrix multiplication


view this post on Zulip Bas Spitters (Jan 04 2023 at 14:37):

I want to multiply two concrete matrices in math-comp, would I need coqeal for that?

view this post on Zulip Pierre Roux (Jan 04 2023 at 15:15):

I would say yes

view this post on Zulip Pierre Roux (Jan 04 2023 at 15:17):

(of course if it's like two 2x2 matrices once in a proof, you can probably do the proof by hand but if you want to perform actual computations, the refinement offered by CoqEAL seems a much better solution) Note that CoqEAL offers a refinement of MathComp matrices by lists of lists. Using arrays of arrays could probably be more efficient now that Coq offers primitive arrays.

view this post on Zulip Laurent Théry (Jan 04 2023 at 16:55):

Just a naive question, how big are the matrices you want to multiply?

view this post on Zulip Karl Palmskog (Jan 04 2023 at 20:47):

is there a refinement in CoqEAL from MathComp matrices to primitive arrays of arrays? Might be worth opening an issue otherwise (as memento for someone to implement it)

view this post on Zulip Pierre Roux (Jan 04 2023 at 21:19):

no (and yes, why not)

view this post on Zulip Karl Palmskog (Jan 04 2023 at 21:55):

https://github.com/coq-community/coqeal/issues/72

view this post on Zulip Bas Spitters (Jan 05 2023 at 07:31):

Laurent Théry said:

Just a naive question, how big are the matrices you want to multiply?

Just 4 by 4, so an ad hoc solution might also work.

view this post on Zulip Laurent Théry (Jan 05 2023 at 09:11):

something like that?

From mathcomp Require Import all_ssreflect all_algebra.

Section Matrix44.

Local Open Scope ring_scope.
Import GRing.Theory.

Variable R: ringType.

Definition seq2matrix m n (l: seq (seq R)) :=
  \matrix_(i<m,j<n) nth 1 (nth [::] l i) j.

Local Notation "''M{' l } " := (seq2matrix _ _ l).

Lemma matrix_example :
  'M{[:: [:: 1%:R;  2%:R;  3%:R;  4%:R];
         [:: 5%:R;  6%:R;  7%:R;  8%:R];
         [:: 9%:R; 10%:R; 11%:R; 12%:R];
         [::13%:R; 14%:R; 15%:R; 16%:R]]}
*
  'M{[:: [::17%:R; 18%:R; 19%:R; 20%:R];
         [::21%:R; 22%:R; 23%:R; 24%:R];
         [::25%:R; 26%:R; 27%:R; 28%:R];
         [::29%:R; 30%:R; 31%:R; 32%:R]]}
=
  'M{[:: [:: 250%:R;  260%:R;  270%:R;  280%:R];
         [:: 618%:R;  644%:R;  670%:R;  696%:R];
         [:: 986%:R; 1028%:R; 1070%:R; 1112%:R];
         [::1354%:R; 1412%:R; 1470%:R; 1528%:R]]}
          :> 'M[R]_(4,4).
Proof.
apply/matrixP=> i j.
rewrite !mxE !big_ord_recl big_ord0 !mxE.
by case: i=> [[|[|[|[|]]]]] //= _;
  case: j=> [[|[|[|[|]]]]] //= _;
  rewrite -!natrM -(natrD _ _ 0) -!natrD.
Qed.

view this post on Zulip Bas Spitters (Jan 05 2023 at 11:04):

Yes. Thanks!

view this post on Zulip Maxime Dénès (Jan 05 2023 at 13:04):

Karl Palmskog said:

is there a refinement in CoqEAL from MathComp matrices to primitive arrays of arrays? Might be worth opening an issue otherwise (as memento for someone to implement it)

I had done it in my thesis. But it probably required the fold operator on primitive arrays which didn't pass the review phase when integrating primitive arrays.

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:07):

@Maxime Dénès do you have some code we could link to in the issue? Sounds like this might be a nice student course/internship project to pick up the old code and get it into CoqEAL properly.

view this post on Zulip Maxime Dénès (Jan 05 2023 at 13:07):

It's in Coq that you'd need to get it.

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:09):

you mean, it's in some old Coq commit?

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:10):

or you mean the whole refinement can't be done without the folds?

view this post on Zulip Maxime Dénès (Jan 05 2023 at 13:12):

it can probably be done through unary integers and casting to primitive integers, but that wouldn't make much sense IMHO

view this post on Zulip Maxime Dénès (Jan 05 2023 at 13:13):

not efficient, and you'd need to redo it once the proper iterators land

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:14):

is there some issue in Coq open about iterators for primitive arrays? Should we create one? At least it would be nice to document these things.

view this post on Zulip Pierre-Yves Strub (Jan 05 2023 at 13:23):

yes, having a fold operator over primitives array would be great. I ended up with that in some dev:

Fixpoint iofold_ {T} (f : int -> T -> T * bool) (i M : int) (x : T) (n : nat) :=
  if n is n.+1 then
    if M <=? i then (i, (x, false)) else
      let maybe b i x :=
        if b then iofold_ f i M x n else (i, (x, b)) in

      let: (i, (x, b)) := ((i + 1)%uint63, f i x) in
      let: (i, (x, b)) := maybe b i x in
      let: (i, (x, b)) := maybe b i x in
      (i, (x, b))
  else (i, (x, true)).


Lemma fuel : {x : nat | x = Uint63.size}.
Proof. by exists Uint63.size. Qed.

Definition iofold {T} (f : int -> T -> T * bool) (M : int) (x : T) :=
  (iofold_ f 0 M x (tag fuel)).2.1.

view this post on Zulip Pierre-Yves Strub (Jan 05 2023 at 13:24):

At least, you only need 63 sticks, not 2^63 :)

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:30):

https://github.com/coq/coq/issues/17067

view this post on Zulip Pierre-Yves Strub (Jan 05 2023 at 13:32):

I can make this code public if needed.

view this post on Zulip Pierre-Yves Strub (Jan 05 2023 at 13:32):

(We have also fold operators for arrays + related lemmas)

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:36):

I don't personally use matrices in Coq projects, but perhaps @Bas Spitters would be interested? The code might go into CoqEAL?

view this post on Zulip Pierre-Yves Strub (Jan 05 2023 at 13:37):

The Coq will be released under a meta approved repo on github under a free license (MIT I think). Then, you'll be free to copy it whereever you want.

view this post on Zulip Karl Palmskog (Jan 05 2023 at 13:40):

if someone is going to rely on it going forward, we'd need to get it into some Coq Platform project. Fine by me to copy from some other repo to say, CoqEAL. Also, CoqEAL is MIT so no license fiddling required.

view this post on Zulip Bas Spitters (Jan 05 2023 at 14:23):

@Laurent Théry One downside of your lemma is that one still needs to supply the answer manually.
Ideally, Coq would compute the answer itself. (Your answer suffices for what we had in mind though).

view this post on Zulip Karl Palmskog (Jan 05 2023 at 14:41):

there was a recent paper on matrices in Coq: https://link.springer.com/chapter/10.1007/978-3-031-21213-0_11 ("Integration of Multiple Formal Matrix Models in Coq")

but, from what I can tell, no code is provided anywhere. :crying_cat:

view this post on Zulip Bas Spitters (Jan 05 2023 at 15:13):

I don't need their code now, but if you feel like it, it wouldn't be too hard to send them an email...

view this post on Zulip Alexander Gryzlov (Jan 05 2023 at 15:22):

Isn't this it? https://github.com/zhengpushi/coq-matrix

view this post on Zulip Pierre Roux (Jan 05 2023 at 15:26):

All comments seem to be in Chinese unfortunately :-(

view this post on Zulip Karl Palmskog (Jan 05 2023 at 15:38):

ah, that seems indeed to be the code, but the next dilemma: GPL-3.0 makes it incompatible with MathComp

view this post on Zulip Laurent Théry (Jan 06 2023 at 10:21):

Yes and no. You can always trick Coq with evars. I am not an expert with manipulating evar but here is an old style solution.

Lemma compute16 l1 l2 A v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16 :
let M : 'M[R]_(4,4) :=   'M{[:: [:: v1%:R;  v2%:R; v3%:R; v4%:R];
         [:: v5%:R;  v6%:R;  v7%:R;  v8%:R];
         [:: v9%:R; v10%:R; v11%:R; v12%:R];
         [:: v13%:R; v14%:R; v15%:R; v16%:R]]} in
  'M{l1} * 'M{l2} = M ->
  ('M{l1} * 'M{l2} = M -> A) -> A.
Proof. by move=> M H1 H2; apply: H2. Qed.

Ltac compute16Tac l1 l2  :=
apply: (compute16 l1 l2);
first by
let i := fresh "i" in
let j := fresh "j" in
let u := fresh "u" in
(apply/matrixP=> i j;
rewrite !mxE !big_ord_recl big_ord0 !mxE;
case: i => [[|[|[|[|]]]]] //= _;
  case: j=> [[|[|[|[|]]]]] //= _;
  rewrite -!natrM -(natrD _ _ 0) -!natrD;
  set u := (_ + _)%N; compute in u; rewrite /u; reflexivity).

Goal True.
compute16Tac
([:: [:: 1%:R;  2%:R;  3%:R;  4%:R];
         [:: 5%:R;  6%:R;  7%:R;  8%:R];
         [:: 9%:R; 10%:R; 11%:R; 12%:R];
         [::13%:R; 14%:R; 15%:R; 16%:R]] : seq (seq R))
([:: [::17%:R; 18%:R; 19%:R; 20%:R];
         [::21%:R; 22%:R; 23%:R; 24%:R];
         [::25%:R; 26%:R; 27%:R; 28%:R];
         [::29%:R; 30%:R; 31%:R; 32%:R]] : seq (seq R)) => H.
Qed.

view this post on Zulip Bas Spitters (Jan 07 2023 at 17:10):

@Laurent Théry Thanks. That's an interesting trick.

view this post on Zulip Bas Spitters (Jan 07 2023 at 17:15):

I've invited them to join:
https://github.com/zhengpushi/coq-matrix/issues/1

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 01:53):

Karl Palmskog 发言道

there was a recent paper on matrices in Coq: https://link.springer.com/chapter/10.1007/978-3-031-21213-0_11 ("Integration of Multiple Formal Matrix Models in Coq")

but, from what I can tell, no code is provided anywhere. :crying_cat:

I am the author. Thank you for your quote. Here is the code : https://github.com/zhengpushi/coq-matrix
This formal matrix project is still alive, I will upgrade a newer version later. (1. setoid equality insted of eq, 2. use Typeclasses to simplify the code, 3. Hierarchy of the matrix element. )

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 01:54):

Pierre Roux 发言道

All comments seem to be in Chinese unfortunately :-(

I'm so sorry. I will fix these poor comment gradually.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 02:03):

Bas Spitters 发言道

Laurent Théry Thanks. That's an interesting trick.

Thank for your invitation.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 04:02):

formal matrix model, different implementations. Part 1.

(**
   purpose : Quick view of the different formal matrix models
   author  : ZhengPu Shi (zhengpushi@nuaa.edu.cn)
   remark  :
   1. there are five formal matrix models I known now.
   (1) DR: dependent record
   (2) DL: dependent list
   (3) DP: dependent pair
   (4) NF: function with natural indexing
   (5) FF: function with fin indexing
 *)

Require List.
Require ZArith.

Declare Scope mat_scope.
Delimit Scope mat_scope with mat.
Open Scope mat_scope.


(** * model 1 : DR : Depencent Record *)
(** Main idea: use (list (list A)) to store the matrix data *)
Module DR.

  Import List ListNotations.

  (** matrix models *)
  Section matrix.
    Context {A : Type}.

    (** All list in a dlist have same length. Note: dlist means list (list A)  *)
    Definition width (dl : list (list A)) (c : nat) : Prop :=
      Forall (fun l => length l = c) dl.

    (** Definition of matrix type *)
    Record mat (r c : nat) := {
        mat_data : list (list A);
        mat_length : length mat_data = r;
        mat_width : width mat_data c
      }.

    Arguments Build_mat {r c}.
    Arguments mat_data {r c}.
    Arguments mat_length {r c}.
    Arguments mat_width {r c}.

    (** Construct a dlist by column with a list and a dlist *)
    Fixpoint consc (l : list A) (dl : list (list A)) : list (list A) :=
      match l, dl with
      | xl :: tl, xdl :: tdl => (xl :: xdl) :: (consc tl tdl)
      | _, _ => []
      end.

    (** a list list that every list is nil, named as dnil *)
    Fixpoint dnil (n : nat) : list (list A) :=
      match n with
      | O => nil
      | S n' => nil :: (dnil n')
      end.

    (** Transposition of the content of a matrix *)
    Fixpoint dltrans (dl : list (list A)) (c : nat) : list (list A) :=
      match dl with
      | [] => dnil c
      | l :: tl => consc l (dltrans tl c)
      end.

    (** dltrans length law *)
    Lemma dltrans_length : forall dl c, width dl c -> length (dltrans dl c) = c.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    (** dltrans width law *)
    Lemma dltrans_width : forall dl r c,
        length dl = r -> width dl c -> width (dltrans dl c) r.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    (** Transposition of a matrix *)
    Definition mtrans {r c} (m : mat r c) : mat c r.
      refine (Build_mat (dltrans (mat_data m) c) _ _).
      apply dltrans_length. apply mat_width.
      apply dltrans_width. apply mat_length. apply mat_width.
    Defined.

    (** map operation to two list *)
    Fixpoint map2 (f : A -> A -> A) (l1 l2 : list A) : list A :=
      match l1, l2 with
      | x1 :: t1, x2 :: t2 => (f x1 x2) :: (map2 f t1 t2)
      | _, _ => []
      end.

    (** Assume that we have a zero number *)
    Context (A0 : A).

    (** list to fixed-length list (vector) *)
    Fixpoint l2v_aux (l : list A) (n : nat) : list A :=
      match n with
      | 0 => []
      | S n' => (hd A0 l) :: (l2v_aux (tl l) n')
      end.

    (** list list to fixed-shape list list (matrix) *)
    Fixpoint l2m_aux (dl : list (list A)) (r c : nat) : list (list A) :=
      match r with
      | 0 => []
      | S n' => (l2v_aux (hd nil dl) c) :: (l2m_aux (tl dl) n' c)
      end.

    Lemma l2m_aux_length : forall (dl : list (list A)) (r c : nat),
        length (l2m_aux dl r c) = r.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    Lemma l2m_aux_width : forall (dl : list (list A)) (r c : nat),
        width (l2m_aux dl r c) c.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    (** dlist to matrix *)
    Definition l2m (dl : list (list A)) (r c : nat) : mat r c.
      refine (Build_mat (l2m_aux dl r c) _ _).
      apply l2m_aux_length.
      apply l2m_aux_width.
    Defined.

    (** matrix to dlist *)
    Definition m2l {r c} (m : mat r c) := mat_data m.

    (** Assume that we have two operations of the element: Addition, Multiplication *)
    Context (Aadd Amul : A -> A -> A).

    (** dot product, marked as l1 . l2 *)
    Definition ldot (l1 l2 : list A) : A :=
      fold_right Aadd A0 (map2 Amul l1 l2).

    (** list dot product to dlist *)
    Fixpoint ldotdl (l : list A) (dl : list (list A)) : list A :=
      match dl with
      | h :: t => (ldot l h) :: (ldotdl l t)
      | [] => []
      end.

    (** dlist dot product *)
    Fixpoint dldotdl (dl1 dl2 : list (list A)) : list (list A) :=
      match dl1 with
      | h1 :: t1 => (ldotdl h1 dl2) :: (dldotdl t1 dl2)
      | [] => []
      end.

    (** dldotdl length law *)
    Lemma dldotdl_length : forall dl1 dl2 r1,
        length dl1 = r1 -> length (dldotdl dl1 dl2) = r1.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    (** dldotdl width law *)
    Lemma dldotdl_width : forall dl1 dl2 r2 c,
        length dl2 = r2 -> width dl1 c -> width dl2 c ->  width (dldotdl dl1 dl2) r2.
    Proof. Admitted. (* this lemma is proved. details see the project *)

    (** Multiplication of two matrices *)
    Definition mmul {r c t : nat} (m1 : mat r c) (m2 : mat c t) : mat r t.
      refine (Build_mat (dldotdl (mat_data m1) (mat_data (mtrans m2))) _ _).
      apply dldotdl_length. apply mat_length.
      apply dldotdl_width with (c:=c). apply mat_length. apply mat_width.
      simpl. apply dltrans_width. apply mat_length. apply mat_width.
    Defined.

  End matrix.

  Section matrix_Z.
    Import ZArith.
    Open Scope Z_scope.
    Open Scope mat_scope.

    Notation l2m := (l2m 0).
    Notation mmul := (mmul 0 Zplus Zmult).

    Let m1 : mat 4 4 :=
          l2m [[1;2;3;4];
               [5;6;7;8];
               [9;10;11;12];
               [13;14;15;16]] 4 4.
    Let m2 : mat 4 4 :=
          l2m [[17;18;19;20];
               [21;22;23;24];
               [25;26;27;28];
               [29;30;31;32]] 4 4.

    Compute m2l (mmul m1 m2).
    (* = [[250; 260; 270; 280]; *)
    (*    [618; 644; 670; 696]; *)
    (*    [986; 1028; 1070; 1112]; *)
    (*    [1354; 1412; 1470; 1528]] *)
    (*  : list (list Z) *)
  End matrix_Z.

End DR.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 05:15):

Formal matrix model, different implementations. Part 2.

(**
   purpose : Quick view of the different formal matrix models
   author  : ZhengPu Shi (zhengpushi@nuaa.edu.cn)
   remark  :
   1. there are five formal matrix models I known now.
   (1) DR: dependent record
   (2) DL: dependent list
   (3) DP: dependent pair
   (4) NF: function with natural indexing
   (5) FF: function with fin indexing
 *)

Require List ZArith.
Require Vectors.Fin Vectors.Vector.

Declare Scope mat_scope.
Delimit Scope mat_scope with mat.
Open Scope mat_scope.


(** * model 2 : DL : Dependent List *)
(** Main idea: use Coq.Vectors.Vector *)
Module DL.
  Import Fin Vector.

  Notation fin := Fin.t.
  Notation vec := Vector.t.
  Arguments vec {A}.

  (** vector model *)
  Section vector.
    Import VectorNotations.
    Context {A : Type} (A0 : A) (Aadd Amul : A -> A -> A).

    (** Build vector with a function [gen: fin n -> A] *)
    Fixpoint vmake {n} : (fin n -> A) -> vec n :=
      match n with
      | O => fun _ => []
      | S n' => fun (gen : fin (S n') -> A) =>
                 (gen F1) :: (vmake (fun (fn':fin n') => gen (FS fn')))
      end.
  End vector.

  (** matrix model *)
  Section matrix.
    Import List ListNotations Vector VectorNotations.
    Context {A : Type} (A0 : A) (Aadd Amul : A -> A -> A).

    (** Definition of matrix type *)
    Definition mat r c := @vec (@vec A c) r.

    (** Convert between vector and list *)
    Fixpoint v2l {n} (v : @vec A n) : list A :=
      match v with
      | []%vector => []%list
      | (x :: v')%vector => (x :: (v2l v'))%list
      end.
    Fixpoint l2v (l : list A) (n : nat) : vec n :=
      match n with
      | 0 => []%vector
      | S n' => (List.hd A0 l) :: (l2v (List.tl l) n')
      end.

    (** Convert between matrix and list list *)
    Fixpoint m2l {r c} (m : mat r c) : list (list A) :=
      match m with
      | []%vector => []%list
      | (x :: v')%vector => (v2l x) :: (m2l v')
      end.
    Fixpoint l2m (dl : list (list A)) (r c : nat) : mat r c :=
      match r with
      | 0 => []
      | S n' => (l2v (List.hd List.nil dl) c) :: (l2m (List.tl dl) n' c)
      end.

    (** dot product of two vectors *)
    Definition vdot {n} (v1 v2 : vec n) := fold_left Aadd A0 (map2 Amul v1 v2).

    (** Get a column *)
    Definition mcol {r c} (m : mat r c) :=
      fun fc : fin c => vmake (fun fr : fin r => nth (nth m fr) fc).

    (** Multiplication of two matrices *)
    Definition mmul {r s c : nat} (m1 : mat r s) (m2 : mat s c) : mat r c :=
      vmake (fun fr : fin r =>
               vmake (fun fc : fin c =>
                        vdot (nth m1 fr) (mcol m2 fc)
        )).
  End matrix.

  Section matrix_Z.
    Import List ListNotations ZArith.
    Open Scope Z_scope.

    Notation l2m := (l2m 0).
    Notation mmul := (mmul 0 Zplus Zmult).

    Let m1 : mat 4 4 :=
          l2m [[1;2;3;4];
               [5;6;7;8];
               [9;10;11;12];
               [13;14;15;16]] 4 4.
    Let m2 : mat 4 4 :=
          l2m [[17;18;19;20];
               [21;22;23;24];
               [25;26;27;28];
               [29;30;31;32]] 4 4.

    Compute m2l (mmul m1 m2).
    (* = [[250; 260; 270; 280]; *)
    (*    [618; 644; 670; 696]; *)
    (*    [986; 1028; 1070; 1112]; *)
    (*    [1354; 1412; 1470; 1528]] *)
    (*  : list (list Z) *)
  End matrix_Z.

End DL.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 05:46):

Formal matrix model, different implementations. Part 3.

(**
   purpose : Quick view of the different formal matrix models
   author  : ZhengPu Shi (zhengpushi@nuaa.edu.cn)
   remark  :
   1. there are five formal matrix models I known now.
   (1) DR: dependent record
   (2) DL: dependent list
   (3) DP: dependent pair
   (4) NF: function with natural indexing
   (5) FF: function with fin indexing
 *)

Require List ZArith.
Require Vectors.Fin Vectors.Vector.

Declare Scope vec_scope.
Delimit Scope vec_scope with vec.
Open Scope vec_scope.

Declare Scope mat_scope.
Delimit Scope mat_scope with mat.
Open Scope mat_scope.


(** * model 3 : DP : Dependent Pair *)
(** Main idea: matrix is a vector of vectors, which with a fixpoint function *)
Module DP.
  Import List ListNotations.

  (** vector model *)
  Section vector.
    Context {A : Type} (A0 : A).

    (** Definition of a vector *)
    Fixpoint vec (n : nat) : Type :=
      match n with
      | O => unit
      | S n => prod A (vec n)
      end.

    (** Convert between vector and list *)
    Fixpoint v2l {n} : vec n -> list A :=
      match n with
      | 0 => fun (v : vec 0) => []
      | S n' => fun (v : vec (S n')) => (fst v) :: (v2l (snd v))
      end.
    Fixpoint l2v (l : list A) (n : nat) {struct n} : vec n :=
      match n as n0 return (vec n0) with
      | 0 => tt
      | S n' => (hd A0 l, l2v (tl l) n')
      end.

    Notation "[ ]" := (tt : vec 0) : vec_scope.
    Notation "[ x ]" := (pair x (vec 0)) : vec_scope.
    Notation "[ x1 ; .. ; xn ]" := (pair x1 .. (pair xn tt) .. ) : vec_scope.

  End vector.

  (** matrix model *)
  Section matrix.
    Context {A : Type} (A0 : A) (Aadd Amul : A -> A -> A).

    (** Definition of matrix type *)
    Definition mat (r c : nat) : Type := @vec (@vec A c) r.

    (** Convert between matrix and list list *)
    Fixpoint m2l {r c} : (mat r c) -> list (list A) :=
      match r with
      | O => fun _ => nil
      | S r' => fun (m : @vec (@vec A c) (S r')) =>
                 cons (v2l (fst m)) (m2l (snd m))
      end.
    Fixpoint l2m (dl : list (list A)) (r c : nat) : mat r c :=
      match r with
      | 0 => tt
      | S n' => (l2v A0 (hd nil dl) c, l2m (tl dl) n' c)
      end.

    (** Fold a vector to an element from left to right *)
    Fixpoint vfoldl {n} (Aadd : A -> A -> A) (A0 : A) : (vec n) -> A :=
      match n with
      | O => fun (v : vec 0) => A0
      | S n' => fun (v : vec (S n')) => Aadd (fst v) (vfoldl Aadd A0 (snd v))
      end.

    (** Maping of two vectors *)
    Fixpoint vmap2 {n} (f : A -> A -> A) : vec n -> vec n -> vec n :=
      match n with
      | O => fun (v1 : vec 0) (v2 : vec 0) => tt
      | S n' => fun (v1 : vec (S n')) (v2 : vec (S n')) =>
                 (f (fst v1) (fst v2), vmap2 f (snd v1) (snd v2))
      end.

    (** dot product of two vectors *)
    Definition vdot {n} (v1 v2 : vec n) := vfoldl Aadd A0 (vmap2 Amul v1 v2).

    (** Inner product a vector and a matrix. v(c) *' m(r,c) = v(r) *)
    Fixpoint vdotm {r c} (v : @vec A c) : mat r c -> @vec A r :=
      match r with
      | O => fun _ => tt
      | S r' => fun (m : mat (S r') c) => (vdot v (fst m), vdotm v (snd m))
      end.

    (** Inner product of two matrices. m(r1,c) *' (m(r2,c) = m(r1,r2) *)
    Fixpoint mdot {r c t} : (mat r c) -> (mat t c) -> (mat r t) :=
      match r with
      | O => fun _ _ => tt
      | S r' => fun (m1 : mat (S r') _) m2 => (vdotm (fst m1) m2, mdot (snd m1) m2)
      end.

    (** Get head column of a matrix *)
    Fixpoint mhdc {r c} : mat r (S c) -> vec r :=
      match r with
      | O => fun _ => tt
      | S r' => fun (m : mat (S r') (S c)) => (fst (fst m), mhdc (snd m))
      end.

    (** Get tail columns of a matrix *)
    Fixpoint mtlc {r c} : mat r (S c) -> mat r c :=
      match r with
      | O => fun _ => tt
      | S r' => fun (m : mat (S r') (S c)) => (snd (fst m), mtlc (snd m))
      end.

    (** Transpose a matrix *)
    Fixpoint mtrans {r c} : mat r c -> mat c r :=
      match c with
      | O => fun (m : mat r 0) => tt
      | S c' => fun (m : mat r (S c')) => (mhdc m, mtrans (mtlc m))
      end.

    (** Multiplication of two matrices *)
    Definition mmul {r c s} (m1 : mat r c) (m2 : mat c s) : mat r s :=
      mdot m1 (mtrans m2).

  End matrix.

  Section matrix_Z.
    Import List ListNotations ZArith.
    Open Scope Z_scope.

    Notation l2m := (l2m 0).
    Notation mmul := (mmul 0 Zplus Zmult).

    Let m1 : mat 4 4 :=
          l2m [[1;2;3;4];
               [5;6;7;8];
               [9;10;11;12];
               [13;14;15;16]] 4 4.
    Let m2 : mat 4 4 :=
          l2m [[17;18;19;20];
               [21;22;23;24];
               [25;26;27;28];
               [29;30;31;32]] 4 4.

    Compute m2l (mmul m1 m2).
    (* = [[250; 260; 270; 280]; *)
    (*    [618; 644; 670; 696]; *)
    (*    [986; 1028; 1070; 1112]; *)
    (*    [1354; 1412; 1470; 1528]] *)
    (*  : list (list Z) *)
  End matrix_Z.

End DP.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 06:44):

Formal matrix model, different implementations. Part 4.

(**
   purpose : Quick view of the different formal matrix models
   author  : ZhengPu Shi (zhengpushi@nuaa.edu.cn)
   remark  :
   1. there are five formal matrix models I known now.
   (1) DR: dependent record
   (2) DL: dependent list
   (3) DP: dependent pair
   (4) NF: function with natural indexing
   (5) FF: function with fin indexing
 *)

Require List ZArith.
Require Vectors.Fin Vectors.Vector.

Declare Scope vec_scope.
Delimit Scope vec_scope with vec.
Open Scope vec_scope.

Declare Scope mat_scope.
Delimit Scope mat_scope with mat.
Open Scope mat_scope.


(** * model 4 : NF : Function with Natural number indexing *)
(** Main idea: matrix is a function from two index to a value *)
Module NF.
  Import List ListNotations.

  (** matrix model *)
  Section matrix.
    Context {A : Type} (A0 : A) (Aadd Amul : A -> A -> A).

    (** Definition of matrix type *)
    Definition mat (r c : nat) : Type := nat -> nat -> A.

    (** Get i-th row of a matrix *)
    Fixpoint mrow_aux {r c : nat} (ri : nat) (cnt : nat) (m : mat r c) : list A :=
      match c with
      | O => nil
      | S c' => m ri cnt :: (@mrow_aux r c' ri (S cnt) m)
      end.
    Definition mrow {r c : nat} (ri : nat) (m : mat r c) := @mrow_aux r c ri 0 m.

    (** Convert between matrix and list list *)
    Fixpoint m2l_aux {r c : nat} (cnt : nat) (m : mat r c) : list (list A) :=
      match r with
      | O => nil
      | S r' => mrow cnt m :: (@m2l_aux r' c (S cnt) m)
      end.
    Definition m2l {r c} (m : mat r c) := m2l_aux 0 m.
    Definition l2m (r c : nat) (dl : list (list A)) : mat r c :=
      (fun x y => nth y (nth x dl []) A0).

    (** Sum of a sequence *)
    Fixpoint seqsum (f : nat -> A) (n : nat) : A :=
      match n with
      | O => A0
      | S n' => Aadd (seqsum f n') (f n')
      end.

    (** Multiplication of two matrices *)
    Definition mmul {r c t : nat} (m1 : mat r c) (m2 : mat c t) : mat r t :=
        (fun x z => seqsum (fun y => Amul (m1 x y) (m2 y z)) c).

  End matrix.

  Section matrix_Z.
    Import List ListNotations ZArith.
    Open Scope Z_scope.

    Notation l2m := (l2m 0).
    Notation mmul := (mmul 0 Zplus Zmult).

    Let m1 : mat 4 4 :=
          l2m 4 4
            [[1;2;3;4];
             [5;6;7;8];
             [9;10;11;12];
             [13;14;15;16]].
    Let m2 : mat 4 4 :=
          l2m 4 4
            [[17;18;19;20];
             [21;22;23;24];
             [25;26;27;28];
             [29;30;31;32]].

    Compute m2l (mmul m1 m2).
    (* = [[250; 260; 270; 280]; *)
    (*    [618; 644; 670; 696]; *)
    (*    [986; 1028; 1070; 1112]; *)
    (*    [1354; 1412; 1470; 1528]] *)
    (*  : list (list Z) *)
  End matrix_Z.

End NF.

view this post on Zulip Julien Puydt (Jan 09 2023 at 07:22):

@ZhengPu Shi You mention you want to use type classes in a later version -- I haven't looked at your code, but MC is moving to hierarchy builder for its structures, so perhaps you should move to that too? (else you risk needing another transition down the line)

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 07:38):

Julien Puydt 发言道

ZhengPu Shi You mention you want to use type classes in a later version -- I haven't looked at your code, but MC is moving to hierarchy builder for its structures, so perhaps you should move to that too? (else you risk needing another transition down the line)

Yes, the code above havn't use Typeclasses. These codes are demonstrate for matrix multiplication computation directly.
The MC project have good hierarchy and more complicated. I havn't compeletely understand the matrix design of it. But the typeclasses technologies and module system have attracted/puzzled me for a long time.

In the project (https://github.com/zhengpushi/coq-matrix) , I just made a simple integration of several different matrix implementations.
The main purpose is to hide the differences of matrix operations, to help the final proof engineers who need a matrix library.
Although the project is available, but the organization of the contents is poor.
Later (maybe less than one month), a new version code will be pulled to github. (I am still working)

view this post on Zulip Karl Palmskog (Jan 09 2023 at 08:20):

@ZhengPu Shi as you may know, the Mathematical Component library is licensed under the CECILL-B license. However, CECILL-B is incompatible with GPL-3.0 according to the authors of GPL. This means that we can't legally combine your work with MathComp and distribute the results (Coq .vo files). Is there any chance you would consider releasing the code under a different license (such as LGPL) to facilitate easier combined development?

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 08:59):

Karl Palmskog 发言道

ZhengPu Shi as you may know, the Mathematical Component library is licensed under the CECILL-B license. However, CECILL-B is incompatible with GPL-3.0 according to the authors of GPL. This means that we can't legally combine your work with MathComp and distribute the results (Coq .vo files). Is there any chance you would consider releasing the code under a different license (such as LGPL) to facilitate easier combined development?

OK. I will change the license later. This project is belong to me, and I don't really care about the license.

view this post on Zulip Karl Palmskog (Jan 09 2023 at 09:01):

if the choice of license is not a concern, then we recommend using the MIT license. This would also facilitate integrating code into CoqEAL, a quite widely used library on refinement of Coq data.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 09:08):

Karl Palmskog 发言道

if the choice of license is not a concern, then we recommend using the MIT license. This would also facilitate integrating code into CoqEAL, a quite widely used library on refinement of Coq data.

Thank you for your advise. I could re-create a project later to use MIT license.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 09:11):

Karl Palmskog 发言道

if the choice of license is not a concern, then we recommend using the MIT license. This would also facilitate integrating code into CoqEAL, a quite widely used library on refinement of Coq data.

In fact, I have no experience to contribute to Coq project. I expect to learn from the Coq community, and if I can contribute any help I will feel very happy.

view this post on Zulip Karl Palmskog (Jan 09 2023 at 10:19):

ZhengPu Shi said:

Thank you for your advise. I could re-create a project later to use MIT license.

There is no need to recreate a project just to change its license, you just need to make a commit to change the LICENSE file in your GitHub repository. For example, you can see here how a maintainer changed his license from GPL to the MIT license. Be sure not to leave any references to GPL inside your .v source files, though.

view this post on Zulip Pierre Roux (Jan 09 2023 at 12:33):

ZhengPu Shi said:

Pierre Roux 发言道

All comments seem to be in Chinese unfortunately :-(

I'm so sorry. I will fix these poor comment gradually.

Don't be ashamed, I think we still have a few comments written in French in Coq :wink:

About your paper, you write that MathComp matrices are not well suited for computations. I fully agree with that but we could have expected a citation of CoqEAL which provides, among others, a refinement of matrix by lists of lists: seqmx.v (CoqEAL is based on typeclasses). That's what we used for instance in our validsdp tactic to perform proofs with mathcomp data-structures while still enjoying effective computations.

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 14:18):

Karl Palmskog 发言道

ZhengPu Shi said:

Thank you for your advise. I could re-create a project later to use MIT license.

There is no need to recreate a project just to change its license, you just need to make a commit to change the LICENSE file in your GitHub repository. For example, you can see here how a maintainer changed his license from GPL to the MIT license. Be sure not to leave any references to GPL inside your .v source files, though.

Thank you for your help. I have switched the project with MIT license. coq-matrix

view this post on Zulip ZhengPu Shi (Jan 09 2023 at 14:24):

Pierre Roux 发言道

ZhengPu Shi said:

Pierre Roux 发言道

All comments seem to be in Chinese unfortunately :-(

I'm so sorry. I will fix these poor comment gradually.

Don't be ashamed, I think we still have a few comments written in French in Coq :wink:

About your paper, you write that MathComp matrices are not well suited for computations. I fully agree with that but we could have expected a citation of CoqEAL which provides, among others, a refinement of matrix by lists of lists: seqmx.v (CoqEAL is based on typeclasses). That's what we used for instance in our validsdp tactic to perform proofs with mathcomp data-structures while still enjoying effective computations.

Thank you for your comment and advise. I also expect the MathComp could be work with computation. I will test the seqmx.v later.

view this post on Zulip Karl Palmskog (Jan 09 2023 at 14:26):

I added some tracking of the general matrix computation issue: https://github.com/coq-community/coqeal/issues/72 https://github.com/coq-community/manifesto/issues/143

This is something we'd like to provide to Coq users for sure in the Coq Platform (translate seamlessly between matrix representations in different libraries, compute quickly). I don't personally use matrices in my Coq projects, but I will try help out from the ecosystem side with packaging and so on.

view this post on Zulip Alexander Gryzlov (Jan 10 2023 at 10:55):

I imagine ultimately you'd want some form of BLAS/LAPACK intergration and some form of mutable state reasoning on top of that for your (quasi)CAS

view this post on Zulip Bas Spitters (Jan 10 2023 at 13:05):

@Vadim Zaliva 's Helix does such numerical computations. I imagine there are matrices in there somewhere.
https://github.com/vzaliva/helix

view this post on Zulip ZhengPu Shi (Jan 11 2023 at 01:57):

I have quickly checked the project and the author's Phd thesis. The Helix is a code generation technology from Coq to C, with some immediate representations.
The matrix they used https://github.com/vzaliva/helix/blob/master/coq/Util/Matrix.v , which is a vector of vectors (come from Coq.Vectors.Vector).
It looked liked that they havn't developed matrix theory too much (only one file, 26 lines). Maybe they mainly focused on Vector first.
This model is also integrated by us. see my demo code.

view this post on Zulip Julin Shaji (Oct 16 2023 at 17:58):

For matrix multiplication with mulmx, should the matrix elements necessarily be of a semiring type?

 HB.about mulmx.
(*
HB: mulmx is canonically equipped with structures:
    - mathcomp.algebra.ssralg.GRing.Nmodule_isSemiRing.axioms_
      (from "./algebra/matrix.v", line 2725)
    - GRing.Linear
      (from "./algebra/matrix.v", line 3447)
    - GRing.Additive
      (from "./algebra/matrix.v", line 3242)
*)

Last updated: Jul 23 2024 at 21:01 UTC