some changes

src/globalization/trs_solvers/solvers/Dogleg.jl

@@ -25,8 +25,6 @@ function (dogleg::Dogleg)(∇f, H, Δ, p, scheme, mstyle; abstol = 1e-10, maxite
     n = length(∇f)
 
     # find the Cauchy point; assumes ∇f is not ≈ 0
-    @show isa(scheme.approx, Direct)
-    @show H
     d_cauchy = if isa(scheme.approx, Direct)
         -∇f * norm(∇f)^2 / (∇f' * H * ∇f)
     else

src/optimize/linesearch/linesearch.jl

@@ -6,7 +6,7 @@ struct LineSearch{S,LS,K}
     linesearcher::LS
     scaling::K
 end
-LineSearch() = LineSearch(DBFGS(), Backtracking(), InitialScaling(ShannoPhua()))
+LineSearch() = LineSearch(DBFGS(Direct()), Backtracking(), InitialScaling(ShannoPhua()))
 LineSearch(m) = LineSearch(m, Backtracking(), InitialScaling(ShannoPhua()))
 LineSearch(m, ls) = LineSearch(m, ls, InitialScaling(ShannoPhua()))
 

src/quasinewton/approximations/DBFGS.jl

@@ -11,74 +11,44 @@ summary(dbfgs::DBFGS{Inverse}) = "Inverse Damped BFGS"
 summary(dbfgs::DBFGS{Direct}) = "Direct Damped BFGS"
 
 function update!(scheme::DBFGS{<:Direct,<:Any,<:Any}, B, s, y)
-    # We could write this as
-    #     B .+= (y*y')/dot(s, y) - (B*s)*(s'*B)/(s'*B*s)
-    #     B .+= (y*y')/dot(s, y) - b*b'/dot(s, b)
-    # where b = B*s
-    # But instead, we split up the calculations. First calculate the denominator
-    # in the first term
+    # Procedure 18.2 (Damped BFGS Updating) Nocedal Wright 2nd edition
     σ = dot(s, y)
     ρ = inv(σ) # scalar
     # Then calculate the vector b
     b = B * s # vector temporary
     sb = dot(s, b)
     if σ ≥ scheme.theta * sb
-        θ = 1.0
-        # Calculate one vector divided by dot(s, b)
-        ρbb = inv(sb) * b
-        # And calculate
-        B .+= (inv(σ) * y) * y' .- ρbb * b'
+        θ = 1
+        r = y
+        σ̂ = σ
     else
         θ = 0.8 * sb / (sb - σ)
         r = y * θ + (1 - θ) * b
-        # Calculate one vector divided by dot(s, b)
-        ρbb = inv(dot(s, b)) * b
-        # And calculate
-        B .+= (inv(dot(s, r)) * r) * r' .- ρbb * b'
+        σ̂ = dot(s, r)
     end
+    B .= B .+ r*r'/σ̂ - b*b'/sb
 end
 function update(scheme::DBFGS{<:Direct,<:Any}, B, s, y)
     # As above, but out of place
-
     σ = dot(s, y)
     b = B * s
     sb = dot(s, b)
     if σ ≥ scheme.theta * sb
-        θ = 1.0
-        # Calculate one vector divided by dot(s, b)
-        ρbb = inv(sb) * b
-        # And calculate
-        return B = B .+ (inv(σ) * y) * y' .- ρbb * b'
+        θ = 1
+        r = y
+        σ̂ = σ
     else
         θ = 0.8 * sb / (sb - σ)
         r = y * θ + (1 - θ) * b
-        # Calculate one vector divided by dot(s, b)
-        ρbb = inv(dot(s, b)) * b
-        # And calculate
-        return B = B .+ (inv(dot(s, r)) * r) * r' .- ρbb * b'
+        σ̂ = dot(s, r)
     end
+    B = B + r*r'/σ̂ - b*b'/sb
 end
 function update(scheme::DBFGS{<:Inverse,<:Any}, H, s, y)
-    σ = dot(s, y)
-    ρ = inv(σ)
-    #   if isfinite(ρ)
-    C = (I - ρ * s * y')
-    H = C * H * C' + ρ * s * s'
-    #   end
-    H
+
 end
 function update!(scheme::DBFGS{<:Inverse,<:Any}, H, s, y)
-    σ = dot(s, y)
-    ρ = inv(σ)
-
-    if isfinite(ρ)
-        Hy = H * y
-        H .= H .+ ((σ + y' * Hy) .* ρ^2) * (s * s')
-        Hys = Hy * s'
-        Hys .= Hys .+ Hys'
-        H .= H .- Hys .* ρ
-    end
-    H
+
 end
 
 function update!(scheme::DBFGS{<:Inverse,<:Any}, A::UniformScaling, s, y)

src/quasinewton/approximations/SR1.jl

@@ -13,18 +13,18 @@ function update(scheme::SR1{<:Inverse}, H, s, y)
     w = s - H * y
     θ = dot(w, y) # angle between residual and change in gradient
     ρy = norm(y)
-    if abs(θ) ≥ T(1e-12) * norm(w) * ρy && !iszero(ρy)
-        H = H + (w * w') / θ
+    if abs(θ) ≥ T(1e-12) * norm(w) * ρy && !iszero(ρy) # why iszero here?
+        H = H + (w * w') / θ # N&W2E
     end
     H
 end
 function update(scheme::SR1{<:Direct}, B, s, y)
     T = real(eltype(s))
     res = y - B * s # resesidual in secant equation
     θ = dot(res, s) # angle between residual and change in state
-    if abs(inv(θ)) ≥ T(1e-12) * norm(res) * norm(s)
+    if abs(θ) ≥ T(1e-12) * norm(res) * norm(s)
         if true #abs(θ) ≥ 1e-12
-            B = B + (res * res') / θ
+            B = B + (res * res') / θ # checked against N&W2E
         end
     end
     B

test/optimize/interface.jl

@@ -863,7 +863,7 @@ end
         GradientDescent(),
         BFGS(Inverse()),
         BFGS(Direct()),
-        DBFGS(),
+        DBFGS(Direct()),
         SR1(Inverse()),
         SR1(Direct()),
         DFP(),