% $Id: SCG.lhs,v 1.6 1997/06/26 16:39:06 tjchol01 Exp tjchol01 $

\begin{code}
module SCG (
  F, DF, ScgInput (ScgInput),
  scg,
  -- for debugging only
  ScgData (ScgData), scgFull, showsScg
) where
import Arr
\end{code}

\begin{code}
type F a = Vector a -> a
type DF a = Vector a -> Vector a
\end{code}

\begin{code}
data ScgData a = ScgData {k :: Int, err :: a,
                          w, p, r :: Vector a,
                          delta, pnorm2, lambda, lambdabar :: a,
                          success :: Bool}
\end{code}

\begin{code}
calculate2order :: Floating a => ScgData a -> a -> DF a -> ScgData a
calculate2order d sigma1 df =
  let pnorm2' = vectorNorm2 (p d)
      sigma = sigma1 / (sqrt pnorm2')
      s = scaleVector (1/sigma) (df ((w d) + (scaleVector sigma (p d))) - df (w d))
  in d {pnorm2 = pnorm2', delta = inner (p d) s}
\end{code}

\begin{code}
hessPosDef :: (Floating a, Ord a) => ScgData a -> ScgData a
hessPosDef d =
  let delta' = delta d + (lambda d - lambdabar d) * pnorm2 d {- step 3 -}
  in if delta' <= 0                             {- step 4 -}
     then let lambdabar' = 2.0 * (lambda d - delta' / pnorm2 d)
          in d {delta = -delta' + lambda d * pnorm2 d, lambda = lambdabar', lambdabar = lambdabar'}
     else d {delta = delta'}
\end{code}

\begin{code}
reduceError :: (Floating a, Ord a) => ScgData a -> DF a -> Bool -> a -> a -> ScgData a
reduceError d df restart bdelta mu =
  let r' = negate (df (w d))
      p' = if restart
           then r'
           else let beta = (vectorNorm2 r' - inner r' (r d)) / mu
                in r' + scaleVector beta (p d)
  in d {p = p', r = r', lambda = if bdelta >= 0.75 then lambda d / 4 else lambda d
    }
\end{code}

\begin{code}
data ScgInput a = ScgInput (F a) (DF a) (Vector a)
\end{code}

\begin{code}
scgIter :: (Floating a, Ord a) => ScgInput a -> [ScgData a]
scgIter (ScgInput f df w1) =
    let p1 = negate (df w1)                     {- step 1 -}
        r1 = p1
        pnorm21 = vectorNorm2 p1
        n = arraySize w1
        sigma1 = 1.0e-4
        lambda1 = 1.0e-6
        err1 = f w1
    in iterate (\d ->
           let d1 = if success d                {- step 2 -}
                    then calculate2order d sigma1 df
                    else d
               d2 = hessPosDef d1
               mu = inner (p d2) (r d2)         {- step 5 -}
               alpha = mu / (delta d2)
               w' = (w d2) + scaleVector alpha (p d2)
               err' = f w'
               bdelta = 2 * (delta d2) * ((err d2) - err') / (mu^2) {- step 6 -}
               success' = (bdelta >= 0)         {- step 7 -}
               restart = ((k d) `mod` n == 0)
               d3 = if success'
                    then (reduceError (d2 {w = w'}) df restart bdelta mu)
                            {err = err', lambdabar = 0}
                    else d2 {lambdabar = lambda d2}
               d4 = if bdelta < 0.25          {- step 8 -}
                    then d3 {lambda = (lambda d3) + (delta d3) * (1 - bdelta) / (pnorm2 d3)}
                    else d3
           in d4 {k = k d4 + 1, success = success'}
       )
       (ScgData 1 err1 w1 p1 r1 0.0 pnorm21 lambda1 0.0 True)
\end{code}

Display internal SCG progress state
\begin{code}
showsScg :: (Show a) => ScgData a -> ShowS
showsScg (ScgData {k, err, w, p, r, delta, pnorm2, lambda, lambdabar, success}) =
  showString "Iter: " . shows k .
  showString ", err: " . shows err .
  showString ", w: " . shows w .
  showString ", p: " . shows p .
  showString ", r: " . shows r .
  showString ", delta: " . shows delta .
  showString ", pnorm2: " . shows pnorm2 .
  showString ", lambda: " . shows lambda .
  showString ", lambdabar: " . shows lambdabar .
  showString ", success: " . shows success .
  showString "\n"
\end{code}

Return data of the first iteration satisfying termination criteria
\begin{code}
scgFull :: (Floating a, Ord a) => ScgInput a -> a -> ScgData a
scgFull d eps = terminate (scgIter d) 0 where
  terminate (x:x':xs) underTol =
    if vectorNorm2 (r x') <= eps^2 || underTol > 2
    then x'
    else if 2 * abs (err x - err x') <= eps * abs (err x + err x' + eps) --!!!
         then terminate (x':xs) (underTol + 1)
         else if success x' -- success resets underTol
              then terminate (x':xs) 0
              else terminate (x':xs) underTol
\end{code}

Just the weights
\begin{code}
scg :: (Floating a, Ord a) => ScgInput a -> a -> Vector a
scg d eps = w (scgFull d eps)
\end{code}

Test cases:
\begin{itemize}
\item Rosenbrock valley $f(1,1) = 0$
\begin{code}
rosenbrock = ScgInput
  (\x -> 100 * (x!2 - x!1^2)^2 + (1 - x!1)^2)
  (\x -> listVector1 [-2 * (1 - x!1) - 400 * x!1 * (x!2 - x!1^2),
                              200 * (x!2 -x!1^2)])
  (listVector1 [-1,-1.0])
\end{code}

\item Hump camel $f(-1.7470015, 0.8670916) = -1.6526237$
\begin{code}
humpCamel = ScgInput
  (\x -> 4.0*x!1^2 - 2.1*x!1^4 + (x!1^6)/3.0 + 2*x!1*x!2 - 4*x!2^2 + 4*x!2^4)
  (\x -> listVector1 [8*x!1 - 8.4*x!1^3 + 2*x!1^5 + 2*x!2,
                              2*x!1 - 8*x!2 + 16*x!2^3])
  (listVector1 [1,1.0])
\end{code}

\item n-dimensional square error
\begin{code}
squareError n = ScgInput
  (sumVector . (map (^2)))
  (map (2*))
  (listVector1 (replicate n 1))
\end{code}

\item Levy function
\begin{code}
levy = ScgInput
  (\x -> sum [i * cos ((i+1) * x!1 + i) | i <- [1..5]] *
                   sum [i * cos ((i+1) * x!2 + i) | i <- [1..6]])
  (\x -> listVector1 [-sum [i * sin ((i+1) * x!1 + i) | i <- [1..5]] *
                                sum [i * cos ((i+1) * x!2 + i) | i <- [1..6]],
                               -sum [i * cos ((i+1) * x!1 + i) | i <- [1..5]] *
                                sum [i * sin ((i+1) * x!2 + i) | i <- [1..6]]])
  (listVector1 [1,1])
\end{code}
\end{itemize}
\begin{code}
--main = print (hdf (listVector1 [-1.74700153544529, 0.86709162557415]))
--main = print ((listVector1 [1,1]))
--main = sequence (map (\s -> putStr (showsScg s "")) (take 400 (scgIter rosenbrock)))
--main = putStr (showsScg (scgFull levy 1.0e-8) "")
\end{code}

hugs: tuples defined only up to 5 elements
hugs: -1^2 = 1

Based on:
@Article{Moller93,
  author = 	 "Martin F. M\o{}ller",
  title = 	 "A Scaled Conjugate Gradient Algorithm for Fast
		  Supervised Learning",
  journal =	 "Neural Networks",
  year =	 1993,
  volume =	 6,
  pages =	 "525-533"
}