1

I want to calculate a product of huge (specific) matrices. From a point complexity of view, the product should be taken the form of an elementwise expression.

I tried to "vectorize" the matrices with mxvec / vec_mx and calculate the product via one dimensional streams. But indices access was blocked by the term of enum ('I_p * 'I_q).

I want to know a nth value of enum ('I_p * 'I_q) because I want to decscribe a multiplication of matrices in the form of a primitive expression in an underlying field.

How do I do this? In particular, how do I prove this statement?

From mathcomp Require Import all_ssreflect.

Lemma nth_enum_prod p q (a : 'I_q) :
  val a = index (ord0, a) (enum (prod_finType (ordinal_finType p.+1) (ordinal_finType q))).
John Kugelman
  • 349,597
  • 67
  • 533
  • 578

1 Answers1

0

I'm surprised you need to vectorize the matrices if your definition is point-wise, usually you should be able to define your result as \matrix_(i, j) op, for example the standard definition of matrix multiplication is:

\matrix_(i, k) \sum_j (A i j * B j k).

By the way, a quick "dirty" proof of your lemma is:

Lemma nth_enum_prod p q (a : 'I_q) : val a = index (@ord0 p, a) (enum predT).
Proof.
have /(_ _ 'I_q) pair_snd_inj: injective [eta pair ord0] by move => n T i j [].
have Hfst : (ord0, a) \in [seq (ord0, x2) | x2 <- enum 'I_q].
  by move=> n; rewrite mem_map /= ?mem_enum.
rewrite enumT !unlock /= /prod_enum enum_ordS /= index_cat {}Hfst.
by rewrite index_map /= ?index_enum_ord.
Qed.

but indeed if you find yourself using this it means you are into a different kind of problem. I just posted it as an illustration on how to manipulate this kind of expressions.

edit: based on your comment, a more principled way to manipulate the above is to define a lemma about index and products; I've left the full proof as an exercise, but the outline is:

Lemma index_allpairs (T U : eqType) (x : T) (y : U) r s :
  (* TODO: Some conditions are missing here *)
  index (x,y) [seq (x,y) | x <- r , y <- s] =
  size s * (index x r) + index y s.
Proof.
Admitted.

Lemma index_ord_allpairs p q (x : 'I_p) (y : 'I_q) :
  index (x,y) [seq (x,y) | x <- enum 'I_p , y <- enum 'I_q] = q * x + y.
Proof. by rewrite index_allpairs ?mem_enum ?size_enum_ord ?index_enum_ord. Qed.

Lemma nth_enum_prod p q (a : 'I_q) : val a = index (@ord0 p, a) (enum predT).
Proof. by rewrite enumT unlock index_ord_allpairs muln0. Qed.
ejgallego
  • 6,709
  • 1
  • 14
  • 29
  • The reason why I need this solution is that enormous zero will appear in the matrix expression. In addition, these matrix expressions will be extracted to optimized codes. Your proof have been very helpful for me! Thank you – Kazunari Tanaka Mar 08 '20 at 02:02
  • @KazunariTanaka still this solution looks a bit messy to me, I'd say that maybe sparsity should be handled at a lower-level once you have done some refinement of the matrix itself, _à la CoqEAL_ but indeed hard to say without discussing more. – ejgallego Mar 08 '20 at 05:14
  • End of fighting couple hours, I finally got it what you said. I'm going to try to read the paper of refinements... Thank you again! – Kazunari Tanaka Mar 08 '20 at 06:03