First version of LBFGS operators.

src/LinearOperators.jl

@@ -9,6 +9,7 @@ using Docile
 
 export LinearOperator, opEye, opOnes, opZeros, opDiagonal,
        opInverse, opCholesky, opHouseholder, opHermitian,
+       LBFGSOperator, InverseLBFGSOperator,
        check_ctranspose, check_hermitian, check_positive_definite,
        shape, hermitian, symmetric
 
@@ -474,5 +475,175 @@ function opHermitian(T :: KindOfMatrix)
   return opHermitian(d, T);
 end
 
+
+@doc "A data type to hold information relative to LBFGS operators." ->
+type LBFGSData
+  mem :: Int;
+  scaling :: Bool;
+  s   :: Array;
+  y   :: Array;
+  ys  :: Vector;
+  α   :: Vector;
+  a   :: Array;
+  b   :: Array;
+  insert :: Int;
+
+  function LBFGSData(n :: Int, mem :: Int;
+                     dtype :: DataType=Float64, scaling :: Bool=false, inverse :: Bool=true)
+    return new(max(mem, 1),
+               scaling,
+               zeros(dtype, n, mem),
+               zeros(dtype, n, mem),
+               zeros(dtype, mem),
+               zeros(dtype, mem),
+               inverse ? dtype[] : zeros(dtype, n, mem),
+               inverse ? dtype[] : zeros(dtype, n, mem),
+               1)
+  end
+end
+
+
+@doc "A type for limited-memory BFGS approximations." ->
+type LBFGSOperator <: AbstractLinearOperator
+  nrow   :: Int
+  ncol   :: Int
+  dtype   :: DataType
+  symmetric :: Bool
+  hermitian :: Bool
+  prod   :: Function           # apply the operator to a vector
+  tprod  :: Nullable{Function} # apply the transpose operator to a vector
+  ctprod :: Nullable{Function} # apply the transpose conjugate operator to a vector
+  inverse :: Bool
+  data :: LBFGSData
+end
+
+@doc "Construct a limited-memory BFGS approximation in inverse form." ->
+function InverseLBFGSOperator(n, mem :: Int=5; dtype :: DataType=Float64, scaling :: Bool=false)
+  lbfgs_data = LBFGSData(n, mem, dtype=dtype, scaling=scaling);
+  insert = 1;
+
+  function lbfgs_multiply(data :: LBFGSData, x :: Array)
+    # Multiply operator with a vector.
+    # See, e.g., Nocedal & Wright, 2nd ed., Procedure 7.4, p. 178.
+
+    if dtype == typeof(x[1])
+      q = copy(x);
+    else
+      result_type = promote_type(dtype, typeof(x[1]))
+      q = convert(Array{result_type}, x);
+    end
+
+    for i = 1 : data.mem
+      k = mod(data.insert - i - 1, data.mem) + 1;
+      if data.ys[k] != 0
+        data.α[k] = dot(data.s[:,k], q) / data.ys[k];
+        q -= data.α[k] * data.y[:,k];
+      end
+    end
+
+    r = q;
+    if data.scaling
+      last = mod(data.insert -1, data.mem) + 1;
+      if data.ys[last] != 0
+        γ = data.ys[last] / dot(data.y[:,last], data.y[:,last]);
+        r *= γ
+      end
+    end
+
+    for i = 1 : data.mem
+      k = mod(data.insert + i - 2, data.mem) + 1;
+      if data.ys[k] != 0
+        β = dot(data.y[:,k], r) / data.ys[k];
+        r += (data.α[k] - β) * data.s[:,k];
+      end
+    end
+
+    return r
+  end
+
+  return LBFGSOperator(n, n, dtype, true, true,
+                       x -> lbfgs_multiply(lbfgs_data, x),
+                       Nullable{Function}(),
+                       Nullable{Function}(),
+                       true,
+                       lbfgs_data)
+end
+
+@doc "Construct a limited-memory BFGS approximation in forward form." ->
+function LBFGSOperator(n, mem :: Int=5; dtype :: DataType=Float64, scaling :: Bool=false)
+  lbfgs_data = LBFGSData(n, mem, dtype=dtype, scaling=scaling, inverse=false);
+  insert = 1;
+
+  function lbfgs_multiply(data :: LBFGSData, x :: Array)
+    # Multiply operator with a vector.
+    # See, e.g., Nocedal & Wright, 2nd ed., Procedure 7.6, p. 184.
+
+    if dtype == typeof(x[1])
+      q = copy(x);
+    else
+      result_type = promote_type(dtype, typeof(x[1]))
+      q = convert(Array{result_type}, x);
+    end
+
+    # B = B₀ + Σᵢ (bᵢbᵢ' - aᵢaᵢ').
+    for i = 1 : data.mem
+      k = mod(data.insert + i - 2, data.mem) + 1;
+      if data.ys[k] != 0
+        q += dot(data.b[:, k], x) * data.b[:, k] - dot(data.a[:, k], x) * data.a[:, k];
+      end
+    end
+    return q
+  end
+
+  return LBFGSOperator(n, n, dtype, true, true,
+                       x -> lbfgs_multiply(lbfgs_data, x),
+                       Nullable{Function}(),
+                       Nullable{Function}(),
+                       false,
+                       lbfgs_data)
+end
+
+import Base.push!
+
+@doc "Push a new {s,y} pair into a L-BFGS operator." ->
+function push!(op :: LBFGSOperator, s :: Vector, y :: Vector)
+
+  ys = dot(y, s);
+  if ys <= 1.0e-20
+    warn("Rejecting L-BFGS {s,y} pair")
+    return
+  end
+  data = op.data;
+  insert = data.insert;
+
+  data.s[:, insert] = s;
+  data.y[:, insert] = y;
+  data.ys[insert] = ys;
+
+  # Update arrays a and b used in forward products.
+  if !op.inverse
+    data.b[:, insert] = y / sqrt(ys);
+
+    for i = 1 : data.mem
+      k = mod(insert + i - 1, data.mem) + 1;
+      if data.ys[k] != 0
+        data.a[:, k] = data.s[:, k];   # B₀ = I.
+
+        for j = 1 : i - 1
+          l = mod(insert + j - 1, data.mem) + 1;
+          if data.ys[l] != 0
+            data.a[:, k] += dot(data.b[:, l], data.s[:, k]) * data.b[:, l];
+            data.a[:, k] -= dot(data.a[:, l], data.s[:, k]) * data.a[:, l];
+          end
+        end
+        data.a[:, k] /= sqrt(dot(data.s[:, k], data.a[:, k]));
+      end
+    end
+  end
+
+  op.data.insert = mod(insert, data.mem) + 1;
+  return
+end
+
 end  # module