src/Linear/Iteration.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2 {-# LANGUAGE TypeFamilies #-}
   3
   4 -- | Classical iterative methods to solve the system Ax = b.
   5
   6 module Linear.Iteration
   7 where
   8
   9 import Data.List (find)
  10 import Data.Maybe (fromJust)
  11 import Data.Vector.Fixed (Arity, N1, S)
  12 import NumericPrelude hiding ((*))
  13 import qualified Algebra.Algebraic as Algebraic
  14 import qualified Algebra.Field as Field
  15 import qualified Algebra.RealField as RealField
  16 import qualified Algebra.ToRational as ToRational
  17 import qualified Prelude as P
  18
  19 import Linear.Matrix
  20 import Linear.System
  21 import Normed
  22
  23 -- | A generalized implementation for Jacobi, Gauss-Seidel, etc. All
  24 --   that we really need to know is how to construct the matrix M, so we
  25 --   take a function that does it as an argument.
  26 classical_iteration :: (Field.C a, Arity m)
  27                  => (Mat m m  a -> Mat m m  a)
  28                  -> Mat m m  a
  29                  -> Mat m N1 a
  30                  -> Mat m N1 a
  31                  -> Mat m N1 a
  32 classical_iteration m_function matrix b x_current =
  33   x_next
  34   where
  35     big_m = m_function matrix
  36     big_n = big_m - matrix
  37     rhs = big_n*x_current + b
  38     -- TODO: Should be solve below! M might not be lower-triangular.
  39     x_next = forward_substitute big_m rhs
  40
  41
  42 -- | Perform one iteration of successive over-relaxation.
  43 --
  44 sor_iteration :: forall m a.
  45                  (Field.C a, Arity m)
  46               => a -- ^ Omega
  47               -> Mat m m  a -- ^ Matrix A
  48               -> Mat m N1 a -- ^ Vector b
  49               -> Mat m N1 a -- ^ Vector x_current
  50               -> Mat m N1 a -- ^ Output vector x_next
  51 sor_iteration omega =
  52   classical_iteration m_function
  53   where
  54     m_function :: Mat m m a -> Mat m m a
  55     m_function matrix =
  56       let diag = (recip omega) *> (diagonal_part matrix)
  57           lt   = lt_part_strict matrix
  58       in
  59         diag + lt
  60
  61
  62 -- | Compute an infinite list of SOR iterations starting with the
  63 --   vector x0.
  64 sor_iterations :: (Field.C a, Arity m)
  65                => a
  66                -> Mat m m  a
  67                -> Mat m N1 a
  68                -> Mat m N1 a
  69                -> [Mat m N1 a]
  70 sor_iterations omega matrix b =
  71   iterate (sor_iteration omega matrix b)
  72
  73
  74 -- | Perform one iteration of Gauss-Seidel.
  75 gauss_seidel_iteration :: (Field.C a, Arity m)
  76                        => Mat m m  a
  77                        -> Mat m N1 a
  78                        -> Mat m N1 a
  79                        -> Mat m N1 a
  80 gauss_seidel_iteration = sor_iteration one
  81
  82
  83 -- | Compute an infinite list of Gauss-Seidel iterations starting with
  84 --   the vector x0.
  85 gauss_seidel_iterations :: (Field.C a, Arity m)
  86                         => Mat m m  a
  87                         -> Mat m N1 a
  88                         -> Mat m N1 a
  89                         -> [Mat m N1 a]
  90 gauss_seidel_iterations matrix b =
  91   iterate (gauss_seidel_iteration matrix b)
  92
  93
  94 -- | Perform one Jacobi iteration,
  95 --
  96 --     x1 = M^(-1) * (N*x0 + b)
  97 --
  98 --   Examples:
  99 --
 100 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 101 --   >>> let x0 = vec2d (0, 0::Double)
 102 --   >>> let b  = vec2d (1, 1::Double)
 103 --   >>> jacobi_iteration m b x0
 104 --   ((0.25),(0.5))
 105 --   >>> let x1 = jacobi_iteration m b x0
 106 --   >>> jacobi_iteration m b x1
 107 --   ((0.0),(0.25))
 108 --
 109 jacobi_iteration :: (Field.C a, Arity m)
 110                  => Mat m m  a
 111                  -> Mat m N1 a
 112                  -> Mat m N1 a
 113                  -> Mat m N1 a
 114 jacobi_iteration =
 115   classical_iteration diagonal_part
 116
 117
 118 -- | Compute an infinite list of Jacobi iterations starting with the
 119 --   vector x0.
 120 jacobi_iterations :: (Field.C a, Arity m)
 121                   => Mat m m  a
 122                   -> Mat m N1 a
 123                   -> Mat m N1 a
 124                   -> [Mat m N1 a]
 125 jacobi_iterations matrix b =
 126   iterate (jacobi_iteration matrix b)
 127
 128
 129 -- | Solve the system Ax = b using the Jacobi method. This will run
 130 --   forever if the iterations do not converge.
 131 --
 132 --   Examples:
 133 --
 134 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 135 --   >>> let x0 = vec2d (0, 0::Double)
 136 --   >>> let b  = vec2d (1, 1::Double)
 137 --   >>> let epsilon = 10**(-6)
 138 --   >>> jacobi_method m b x0 epsilon
 139 --   ((0.0),(0.4999995231628418))
 140 --
 141 jacobi_method :: (RealField.C a,
 142                   Algebraic.C a, -- Normed instance
 143                   ToRational.C a, -- Normed instance
 144                   Algebraic.C b,
 145                   RealField.C b,
 146                   Arity m,
 147                   Arity n, -- Normed instance
 148                   m ~ S n)
 149               => Mat m m  a
 150               -> Mat m N1 a
 151               -> Mat m N1 a
 152               -> b
 153               -> Mat m N1 a
 154 jacobi_method =
 155   classical_method jacobi_iterations
 156
 157
 158 -- | Solve the system Ax = b using the Gauss-Seidel method. This will
 159 --   run forever if the iterations do not converge.
 160 --
 161 --   Examples:
 162 --
 163 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 164 --   >>> let x0 = vec2d (0, 0::Double)
 165 --   >>> let b  = vec2d (1, 1::Double)
 166 --   >>> let epsilon = 10**(-12)
 167 --   >>> gauss_seidel_method m b x0 epsilon
 168 --   ((4.547473508864641e-13),(0.49999999999954525))
 169 --
 170 gauss_seidel_method :: (RealField.C a,
 171                         Algebraic.C a, -- Normed instance
 172                         ToRational.C a, -- Normed instance
 173                         Algebraic.C b,
 174                         RealField.C b,
 175                         Arity m,
 176                         Arity n, -- Normed instance
 177                         m ~ S n)
 178                     => Mat m m  a
 179                     -> Mat m N1 a
 180                     -> Mat m N1 a
 181                     -> b
 182                     -> Mat m N1 a
 183 gauss_seidel_method =
 184   classical_method gauss_seidel_iterations
 185
 186
 187 -- | Solve the system Ax = b using the Successive Over-Relaxation
 188 --   (SOR) method. This will run forever if the iterations do not
 189 --   converge.
 190 --
 191 --   Examples:
 192 --
 193 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 194 --   >>> let x0 = vec2d (0, 0::Double)
 195 --   >>> let b  = vec2d (1, 1::Double)
 196 --   >>> let epsilon = 10**(-12)
 197 --   >>> sor_method 1.5 m b x0 epsilon
 198 --   ((6.567246746413957e-13),(0.4999999999993727))
 199 --
 200 sor_method :: (RealField.C a,
 201                Algebraic.C a, -- Normed instance
 202                ToRational.C a, -- Normed instance
 203                Algebraic.C b,
 204                RealField.C b,
 205                Arity m,
 206                Arity n, -- Normed instance
 207                m ~ S n)
 208            => a
 209            -> Mat m m  a
 210            -> Mat m N1 a
 211            -> Mat m N1 a
 212            -> b
 213            -> Mat m N1 a
 214 sor_method omega =
 215   classical_method (sor_iterations omega)
 216
 217
 218 -- | General implementation for all classical iteration methods. For
 219 --   its first argument, it takes a function which generates the
 220 --   sequence of iterates when supplied with the remaining arguments
 221 --   (except for the tolerance).
 222 --
 223 classical_method :: forall m n a b.
 224                     (RealField.C a,
 225                      Algebraic.C a, -- Normed instance
 226                      ToRational.C a, -- Normed instance
 227                      Algebraic.C b,
 228                      RealField.C b,
 229                      Arity m,
 230                      Arity n, -- Normed instance
 231                      m ~ S n)
 232                  => (Mat m m  a -> Mat m N1 a -> Mat m N1 a -> [Mat m N1 a])
 233                  -> Mat m m  a
 234                  -> Mat m N1 a
 235                  -> Mat m N1 a
 236                  -> b
 237                  -> Mat m N1 a
 238 classical_method iterations_function matrix b x0 epsilon =
 239   -- fromJust is "safe," because the list is infinite. If the
 240   -- algorithm doesn't converge, 'find' will search forever and never
 241   -- return Nothing.
 242   fst' $ fromJust $ find error_small_enough diff_pairs
 243   where
 244     x_n = iterations_function matrix b x0
 245
 246     pairs :: [(Mat m N1 a, Mat m N1 a)]
 247     pairs = zip (tail x_n) x_n
 248
 249     append_diff :: (Mat m N1 a, Mat m N1 a)
 250                 -> (Mat m N1 a, Mat m N1 a, b)
 251     append_diff (cur,prev) =
 252       (cur,prev,diff)
 253       where
 254         diff = norm (cur - prev)
 255
 256     diff_pairs :: [(Mat m N1 a, Mat m N1 a, b)]
 257     diff_pairs = map append_diff pairs
 258
 259     fst' :: (c, d, e) -> c
 260     fst' (x,_,_) = x
 261
 262     error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
 263     error_small_enough (_,_,err) = err < epsilon
 264
 265
 266
 267 -- | Compute the Rayleigh quotient of  @matrix@ and @vector@.
 268 --
 269 --   Examples:
 270 --
 271 --   >>> let m = fromList [[3,1],[1,2]] :: Mat2 Rational
 272 --   >>> let v = vec2d (1, 1::Rational)
 273 --   >>> rayleigh_quotient m v
 274 --   7 % 2
 275 --
 276 rayleigh_quotient :: (RealField.C a,
 277                       Arity m,
 278                       Arity n,
 279                       m ~ S n)
 280                   => (Mat m m a)
 281                   -> (Mat m N1 a)
 282                   -> a
 283 rayleigh_quotient matrix vector =
 284   (vector `dot` (matrix * vector)) / (norm_squared vector)
 285   where
 286     -- We don't use the norm function here to avoid the algebraic
 287     -- requirement on our field.
 288     norm_squared v = ((transpose v) * v) !!! (0,0)