Expressive Linear Algebra in Haskell Henning Thielemann 2019-08-21 - PowerPoint PPT Presentation

Expressive Linear Algebra in Haskell Expressive Linear Algebra in Haskell Henning Thielemann 2019-08-21

Expressive Linear Algebra in Haskell Motivation 1 Motivation 2 Solution 3 More realistic problem 4 More features 5 Closing

Expressive Linear Algebra in Haskell Motivation Resistor cube R X , 1 , 1 R Z , 0 , 1 R Z , 1 , 1 R Y , 0 , 1 R Y , 1 , 1 R X , 1 , 0 R X , 0 , 1 R Y , 0 , 0 R Y , 1 , 0 R Z , 0 , 0 R Z , 1 , 0 Wanted: R X , 0 , 0 Total resistance

Expressive Linear Algebra in Haskell Motivation Ohm’s law R X , 1 , 1 R Z , 0 , 1 R Z , 1 , 1 R Y , 0 , 1 R Y , 1 , 1 R X , 1 , 0 R X , 0 , 1 R Y , 0 , 0 R Y , 1 , 0 R Z , 0 , 0 R X , 0 , 0 = U X , 0 , 0 I X , 0 , 0 , R Z , 1 , 0 R Y , 1 , 0 = U Y , 1 , 0 I Y , 1 , 0 , R X , 0 , 0 . . .

Expressive Linear Algebra in Haskell Motivation Kirchhoff’s current law I X , 1 , 1 I Z , 0 , 1 I Z , 1 , 1 I Y , 0 , 1 I Y , 1 , 1 I X , 1 , 0 I X , 0 , 1 I Y , 0 , 0 I Y , 1 , 0 I Z , 0 , 0 I Z , 1 , 0 0 = − I X , 0 , 0 + I Y , 1 , 0 + I Z , 1 , 0 , 0 = − I X , 0 , 1 + I Y , 1 , 1 − I Z , 1 , 0 , I X , 0 , 0 . . .

Expressive Linear Algebra in Haskell Motivation Kirchhoff’s voltage law V 0 , 1 , 1 U X , 1 , 1 V 1 , 1 , 1 U Z , 0 , 1 U Z , 1 , 1 U Y , 0 , 1 V 0 , 1 , 0 U Y , 1 , 1 V 1 , 1 , 0 U X , 1 , 0 U X , 0 , 1 V 0 , 0 , 1 U Y , 0 , 0 V 1 , 0 , 1 U Y , 1 , 0 U Z , 0 , 0 U Z , 1 , 0 U X , 0 , 0 = − V 0 , 0 , 0 + V 1 , 0 , 0 , V 0 , 0 , 0 U Y , 1 , 0 = − V 1 , 0 , 0 + V 1 , 1 , 0 , U X , 0 , 0 V 1 , 0 , 0 . . .

Expressive Linear Algebra in Haskell Motivation 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 I       0 R 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 − 1 0 0 0 I 0 0 0 R 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 − 1 0 0 I 0       0 0 0 R 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 − 1 0 I 0       0 0 0 0 R 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 − 1 I 0       0 0 0 0 0 R 0 0 0 0 0 0 0 1 0 − 1 0 0 0 0 0 I 0        0 0 0 0 0 0 R 0 0 0 0 0 0 0 1 0 − 1 0 0 0 0   I   0        0 0 0 0 0 0 0 R 0 0 0 0 0 0 0 0 0 1 0 − 1 0 I 0       0 0 0 0 0 0 0 0 R 0 0 0 0 0 0 0 0 0 1 0 − 1 I 0       0 0 0 0 0 0 0 0 0 R 0 0 0 1 − 1 0 0 0 0 0 0 I 0       = 0 0 0 0 0 0 0 0 0 0 R 0 0 0 0 1 − 1 0 0 0 0 I 0 ·       0 0 0 0 0 0 0 0 0 0 0 R 0 0 0 0 0 1 − 1 0 0 I 0       0 0 0 0 0 0 0 0 0 0 0 0 R 0 0 0 0 0 0 1 − 1 0    I          1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 V 0       0 0 1 0 0 0 1 0 0 − 1 0 0 0 0 0 0 0 0 0 0 0 V 0       0 0 0 1 0 − 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 V 0       0 0 0 0 1 0 − 1 0 0 0 − 1 0 0 0 0 0 0 0 0 0 0 V 0       0 − 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 V 0        0 0 − 1 0 0 0 0 0 1 0 0 − 1 0 0 0 0 0 0 0 0 0   V   0  0 0 0 − 1 0 0 0 − 1 0 0 0 0 1 0 0 0 0 0 0 0 0 V 0 0 0 0 0 − 1 0 0 0 − 1 0 0 0 − 1 0 0 0 0 0 0 0 0 V 1 A V 0 , 0 , 0 = 0 I total + I X , 0 , 0 + I Y , 0 , 0 + I Z , 0 , 0 = 0 R X , 0 , 0 · I X , 0 , 0 + V 0 , 0 , 0 − V 1 , 0 , 0 = 0 − I X , 1 , 1 − I Y , 1 , 1 − I Z , 1 , 1 = 1 A

Expressive Linear Algebra in Haskell Motivation Solution U total = (1) R total I total V 1 , 1 , 1 − V 0 , 0 , 0 = (2) I total

Expressive Linear Algebra in Haskell Motivation Matrix blocks   u T 0 0 0 , 0 , 0 A = 0 R V     u 0 , 0 , 0 I 0 R : encodes Ohm’s law V : encodes Kirchhoff’s voltage law I : encodes Kirchhoff’s current law

Expressive Linear Algebra in Haskell Motivation Matrix symmetry   u T 0 0 0 , 0 , 0 A = 0 R V     0 u 0 , 0 , 0 I   R X , 0 , 0 R X , 0 , 1     R = ...       R Z , 1 , 1 I = V T

Expressive Linear Algebra in Haskell Motivation Matrix features Matrix A is symmetric, composed from blocks, and every block has a special sub-structure. Can we represent this with Haskell’s type system?

Expressive Linear Algebra in Haskell Solution 1 Motivation 2 Solution 3 More realistic problem 4 More features 5 Closing

Expressive Linear Algebra in Haskell Solution Solution of simultaneous linear equations Specialised: (#\|) :: (Shape height , Eq Eq Eq height , Floating Floating Floating a) => => => Matrix.Symmetric height a -> Vector height a -> Vector height a Generic: (#\|) :: (Solve typ , HeightOf typ ˜ height , Eq Floating => Eq Eq height , Floating Floating a) => => Matrix typ a -> Vector height a -> Vector height a Infix operator reads as: “Matrix divides column vector”. Implemented by LAPACK’s SPTRS

Expressive Linear Algebra in Haskell Solution Array shapes and indices comfort-array Index type is a type function of the array shape. class class class Shape.C shape where where where Int size :: shape -> Int Int class => where class class Shape.C shape => => Shape.Indexed shape where where type type type Index shape Read a single element: Array (!) :: Array Array sh a -> Index sh -> a

Expressive Linear Algebra in Haskell Solution Zero-based indexing shape newtype newtype newtype Shape.ZeroBased n = ZeroBased n instance instance (Integral instance Integral Integral n) => => => Shape.Indexed (ZeroBased n) where where where type type type Index (ZeroBased n) = n Array Int Int (!) :: Array Array (ZeroBased Int Int) a -> Int Int -> a array array array ! 0 Classical zero-based indexing scheme as in hmatrix In contrast to array lower bound is statically fixed to zero

Expressive Linear Algebra in Haskell Solution Enumeration shape data data data Shape.Enumeration enum = Enumeration instance instance (Enum instance Enum Enum enum , Bounded Bounded Bounded enum) => => => Shape.Indexed (Enumeration enum) where where where type type type Index (Enumeration enum) = enum Array Ordering Ordering (!) :: Array Array (Enumeration Ordering Ordering) a -> Ordering Ordering -> a array compare array array ! compare compare x y Shape statically determined by an Enum Enum Enum type

Expressive Linear Algebra in Haskell Solution Cartesian product shape instance instance instance => (Shape.Indexed sha , Shape.Indexed shb) => => Shape.Indexed (sha , shb) where where where type type type Index (sha , shb) = (Index sha , Index shb) type type type E = Enumeration (!) :: Array Ordering Bool Ordering Bool Array Array (E Ordering Ordering , E Bool Bool) a -> (Ordering Ordering , Bool Bool) -> a array compare odd array array ! (compare compare x y, odd odd x) Represents a two-dimensional array of rectangular shape

Expressive Linear Algebra in Haskell Solution Sum shape instance instance instance (Shape.Indexed sha , Shape.Indexed shb) => => => Shape.Indexed (sha :+: shb) where where where type type type Index (sha :+: shb) = Either Either Either (Index sha) (Index shb) type type type E = Enumeration (!) :: Array Array Array (E Ordering Ordering Ordering :+: E Bool Bool Bool) a -> Either Ordering Bool Either Either Ordering Ordering Bool Bool -> a array Right False array array ! Right Right False False Useful for block matrices

Expressive Linear Algebra in Haskell Solution Array shapes for corners and edges data data data Coord = C0 | C1 deriving deriving deriving (Eq Eq Eq , Ord Ord Ord , Show Show Show , Enum Enum Enum , Bounded Bounded Bounded) data data data Dim = DX | DY | DZ deriving deriving deriving (Eq Eq Eq , Ord Ord Ord , Show Show Show , Enum Enum , Bounded Enum Bounded Bounded) type type type Corner = (Coord ,Coord ,Coord) type type type Edge = (Dim ,Coord ,Coord) type type type CoordSh = Shape.Enumeration Coord type type type DimSh = Shape.Enumeration Dim type type type CornerShape = (CoordSh , CoordSh , CoordSh) type type type EdgeShape = (DimSh , CoordSh , CoordSh) type type type BlockShape = ():+: EdgeShape :+: CornerShape

Expressive Linear Algebra in Haskell Solution Achievement no need to write index flattening function yourself consistent structure of matrix and vectors no fight over zero- vs. one-based index counting no off-by-one errors

Expressive Linear Algebra in Haskell Solution Voltage matrix voltageMatrix :: Matrix EdgeShape CornerShape a voltageMatrix = Matrix.fromRowArray cornerShape $ fmap (\e -> Array Array Array. fromAssociations cornerShape 0 [( edgeCorner e C0 , 1), (edgeCorner e C1 , -1)]) $ indices BoxedArray.indices indices edgeShape edgeCorner :: Edge -> Coord -> Corner edgeCorner (ed ,e0 ,e1) coord = case case case ed of of of DX -> (coord ,e0 ,e1) DY -> (e0 ,coord ,e1) DZ -> (e0 ,e1 ,coord)