Over on twitter there was a blog post about whether you could use dependent types not just for tracking the sizes of matrices, but the units-of-measure in their entries. The poster seemed pessimistic, but I think it probably works out nicely as long as you could get the domain-specific equation-solving necessary for units-of-measure in the first place (ACU plus inverses) into some dependently typed language.

Then an n x m matrix M just needs to also carry a list of units-of-measure T of length n, and another such list U of length m, such that the unit-of-measure of M[i][j] is T[i]U[j]. Matrix multiplication goes like
Mat T U ⇒ Mat U-1 S ⇒ Mat T S
and inverses go
Mat T U ⇒ Mat U-1 T-1
(where --1 is pointwise inversion of every unit-of-measure in the list) and the determinant of a square matrix Mat T U has type (ΠT)(ΠU) --- that is, the product of all the units of measure from both lists.

I started writing out Gauss-Jordan in sml to see how you'd actually type algorithms, and I think it'd be pretty cute:

```(*

Vec [] = unit
Vec (u :: U) = u * Vec U

Bivec U1 U2 = Vec U1 * Vec U2

Bimat [] U1 U2 = unit
Bimat (t :: T) U1 U2 = Bivec (t . U1) (t . U2) * Bimat T U1 U2

t . [] = []
t . (u :: U) = (t . u) :: (t . U)

Mat [] U = unit
Mat (t :: T) U = Vec (t . U) * Mat T U

*)

(* aug_row : {n m : Nat} -> Vec T -> Vec 1s *)
fun aug_row n m [] = []
| aug_row n m (h::tl) = (if n = m then 1.0 else 0.0) :: aug_row n (m+1) tl;

(* aug : {n : Nat} {Z : Vec n *} -> Mat T U -> Bimat T U (Z @ (1 / T)) *)
fun aug' n [] = []
| aug' n (h::tl) = (h, aug_row n 0 h) :: aug' (n+1) tl;

(* aug : Mat T U -> Bimat T U (1 / T) *)
fun aug m = aug' 0 m

(* smul : t -> Vec T -> Vec (t * T)  *)
fun smul s [] = []
| smul (s:real) (h::tl) = s * h :: smul s tl

(* vplus : Vec T -> Vec T -> Vec T *)
fun vplus [] [] = []
| vplus ((h1:real)::tl1) (h2::tl2) = h1 + h2 :: (vplus tl1 tl2)
| vplus _ _ = raise Match

(* normalize : {n:Nat} -> Bivec T1 T2 -> Bivec (T1 / T1[n]) (T2 / T1[n]) *)
fun normalize n (v1, v2) =
let
val s = 1.0 / List.nth(v1, n)
in
(smul s v1, smul s v2)
end

(* tri_pass : {n:Nat} -> Bivec _1 _2 -> Bimat _3 _4 _5 *)
fun tri_pass _ _ [] = []
| tri_pass m (v1, v2) ((w1, w2)::tl) =
let
val s = ~(List.nth(w1, m))
in
(vplus w1 (smul s v1), vplus w2 (smul s v2)) :: tri_pass m (v1, v2) tl
end

fun tri m (h::tl) =
let
val nv = normalize m h
in
nv :: tri (m+1) (tri_pass m nv tl)
end
| tri m [] = []

val x = aug_row 0 0 [1.0, 2.0];
val x = aug [[1.0, 2.0, 3.0], [5.0, 6.0, 7.0], [8.0, 9.0, 5.0]];
val z = tri 0 x;
```
Tags:
• Post a new comment

#### Error

Anonymous comments are disabled in this journal

default userpic