From 5ca099f864b500fb092797f4ecd14974b524c67b Mon Sep 17 00:00:00 2001
From: Jonathan Stickel <jjstickel@gmail.com>
Date: Fri, 14 Jan 2022 12:37:48 -0700
Subject: [PATCH] allow initial-guess methods to solve for lambda

---
 src/solvers.jl    | 26 ++++++++++++++++++++------
 src/validators.jl |  9 ++++++---
 2 files changed, 26 insertions(+), 9 deletions(-)

diff --git a/src/solvers.jl b/src/solvers.jl
index f2437f8..ac02c8a 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -184,20 +184,34 @@ function solve(
     Ψ::RegularizationProblem,
     b::AbstractVector;
     alg = :gcv_svd,
+    method = Brent(),
     λ₁ = 0.0001,
     λ₂ = 1000.0,
+    λ₀ = 1.0,
+    kwargs...
 )
     b̄ = @>> b to_standard_form(Ψ)
     L1, L2, κ = Lcurve_functions(Ψ, b̄)
 
-    solution = @match alg begin
-        :gcv_tr  => @_ optimize(gcv_tr(Ψ, b̄, _), λ₁, λ₂) 
-        :gcv_svd => @_ optimize(gcv_svd(Ψ, b̄, _), λ₁, λ₂) 
-        :L_curve => @_ optimize(1.0 - κ(_), λ₁, λ₂) 
+    optfunc(x) = @match alg begin
+        :gcv_tr => gcv_tr(Ψ, b̄, x)
+        :gcv_svd => gcv_svd(Ψ, b̄, x)
+        :L_curve => 1.0 - κ(x)
         _ => throw("Unknown algorithm, use :gcv_tr, :gcv_svd, or :L_curve")
     end
-
-    λ = @> solution Optim.minimizer
+    if method in [Brent(), GoldenSection()] # univariate methods with bounds
+        solution = optimize(optfunc, log10(λ₁), log10(λ₂), method; kwargs...)
+        log10λ = @> solution Optim.minimizer
+    else # methods that use an intitial guess; gradient-free and gradient methods via
+        # default finite-difference work Ok, but `autodiff=:forward` is not working
+        solution = optimize(x->optfunc(first(x)), [log10(λ₀)], method,
+                            Optim.Options(; kwargs...))
+        # alternative method call, but cautioned against in the docs
+        #solution = optimize(x->optfunc(first(x)), [log10(λ₀)]; method=method, kwargs...)
+        log10λ = @> solution Optim.minimizer first
+    end
+    λ = 10^log10λ
+    
     x̄ = solve(Ψ, b̄, λ)
     x = @>> x̄ to_general_form(Ψ, b) 
      
diff --git a/src/validators.jl b/src/validators.jl
index 1710faf..7b20085 100644
--- a/src/validators.jl
+++ b/src/validators.jl
@@ -19,7 +19,8 @@ Vλ = gcv_tr(Ψ, b̄, 0.1)                         # V(λ) single λ value
 Vλ = @_ map(gcv_tr(Ψ, b̄, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ
 ```
 """
-function gcv_tr(Ψ::RegularizationProblem, b̄::AbstractVector, λ::AbstractFloat)
+function gcv_tr(Ψ::RegularizationProblem, b̄::AbstractVector, log10λ::AbstractFloat)
+    λ = 10^log10λ
     n = size(Ψ.Ā, 1)
     Aλ = Ψ.Ā * inv(Ψ.ĀĀ + λ^2.0 * Ψ.Iₚ) * Ψ.Ā'
     Iₙ = Matrix{Float64}(I, length(b̄), length(b̄)) 
@@ -82,7 +83,8 @@ Vλ = gcv_svd(Ψ, b̄, x̄₀, 0.1)                     # V(λ) single λ value
 Vλ = @_ map(gcv_svd(Ψ, b̄, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ
 ```
 """
-function gcv_svd(Ψ::RegularizationProblem, b̄::AbstractVector, λ::AbstractFloat)
+function gcv_svd(Ψ::RegularizationProblem, b̄::AbstractVector, log10λ::AbstractFloat)
+    λ = 10^log10λ
     F = Ψ.F̄
     n, p = size(F.U)
     z = F.U' * b̄
@@ -148,7 +150,8 @@ function curvature_functions(x̄::Function, L1::Function, L2::Function)
     ρ⁰(λ::AbstractFloat) = (log.(L1.(λ) .^ 2.0))[1]
     ηᵖ(λ::AbstractFloat) = (derivative(η⁰, λ))[1]
 
-    function κ(λ::AbstractFloat)
+    function κ(log10λ::AbstractFloat)
+        λ = 10^log10λ
         nᵖ = ηᵖ(λ)
         η = η⁰(λ)
         ρ = ρ⁰(λ)