diff --git a/numhask-array.cabal b/numhask-array.cabal
--- a/numhask-array.cabal
+++ b/numhask-array.cabal
@@ -1,6 +1,6 @@
 cabal-version:      2.4
 name:               numhask-array
-version:            0.10.0
+version:            0.10.1
 synopsis:           Multi-dimensional arrays.
 description:
   This package provides an interface into the [numhask](https://hackage.haskell.org/package/numhask) API, and both type- and value-level shape manipulation routines.
diff --git a/src/NumHask/Array/Fixed.hs b/src/NumHask/Array/Fixed.hs
--- a/src/NumHask/Array/Fixed.hs
+++ b/src/NumHask/Array/Fixed.hs
@@ -5,6 +5,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE OverloadedLists #-}
 {-# LANGUAGE PolyKinds #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE RebindableSyntax #-}
@@ -14,6 +15,7 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE NoStarIsType #-}
+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}
 {-# OPTIONS_GHC -Wno-redundant-constraints #-}
 {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
 {-# OPTIONS_GHC -fno-warn-unused-imports #-}
@@ -29,10 +31,13 @@
     toDynamic,
 
     -- * Operators
+    takes,
     reshape,
     transpose,
     diag,
+    undiag,
     ident,
+    sequent,
     singleton,
     selects,
     selectsExcept,
@@ -46,6 +51,7 @@
     append,
     reorder,
     expand,
+    expand',
     apply,
     contract,
     dot,
@@ -66,6 +72,7 @@
     --
     -- Vector specialisations.
     Vector,
+    sequentv,
 
     -- * Matrix
 
@@ -77,6 +84,9 @@
     safeCol,
     safeRow,
     mmult,
+    chol,
+    inv,
+    invtri,
   )
 where
 
@@ -90,7 +100,7 @@
 import GHC.TypeLits
 import qualified NumHask.Array.Dynamic as D
 import NumHask.Array.Shape
-import NumHask.Prelude as P hiding (toList)
+import NumHask.Prelude as P hiding (sequence, toList)
 
 -- $setup
 --
@@ -316,6 +326,17 @@
     isDiag [x, y] = x == y
     isDiag (x : y : xs) = x == y && isDiag (y : xs)
 
+-- | A sequential array of Ints
+--
+-- >>> sequent :: Array '[3] Int
+-- [0, 1, 2]
+sequent :: forall s. (HasShape s) => Array s Int
+sequent = tabulate go
+  where
+    go [] = zero
+    go [i] = i
+    go (i : j : _) = bool zero i (i == j)
+
 -- | Extract the diagonal of an array.
 --
 -- >>> diag (ident :: Array '[3,2] Int)
@@ -643,6 +664,20 @@
   where
     r = rank (shape a)
 
+expand' ::
+  forall s s' a b c.
+  ( HasShape s,
+    HasShape s',
+    HasShape ((++) s s')
+  ) =>
+  (a -> b -> c) ->
+  Array s a ->
+  Array s' b ->
+  Array ((++) s s') c
+expand' f a b = tabulate (\i -> f (index a (drop r i)) (index b (take r i)))
+  where
+    r = rank (shape a)
+
 -- | Apply an array of functions to each array of values.
 --
 -- This is in the spirit of the applicative functor operation (<*>).
@@ -858,6 +893,15 @@
     go s = index a (zipWith (!!) pss' s)
     pss' = natValss pss
 
+takes ::
+  forall s s' a.
+  ( HasShape s,
+    HasShape s'
+  ) =>
+  Array s a ->
+  Array s' a
+takes a = tabulate $ \s -> index a s
+
 -- | Remove single dimensions.
 --
 -- >>> let a = [1..24] :: Array '[2,1,3,4,1] Int
@@ -929,6 +973,9 @@
 -- | <https://en.wikipedia.org/wiki/Vector_(mathematics_and_physics) Wiki Vector>
 type Vector s a = Array '[s] a
 
+sequentv :: forall n. (KnownNat n) => Vector n Int
+sequentv = sequent
+
 -- | <https://en.wikipedia.org/wiki/Matrix_(mathematics) Wiki Matrix>
 type Matrix m n a = Array '[m, n] a
 
@@ -944,6 +991,73 @@
   (*) = mmult
 
   one = ident
+
+instance
+  ( Multiplicative a,
+    P.Distributive a,
+    Subtractive a,
+    Eq a,
+    ExpField a,
+    KnownNat m,
+    HasShape '[m, m]
+  ) =>
+  Divisive (Matrix m m a)
+  where
+  recip a = invtri (transpose (chol a)) * invtri (chol a)
+
+-- | inverse of a triangular matrix
+invtri :: forall a n. (KnownNat n, ExpField a, Eq a) => Array '[n, n] a -> Array '[n, n] a
+invtri a = sum (fmap (l ^) (sequentv :: Vector n Int)) * ti
+  where
+    ti = undiag (fmap recip (diag a))
+    tl = a - undiag (diag a)
+    l = negate (ti * tl)
+
+-- | Expand the array to form a diagonal array
+--
+-- >>> undiag ([1,1] :: Array '[2] Int)
+-- [[1, 0],
+--  [0, 1]]
+undiag ::
+  forall a s.
+  ( HasShape s,
+    Additive a,
+    HasShape ((++) s s)
+  ) =>
+  Array s a ->
+  Array ((++) s s) a
+undiag a = tabulate go
+  where
+    go [] = throw (NumHaskException "Rank Underflow")
+    go xs@(x : xs') = bool zero (index a xs) (all (x ==) xs')
+
+-- | cholesky decomposition
+chol :: (KnownNat n, ExpField a) => Array '[n, n] a -> Array '[n, n] a
+chol a =
+  let l =
+        tabulate
+          ( \[i, j] ->
+              bool
+                ( one
+                    / index l [j, j]
+                    * ( index a [i, j]
+                          - sum
+                            ( (\k -> index l [i, k] * index l [j, k])
+                                <$> ([zero .. (j - one)] :: [Int])
+                            )
+                      )
+                )
+                ( sqrt
+                    ( index a [i, i]
+                        - sum
+                          ( (\k -> index l [j, k] ^ 2)
+                              <$> ([zero .. (j - one)] :: [Int])
+                          )
+                    )
+                )
+                (i == j)
+          )
+   in l
 
 -- | Extract specialised to a matrix.
 --