diff --git a/Project.toml b/Project.toml
index 3262bae..4db2e05 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "LineSearches"
 uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.5.0"
+version = "7.5.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/LineSearches.jl b/src/LineSearches.jl
index 5fed4a7..4020dc9 100644
--- a/src/LineSearches.jl
+++ b/src/LineSearches.jl
@@ -29,10 +29,10 @@ function make_ϕdϕ(df, x_new, x, s)
         x_new .= x .+ α.*s
 
         # Evaluate ∇f(x+α*s)
-        NLSolversBase.value_gradient!(df, x_new)
+        f_x_new, g_x_new = NLSolversBase.value_gradient!(df, x_new)
 
         # Calculate ϕ(a_i), ϕ'(a_i)
-        NLSolversBase.value(df), real(dot(NLSolversBase.gradient(df), s))
+        return f_x_new, real(dot(g_x_new, s))
     end
     ϕdϕ
 end
@@ -42,48 +42,18 @@ function make_ϕ_dϕ(df, x_new, x, s)
         x_new .= x .+ α.*s
 
         # Evaluate ∇f(x+α*s)
-        NLSolversBase.gradient!(df, x_new)
+        g_x_new = NLSolversBase.gradient!(df, x_new)
 
         # Calculate ϕ'(a_i)
-        real(dot(NLSolversBase.gradient(df), s))
+        return real(dot(g_x_new, s))
     end
     make_ϕ(df, x_new, x, s), dϕ
 end
 function make_ϕ_dϕ_ϕdϕ(df, x_new, x, s)
-    function dϕ(α)
-        # Move a distance of alpha in the direction of s
-        x_new .= x .+ α.*s
-
-        # Evaluate f(x+α*s) and ∇f(x+α*s)
-        NLSolversBase.gradient!(df, x_new)
-
-        # Calculate ϕ'(a_i)
-        real(dot(NLSolversBase.gradient(df), s))
-    end
-    function ϕdϕ(α)
-        # Move a distance of alpha in the direction of s
-        x_new .= x .+ α.*s
-
-        # Evaluate ∇f(x+α*s)
-        NLSolversBase.value_gradient!(df, x_new)
-
-        # Calculate ϕ'(a_i)
-        NLSolversBase.value(df), real(dot(NLSolversBase.gradient(df), s))
-    end
-    make_ϕ(df, x_new, x, s), dϕ, ϕdϕ
+    make_ϕ_dϕ(df, x_new, x, s)..., make_ϕdϕ(df, x_new, x, s)
 end
 function make_ϕ_ϕdϕ(df, x_new, x, s)
-    function ϕdϕ(α)
-        # Move a distance of alpha in the direction of s
-        x_new .= x .+ α.*s
-
-        # Evaluate ∇f(x+α*s)
-        NLSolversBase.value_gradient!(df, x_new)
-
-        # Calculate ϕ'(a_i)
-        NLSolversBase.value(df), real(dot(NLSolversBase.gradient(df), s))
-    end
-    make_ϕ(df, x_new, x, s), ϕdϕ
+    make_ϕ(df, x_new, x, s), make_ϕdϕ(df, x_new, x, s)
 end
 
 include("types.jl")
diff --git a/src/initialguess.jl b/src/initialguess.jl
index 18d1dcf..282f2d6 100644
--- a/src/initialguess.jl
+++ b/src/initialguess.jl
@@ -74,7 +74,7 @@ function (is::InitialQuadratic{T})(ls, state, phi_0, dphi_0, df) where T
         # If we're at the first iteration
         αguess = is.α0
     else
-        αguess = 2 * (NLSolversBase.value(df) - state.f_x_previous) / dphi_0
+        αguess = 2 * (phi_0 - state.f_x_previous) / dphi_0
         αguess = NaNMath.max(is.αmin, state.alpha*is.ρ, αguess)
         αguess = NaNMath.min(is.αmax, αguess)
         # if αguess ≈ 1, then make it 1 (Newton-type behaviour)
@@ -180,8 +180,10 @@ function (is::InitialHagerZhang)(ls::Tls, state, phi_0, dphi_0, df) where Tls
         # and the user has not provided an initial step size (is.α0 is NaN),
         # then we
         # pick the initial step size according to HZ #I0
-        state.alpha = _hzI0(state.x, NLSolversBase.gradient(df),
-                            NLSolversBase.value(df),
+        # phi_0 is (or should be) equal to NLSolversBase.value(df, state.x) 
+        # TODO: Make the gradient available as part of the state?
+        g_x = NLSolversBase.gradient!(df, state.x)
+        state.alpha = _hzI0(state.x, g_x, phi_0,
                             is.αmax,
                             convert(eltype(state.x), is.ψ0)) # Hack to deal with type instability between is{T} and state.x
         if Tls <: HagerZhang
diff --git a/test/initial.jl b/test/initial.jl
index 15de7c3..4fc62f0 100644
--- a/test/initial.jl
+++ b/test/initial.jl
@@ -90,14 +90,14 @@
     state.f_x_previous = 2*phi_0
     is = InitialQuadratic(snap2one=(0.9,Inf))
     is(ls, state, phi_0, dphi_0, df)
-    @test state.alpha == 1.0
+    @test state.alpha == 0.5
     @test ls.mayterminate[] == false
 
     # Test Quadratic snap2one
     ls = HagerZhang()
     state = getstate()
     state.f_x_previous = 2*phi_0
-    is = InitialQuadratic(snap2one=(0.75,Inf))
+    is = InitialQuadratic(snap2one=(0.5,Inf))
     is(ls, state, phi_0, dphi_0, df)
     @test state.alpha == 1.0
     @test ls.mayterminate[] == false