Add complete interval-union Gauss–Seidel step method

brush-strokes/src/cusps/bench/Main.hs

@@ -69,7 +69,7 @@ benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOp
       ( dunno, sols ) =
         foldMap
           ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) ->
-              ( map ( ( i , ) . snd ) mbCusps, map ( i, ) defCusps ) ) $
+              ( map ( ( i , ) . fst ) mbCusps, map ( i, ) defCusps ) ) $
           IntMap.toList $
             findCuspsIn testCuspOptions testStrokeFnI $
               IntMap.fromList

brush-strokes/src/lib/Math/Linear.hs

@@ -128,8 +128,17 @@ pattern V4 x y z w = T ( ℝ4 x y z w )
 
 type Vec :: Nat -> Type -> Type
 newtype Vec n a = Vec { vecList :: [ a ] }
-  deriving newtype ( Show, Eq, Ord, Functor, Foldable )
-  deriving Applicative via ZipList
+type role Vec nominal representational
+
+deriving newtype instance Show a => Show ( Vec n a )
+deriving newtype instance   Eq a =>   Eq ( Vec n a )
+deriving newtype instance  Ord a =>  Ord ( Vec n a )
+deriving newtype instance Functor  ( Vec n )
+deriving newtype instance Foldable ( Vec n )
+deriving via ZipList
+  instance Applicative ( Vec n )
+instance Traversable ( Vec n ) where
+  traverse f ( Vec as ) = Vec <$> traverse f as
 
 universe :: forall n. KnownNat n => Vec n ( Fin n )
 universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ]

brush-strokes/src/lib/Math/Root/Isolation.hs

@@ -138,9 +138,9 @@ defaultRootIsolationAlgorithms minWidth ε_eq history box =
   where
     verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box
 
-    _bisOptions = defaultBisectionOptions minWidth ε_eq box
-    _gsOptions = defaultGaussSeidelOptions history
-    _box1Options = defaultBox1Options minWidth ε_eq
+    _bisOptions = defaultBisectionOptions @n @d minWidth ε_eq box
+    _gsOptions = defaultGaussSeidelOptions @n @d history
+    _box1Options = defaultBox1Options @n @d minWidth ε_eq
 
     -- Did we reduce the box width by at least ε_eq
     -- in at least one of the coordinates?

brush-strokes/src/lib/Math/Root/Isolation/Core.hs

@@ -89,7 +89,6 @@ type BoxCt n d =
   , Vars ( D 1 ( ℝ n ) ) ~ n
   , Module Double ( T ( ℝ n ) )
   , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
-  , Applicative ( Vec n )
 
   , Ord ( ℝ d )
   , Module Double ( T ( ℝ d ) )

brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs

@@ -12,6 +12,7 @@ module Math.Root.Isolation.GaussSeidel
   where
 
 -- base
+import Prelude hiding ( unzip )
 import Control.Arrow
   ( first )
 import Data.Bifunctor
@@ -24,6 +25,8 @@ import Data.Foldable
  ( toList )
 import Data.List
   ( partition )
+import Data.List.NonEmpty
+  ( unzip )
 import Data.Proxy
   ( Proxy(..) )
 import Data.Type.Ord
@@ -71,12 +74,21 @@ instance RootIsolationAlgorithm GaussSeidel where
 type GaussSeidelOptions :: Nat -> Nat -> Type
 data GaussSeidelOptions n d =
   GaussSeidelOptions
-    { -- | Which preconditioner to user?
+    { -- | Which preconditioner to use?
       gsPreconditioner :: !Preconditioner
       -- | Function that projects over the equations we will consider
       -- (the identity for a well-determined problem, or a projection for
       -- an overdetermined system).
-    , gsPickEqs :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
+    , gsPickEqs :: !( 𝕀ℝ d -> 𝕀ℝ n )
+      -- | Whether to use a partial or a complete Gauss–Seidel update
+    , gsUpdate :: !GaussSeidelUpdateMethod
+    }
+
+-- | Whether to use a partial or a complete Gauss–Seidel update.
+data GaussSeidelUpdateMethod
+  = GS_Partial
+  | GS_Complete
+  deriving stock Show
 
 -- | Default options for the interval Gauss–Seidel method.
 defaultGaussSeidelOptions
@@ -90,7 +102,7 @@ defaultGaussSeidelOptions
   -> GaussSeidelOptions n d
 defaultGaussSeidelOptions history =
   GaussSeidelOptions
-    { gsPreconditioner = InverseMidJacobian
+    { gsPreconditioner = InverseMidpoint
     , gsPickEqs =
         case cmpNat @n @d Proxy Proxy of
           EQI -> id
@@ -103,13 +115,14 @@ defaultGaussSeidelOptions history =
                 choice :: Vec n ( Fin d )
                 choice = choices !! ( length history `mod` length choices )
             in \ u -> tabulate \ i -> index u ( choice ! i )
+    , gsUpdate = GS_Complete
     }
 {-# INLINEABLE defaultGaussSeidelOptions #-}
 
 -- | Preconditioner to use with the interval Gauss–Seidel method.
 data Preconditioner
   = NoPreconditioning
-  | InverseMidJacobian
+  | InverseMidpoint
   deriving stock Show
 
 -- | Interval Newton method with Gauss–Seidel step.
@@ -123,44 +136,66 @@ intervalGaussSeidel
      -- ^ box
   -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ]
 intervalGaussSeidel
-  ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } )
+  ( GaussSeidelOptions
+    { gsPreconditioner = precondMeth
+    , gsPickEqs = pickEqs
+    , gsUpdate
+    } )
   eqs
-  box
-  | let boxMid = singleton $ boxMidpoint box
-        f' :: Vec n ( 𝕀ℝ n )
-        f' = fmap ( \ i -> pickEqs $ eqs box `monIndex` linearMonomial i ) ( universe @n )
-        f_mid = pickEqs $ eqs boxMid `monIndex` zeroMonomial
+  x
+  | let x_mid = singleton $ boxMidpoint x
+        f'_x :: Vec n ( 𝕀ℝ n )
+        f'_x = fmap ( \ i -> pickEqs $ eqs x `monIndex` linearMonomial i ) ( universe @n )
+        f_x_mid = pickEqs $ eqs x_mid `monIndex` zeroMonomial
 
   = let -- Interval Newton method: take one Gauss–Seidel step
         -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
+        minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint f_x_mid )
+
+        -- Precondition the above linear system into A ( x - x_mid ) = B.
         ( a, b ) = precondition precondMeth
-                      ( fmap boxMidpoint f' )
-                      f' ( singleton $ unT $ -1 *^ T ( boxMidpoint f_mid ) )
+                      f'_x ( singleton minus_f_x_mid )
 
         -- NB: we have to change coordinates, putting the midpoint of the box
         -- at the origin, in order to take a Gauss–Seidel step.
-        gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) )
-                  $ gaussSeidelStep a b ( T box ^-^ T boxMid )
+        gsGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) )
+                  $ gaussSeidelUpdate gsUpdate a b ( T x ^-^ T x_mid )
     in
          -- If the Gauss–Seidel step was a contraction, then the box
          -- contains a unique solution (by the Banach fixed point theorem).
          --
          -- These boxes can thus be directly added to the solution set:
          -- Newton's method is guaranteed to converge to the unique solution.
-         let !(done, todo) = bimap ( map fst ) ( map fst )
+         let ( done, todo ) = bimap ( map fst ) ( map fst )
                             $ partition snd gsGuesses
          in do tell $ noDoneBoxes { doneSolBoxes = done }
                return todo
   where
 {-# INLINEABLE intervalGaussSeidel #-}
 
+-- | A partial or complete Gauss–Seidel step for the equation \( A X = B \),
+-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
+gaussSeidelUpdate
+  :: forall n
+  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) )
+  => GaussSeidelUpdateMethod -- ^ which step method to use
+  -> Vec n ( 𝕀ℝ n )    -- ^ columns of \( A \)
+  -> 𝕀ℝ n              -- ^ \( B \)
+  -> T ( 𝕀ℝ n )        -- ^ initial box \( X \)
+  -> [ ( T ( 𝕀ℝ n ), Bool ) ]
+gaussSeidelUpdate upd as b x =
+  case upd of
+    GS_Partial  -> gaussSeidelStep          as b x
+    GS_Complete -> gaussSeidelStep_Complete as b x
+{-# INLINEABLE gaussSeidelUpdate #-}
+
 -- | Take one interval Gauss–Seidel step for the equation \( A X = B \),
 -- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
 --
 -- The boolean indicates whether the Gauss–Seidel step was a contraction.
 gaussSeidelStep
   :: forall n
-  .  ( Representable Double ( ℝ n ), Eq ( ℝ n ) )
+  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) )
   => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
   -> 𝕀ℝ n           -- ^ \( B \)
   -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
@@ -168,15 +203,80 @@ gaussSeidelStep
 gaussSeidelStep as b ( T x0 ) = coerce $
   forEachCoord @n ( x0, True ) $ \ i ( x, contraction ) -> do
     -- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
-    x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )
-          `extendedDivide` ( ( as ! i ) `index` i )
-    ( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i )
+    let s = b `index` i - sum [ ( as ! j ) `index` i * x `index` j
+                              | j <- toList ( universe @n ), j /= i ]
+        x_i = x `index` i
+        a_ii = ( as ! i ) `index` i
+    -- Take a shortcut before performing the division if possible.
+    if | not $ 0 `inside` ( s - a_ii * x_i )
+       -- No solutions: don't bother performing a division.
+       -> [ ]
+       | 0 `inside` s && 0 `inside` a_ii
+       -- The division would produce [-oo,+oo]: don't do anything.
+       -> [ ( x, False ) ]
+       -- Otherwise, perform the division.
+       | otherwise
+       -> do
+           x_i'0 <- s `extendedDivide` a_ii
+           ( x_i', sub_i ) <- x_i'0 `intersect` x_i
            return $ ( set i x_i' x, sub_i && contraction )
--- TODO: try implementing the complete interval union Gauss–Seidel algorithm.
--- See "Algorithm 2" in
--- "Using interval unions to solve linear systems of equations with uncertainties"
 {-# INLINEABLE gaussSeidelStep #-}
 
+-- | The complete interval-union Gauss–Seidel step.
+--
+-- Algorithm 2 from:
+--   "Using interval unions to solve linear systems of equations with uncertainties"
+--    (Montanher, Domes, Schichl, Neumaier) (2017)
+gaussSeidelStep_Complete
+  :: forall n
+  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) )
+  => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
+  -> 𝕀ℝ n           -- ^ \( B \)
+  -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
+  -> [ ( T ( 𝕀ℝ n ), Bool ) ]
+gaussSeidelStep_Complete as b ( T x0 ) = coerce $ do
+  ( x', subs ) <-
+    forEachCoord @n ( x0, pure False ) $ \ i ( x, contractions ) -> do
+      let s = b `index` i - sum [ ( as ! k ) `index` i * x `index` k
+                                | k <- toList ( universe @n ) ]
+          x_i = x `index` i
+          a_ii = ( as ! i ) `index` i
+      ( x', subs ) <- fromComponents \ j -> do
+        let x_j = x `index` j
+            a_ij = ( as ! j ) `index` i
+            s_j = s `ominus` ( a_ij * x_j )
+        -- Shortcut division if possible (see gaussSeidelStep for commentary).
+        if | not $ 0 `inside` ( s_j - a_ii * x_i )
+           -> [ ]
+           | 0 `inside` s_j && 0 `inside` a_ij
+           -> [ ( x_j, False ) ]
+           | otherwise
+           -> do
+               x_j'0 <- s_j `extendedDivide` a_ij
+               ( x_j', sub_j ) <- x_j'0 `intersect` x_j
+               return $ ( x_j', sub_j )
+      return ( x', (||) <$> contractions <*> subs )
+  return ( x', and subs )
+{-# INLINEABLE gaussSeidelStep_Complete #-}
+
+fromComponents
+  :: forall n
+  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) )
+  => ( Fin n -> [ ( 𝕀 Double, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ]
+fromComponents f = do
+  ( xs, bs ) <- unzip <$> traverse f ( universe @n )
+  return $ ( tabulate $ \ i -> xs ! i, bs )
+    -- TODO: this could be more efficient.
+{-# INLINEABLE fromComponents #-}
+
+infixl 6 `ominus`
+ominus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double
+ominus a@( 𝕀 lo1 hi1 ) b@( 𝕀 lo2 hi2 )
+  | width a >= width b
+  = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 )
+  | otherwise
+  = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 )
+
 -- | The midpoint of a box.
 boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n
 boxMidpoint box =
@@ -191,16 +291,15 @@ precondition
   :: forall n
   .  ( KnownNat n, Representable Double ( ℝ n ) )
   => Preconditioner -- ^ pre-conditioning method to use
-  -> Vec n ( ℝ n )  -- ^ entry-wise midpoint matrix of the interval Jacobian matrix
   -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
   -> 𝕀ℝ n           -- ^ \( B \)
   -> ( Vec n ( 𝕀ℝ n ), 𝕀ℝ n )
-precondition meth jac_mid as b =
+precondition meth as b =
   case meth of
     NoPreconditioning
       -> ( as, b )
-    InverseMidJacobian
-      | let mat = toEigen jac_mid
+    InverseMidpoint
+      | let mat = toEigen $ fmap boxMidpoint as
             det = Eigen.determinant mat
       , not $ nearZero det
       -- (TODO: a bit wasteful to compute determinant then inverse.)