# 
# L-BFGS-SKIP example with M=1
# 
# Written by 
# Tamara Gibson
# Applied Math Program, University of Maryland, College Park, MD 20742
# gibson@math.umd.edu
# 
# Companion to 
# Tamara Gibson, Dianne P. O'Leary and Larry Nazareth
# "BFGS with Update Skipping and Varying Memory"
# Technical Report CS-TR-3663 UMIACS-TR-96-49
# University of Maryland College Park, 1996
# 
# 1996
# 
--------------------------------------------------------------------------------
> \
# Initialize Problem\
\
b := [1,1,1];\
a := array([[1,0,0],[0,2,0],[0,0,4]]);\
x0 := [0,0,0];\
H0 := array([[1,0,0],[0,1,0],[0,0,1]]);\
\


                                 b := [1, 1, 1]

                                     [ 1  0  0 ]
                                     [         ]
                                a := [ 0  2  0 ]
                                     [         ]
                                     [ 0  0  4 ]

                                 x0 := [0, 0, 0]

                                      [ 1  0  0 ]
                                      [         ]
                                H0 := [ 0  1  0 ]
                                      [         ]
                                      [ 0  0  1 ]
--------------------------------------------------------------------------------
> \
# Initialize Algorithm\
\
x := evalm(x0):\
H := evalm(H0):\
g := evalm(evalm(a &* x) - b) ;\
Eye := array([[1,0,0],[0,1,0],[0,0,1]]);\
\


                               g := [ -1, -1, -1 ]

                                      [ 1  0  0 ]
                                      [         ]
                               Eye := [ 0  1  0 ]
                                      [         ]
                                      [ 0  0  1 ]
--------------------------------------------------------------------------------
> \
# L-BFGS-SKIP Iterations\
\
for k from 0 to 20 do\
  print(k);\
    # Compute Search Direction d_k\
  d := evalm(-H &* g):\
  # Line Search\
  alpha := evalm(-evalm(transpose(d) * g) / evalm(transpose(d) &* a &* d)):\
  # Compute s_k\
  s := evalm(alpha * d):\
  # Compute x_{k+1}\
  x := evalm (x + s):\
  # Compute y_k\
  y := evalm(a &* s):\
  # Compute g_{k+1}\
  g := evalm(g + y):\
  print(g);\
  # Compute H_{k+1}\
  if irem(k,2) = 1 then\
    # Skip Updates on Even Iterations\
    delta := evalm(transpose(y) * s):\
    P := evalm(Eye - evalm(y * transpose(s))/delta):\
    H := evalm(evalm(transpose(P) &* H0 &* P) + evalm(s*transpose(s))/delta):\
   fi:\
 od:\


                                        0

                               [ -4/7, -1/7, 5/7 ]

                                        1

                                 152     17    125
                             [ - ---, - ---, - --- ]
                                 413    413    413

                                        2

                                  75    225   15
                             [ - ----, ----, ---- ]
                                 6664  3332  3332

                                        3

                             8550       225      28125
                       [ - -------, - ------, - ------- ]
                           1533553    360836    6134212

                                        4

                               225       225    1125
                          [ - -----, - ------, ------ ]
                              75803    303212  303212

                                        5

                            8550       225        28125
                      [ - -------, - -------, - -------- ]
                          4472377    1052324    17889508

                                        6

                            16875      50625       3375
                      [ - ---------, ---------, --------- ]
                          288657824  144328912  144328912

                                        7

                       961875         50625          6328125
                [ - -----------, - -----------, - ------------ ]
                    33213690874    15629972176    265709526992

                                        8

                         50625         50625        253125
                  [ - ----------, - -----------, ----------- ]
                      3283482748    13133930992  13133930992

                                        9

                       961875         50625          6328125
                [ - -----------, - -----------, - ------------ ]
                    96862741066    45582466384    774901928528

                                       10

                       3796875        11390625        759375
               [ - --------------, -------------, ------------- ]
                   12503502304384  6251751152192  6251751152192

                                       11

                 216421875           11390625           1423828125
         [ - ----------------, - ---------------, - ----------------- ]
             1438684233898184    677027874775616    11509473871185472

                                       12

                    11390625           11390625         56953125
            [ - ---------------, - ---------------, --------------- ]
                142227338712368    568909354849472  568909354849472

                                       13

                 216421875           11390625            1423828125
         [ - ----------------, - ----------------, - ----------------- ]
             4195706492014856    1974450113889344    33565651936118848

                                       14

                 854296875          2562890625           170859375
        [ - ------------------, ------------------, ------------------ ]
            541601705816697344  270800852908348672  270800852908348672

                                       15

            48694921875             2562890625              320361328125
   [ - --------------------, - --------------------, - --------------------- ]
       62318046275533738144    29326139423780582656    498544370204269905152

                                       16

              2562890625             2562890625            12814453125
     [ - -------------------, - --------------------, -------------------- ]
         6160719403664932288    24642877614659729152  24642877614659729152

                                       17

           48694921875              2562890625              320361328125
  [ - ---------------------, - --------------------, - ---------------------- ]
      181741222408115502496    85525281133230824704    1453929779264924019968

                                       18

           192216796875             576650390625             38443359375
 [ - -----------------------, -----------------------, ----------------------- ]
     23460019489156062152704  11730009744578031076352  11730009744578031076352

                                       19

                    10956357421875                576650390625
          [ - -------------------------, - -------------------------,
              2699368492471019401445504    1270291055280479718327296

                      72081298828125
              - -------------------------- ]
                21594947939768155211564032

                                       20

                     576650390625                 576650390625
           [ - ------------------------, - -------------------------,
               266857721689150206987008    1067430886756600827948032

                     2883251953125
               ------------------------- ]
               1067430886756600827948032
