Matrices in stdlib -- how to represent them?

[LionOfJewdah]

I have a personal project where I have a construction of matrices over (commutative) (semi-) rings. I'd like to contribute some form of it to the standard library, but as I discussed with @JacquesCarette in this discussion in agda-categories, I have to do some design work and seek input from stdlib maintainers and Agda veterans. (I am a neophyte -- I hadn't heard of Agda until January. Love it.)

In my repo, I use something like this, Fin n and Fin m indexed function-based matrices.

module Matrix {ℓ ℓ≈} (R : CommutativeRing ℓ ℓ≈) where
open import Data.Nat
open import Algebra.Construct.Pointwise
open import Algebra.Module.Bundles

-- the bi(semi-)module or even Hadamard (C)(S)Ring of `Mat n m` is just Pointwise.ring (etc) 
--`(Fin n → Fin m) R`. Simplified :
HadamardMatRing (n m : ℕ) : Ring _ _
HadamardMatRing n m = Pointwise.ring (Fin n → Fin m) (CommutativeRing.ring R) -- simplified

Mat : (n m : ℕ) : Module R _ _
-- define s *ₗ A and A *ᵣ s via the Hadamard lift `const const s`; s : R; coincidence is trivial

For matrix multiplication define e.g.

δ : ∀ {n} → Fin n → Fin n → Carrier
δ zero    zero    = 1#
δ zero    (suc _) = 0#
δ (suc _) zero    = 0#
δ (suc i) (suc j) = δ i j

Iₘ : ∀ {n} → Mat n n
Iₘ = mkMat δ

mulSummand : ∀ {n k m : ℕ} → Mat n k → Mat k m → Fin n → Fin m → Fin k → Carrier
mulSummand A B i j = (entry A i -⟨ _*_ ⟩- flip (entry B) j) ˢ id

innerProd : ∀ {n k m : ℕ} → Mat n k → Mat k m → Fin n → Fin m → Carrier
innerProd A B i j = finSum (mulSummand A B i j)

_·ₘ_ : ∀ {n k m : ℕ} → Mat n k → Mat k m → Mat n m
A ·ₘ B = mkMat (innerProd A B)

and you get the (semi)ring structure for SqMat n = Mat n n almost for free. You can of course define trace, determinant, inverse, and lots of other constructs of linear algebra.

That said, there are definite pros and cons of representing matrices as (records wrapping) functions Fin n → Fin m → R.carrier.

As some pros: Finite ordinals are well-ordered, decidable, natural (sometimes literally defined as the natural numbers), and Data.Fin has lots of nice combinators like splitAt, remQuot, join, punchOut that make linear algebra convenient. Not to mention Data.Fin.{Subset,Properties,Permutation} etc which are also a goldmine. There is also only one "5" or any other n, we have UIP, and you can easily define submatrices, horizontal and vertical concatenation, and hence the Kronecker tensor product and direct sum.

I also have an implementation of matrices as Vec n (Vec m) R.carrier (with round trips), which may be better for algorithms like Gaussian elimination and RREF construction.

That said, a lot of the things we want matrices to be don't require Fin n / Fin m, and also I don't know where on the spectrum of generality would be most useful. On one end, you could lift any commutative (additive) semigroup S to matrices of n x m elements that you can add and take trace of and nothing else (including scale). On the other end, you could define them over fields (although stdlib only has RealHeytingField, not a general Field type). I'm inclined to think the "sweet spot" would be over commutative rings, which give the determinant, subtraction, and lots of nice properties, perhaps generalized to commutative semirings which lose the determinant but retain most everything else. However, I think matrices over general (semi)rings are worth considering as long as we're okay losing some nice properties in the most general case. (In my repo, I actually do have matrices over not-necessarily commutative semirings, and commutative semirings, and rings, and commutative rings, but that is hard to follow architecturally.)

95%+ of the time, even in my math-heavy repo, I'm using matrices over commutative rings. You get trace cyclicity, the general linear group, Cramer's rule with a proof that det A is invertible, and lots more nice things.

What are some thoughts on how to best design this?

[MatthewDaggitt]

I did put together a small PR many years ago using the vector of vectors appraoch and received some pushback on my design: #1525

Looking at it again now, whatever the representation of indices we choose I think we should probably try and go for tensors first, and have matrices as a special case of 2D tensor?

[jamesmckinna]

Patrik Jansson (Chalmers) has done a lot of work in Agda on tensor algebra for (teaching) physics, including relating abstract representations and concrete / coordinatised forms . I'll try and find a good link: eg https://research.chalmers.se/en/publication/546169

As for fields, much as for matrices, various attempts have been made, but in a constructive setting, and without necessarily being able to decide the underlying equality, the natural notion is better that of local ring. #2219

[gallais]

Instead of having matrices, I am tempted to suggest defining a universe of finite types

open import Data.Nat.Base using (ℕ)

data Shape : Set where
  `Fin : ℕ → Shape
  _`⊎_ : Shape → Shape → Shape
  _`×_ : Shape → Shape → Shape

and then have a generic notion of Tensor indexed by such a Shape

open import Data.Fin.Base using (Fin)
open import Data.Sum.Base using (_⊎_)
open import Data.Product.Base using (_×_)

Index : Shape → Set
Index (`Fin n) = Fin n
Index (i₁ `⊎ i₂) = Index i₁ ⊎ Index i₂
Index (i₁ `× i₂) = Index i₁ × Index i₂

record Tensor {a} (sh : Shape) (A : Set a) : Set a where
  constructor mkTensor
  field runTensor : Index sh → A

There is of course an isomorphic flat representation as a big vector
which is worth spelling out too.

Matrix multiplication ought to be an instance of a more general operation
for pairs of tensors of the form

mult : Tensor (s `× t) A → Tensor (t `× u) A → Tensor (s `× u) A

[Taneb]

There's another perspective on matrices, that they're (representations of) homomorphisms between certain modules. This provides a reason for the algebraic operations on them.

[jwaldmann]

We have an application where we need to extract code that actually multiplies matrices (the transition relation of weighted automata)

• we cannot use a functional representation of vectors/arrays because it would re-compute everything in a multiplication chain, and this blows up exponentially https://git.imn.htwk-leipzig.de/waldmann/cetera/-/work_items/23
• vector indexed by Fin _ (as in stdlib) is inefficient because it is still a cons-list, and has linear cost
• we got best results with list-of-list and mtimes a b = map (\ row -> map (dot row) (transpose b)) a (avoiding the indexing, but incurring the cost of transpose)
• in the end, we may need sparse matrices (using Data.AVL perhaps)

This is quite orthogonal to how you want to prove things.

[gallais]

we cannot use a functional representation of vectors/arrays because it would re-compute everything in a multiplication chain

A naïve such representation sure, but you could have a tabulate : Tensor s -> Tensor s
function that is semantically the identity but whose behaviour is to memoise the
computation.

we got best results with list-of-list

There should not be a runtime difference between list-of-list and vect-of-vect as they
are compiled to the same type.

[jwaldmann]

tabulate function that is semantically the identity but whose behaviour is to memoise the
computation.

but when I use that memoized function, I pay the cost of indexing? so that must be O(1)

not be a runtime difference between list-of-list and vect-of-vect

yes, if using map (\ row -> map (dot row) (transpose b)) a and there's no indexing, just traversal (via map/fold)