From 48bdfd42e59422b8507d54214ae15d7519810fdf Mon Sep 17 00:00:00 2001
From: Cooper <simpson.r.cooper@gmail.com>
Date: Tue, 6 Aug 2024 22:38:35 -0700
Subject: [PATCH] Changing some regularization parameters in GCK.

---
 src/linesearch.jl |  2 +-
 src/minimize.jl   |  2 +-
 src/solvers.jl    | 22 +++++++++++-----------
 3 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/linesearch.jl b/src/linesearch.jl
index 0abe25c..ffacd01 100644
--- a/src/linesearch.jl
+++ b/src/linesearch.jl
@@ -66,7 +66,7 @@ function search!(opt::SFNOptimizer, stats::Stats, x::S, f::F, fval::T, g::S, g_n
 
     #Update regularization
     # println("Accepted η: ", η)
-    opt.M = max(min(1e8, opt.M/η^2), 1e-15)
+    opt.M = max(min(1e8, opt.M/η^2), 1e-8)
     # println("Updated M: ", opt.M)
 
     return success
diff --git a/src/minimize.jl b/src/minimize.jl
index 899dafb..ffccc8b 100644
--- a/src/minimize.jl
+++ b/src/minimize.jl
@@ -156,7 +156,7 @@ function iterate!(opt::O, x::S, f::F1, fg!::F2, Hv::H, itmax::I, time_limit::T)
         if opt.linesearch && !search!(opt, stats, x, f, fval, grads, g_norm, Hv)
             break
         else
-            println("P-norm: ", norm(opt.solver.p))
+            # println("P-norm: ", norm(opt.solver.p))
             x .+= opt.solver.p
         end
         ##########
diff --git a/src/solvers.jl b/src/solvers.jl
index d75d239..5f09397 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -69,10 +69,11 @@ function step!(solver::GLKSolver, stats::Stats, Hv::H, g::S, g_norm::T, M::T, ti
     # println("Max E: ", maximum(E.values), " Min E: ", minimum(E.values))
 
     #Quadrature scaling factor
-    # c = eigmax(Hv, tol=1e-1) #+ sqrt(λ)
-    c = λ≤1 ? 1. : sqrt(λ)
+    c = eigmax(Hv, tol=1e-6)
+    # c = sqrt(eigmax(Hv, tol=1e-6)) #+ sqrt(λ)
+    # c = λ≤1 ? 1. : sqrt(λ)
     
-    println("c: ", c)
+    # println("c: ", c)
 
     #Shifts
     shifts = c^2*solver.quad_nodes .+ λ
@@ -121,7 +122,7 @@ function hvp_power(solver::GCKSolver)
     return 2
 end
 
-function GCKSolver(dim::I, type::Type{<:AbstractVector{T}}=Vector{Float64}, quad_order::I=25, krylov_order::I=0) where {I<:Integer, T<:AbstractFloat}
+function GCKSolver(dim::I, type::Type{<:AbstractVector{T}}=Vector{Float64}, quad_order::I=31, krylov_order::I=0) where {I<:Integer, T<:AbstractFloat}
 
     #Quadrature
     nodes, weights = gausschebyshevt(quad_order)
@@ -154,14 +155,12 @@ function step!(solver::GCKSolver, stats::Stats, Hv::H, g::S, g_norm::T, M::T, ti
     - Median
     =#
     # β = 1.0
-    # β = λ > 1 ? λ : one(λ) 
-    # β = eigmax(Hv, tol=1e-6)
-    # E = eigen(Matrix(Hv))
-    # println(E.values)
-    # println("λ: ", λ)
+    # β = λ > 1 ? λ : one(λ)
+    β = eigmax(Hv, tol=1e-6)
+    # β = sqrt(eigmax(Hv, tol=1e-6))
     # β = sum(E.values)/length(g) + λ
+    # β = λ≤1 ? 1. : sqrt(λ)
     # println("β: ", β)
-    β = 100
 
     #Shifts and scalings
     shifts = (λ-β) .* solver.quad_nodes .+ (λ+β)
@@ -182,8 +181,9 @@ function step!(solver::GCKSolver, stats::Stats, Hv::H, g::S, g_norm::T, M::T, ti
     push!(stats.krylov_iterations, solver.krylov_solver.stats.niter)
 
     #Update search direction
-    for i in eachindex(scales)
+    for i in eachindex(shifts)
         @inbounds solver.p .+= solver.quad_weights[i]*solver.krylov_solver.x[i]
+        # @inbounds solver.p .-= solver.quad_weights[i]*(scales[i]*Matrix(Hv)+shifts[i]*I)\g
     end
 
     solver.p .*= sqrt(β)