Skip to content

Commit 58f26a1

Browse files
MaxenceGollierMaxenceGollier
committed
start LM update
1 parent 79c772f commit 58f26a1

3 files changed

Lines changed: 249 additions & 2 deletions

File tree

‎src/LMModel.jl‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
export LMModel
2+
3+
@doc raw"""
4+
LMModel(J, F, v, σ, x0)
5+
6+
Given the unconstrained optimization problem:
7+
```math
8+
\min \tfrac{1}{2} \| F(x) \|^2,
9+
```
10+
this model represents the smooth LM subproblem:
11+
```math
12+
\min_s \ \tfrac{1}{2} \| F(x) + J(x)s \|^2 + \tfrac{1}{2} σ \|s\|^2
13+
```
14+
where `J` is the Jacobian of `F` at `x0` in sparse format or as a linear operator.
15+
`σ > 0` is a regularization parameter and `v` is a vector of the same size as `F(x0)` used for intermediary computations.
16+
"""
17+
mutable struct LMModel{T <: Real, V <: AbstractVector{T}, G <: Union{AbstractMatrix{T}, AbstractLinearOperator{T}}} <:
18+
AbstractNLPModel{T, V}
19+
J::G
20+
F::V
21+
v::V
22+
σ::T
23+
meta::NLPModelMeta{T, V}
24+
counters::Counters
25+
end
26+
27+
function LMModel(J::G, F::V, σ::T, x0::V) where {T, V, G}
28+
@assert length(x0) == size(J, 2)
29+
@assert length(F) == size(J, 1)
30+
meta = NLPModelMeta(
31+
length(x0),
32+
x0 = x0, # Perhaps we should add lvar and uvar as well here.
33+
)
34+
v = similar(F)
35+
return LMModel(J::G, F::V, v::V, σ::T, meta, Counters())
36+
end
37+
38+
function NLPModels.obj(nlp::LMModel, x::AbstractVector{T}) where{T}
39+
@lencheck nlp.meta.nvar x
40+
increment!(nlp, :neval_obj)
41+
nlp.v .= nlp.F
42+
mul!(nlp.v, nlp.J, x, one(T), one(T))
43+
return ( dot(nlp.v, nlp.v) + nlp.σ * dot(x, x) ) / 2
44+
end
45+
46+
function NLPModels.grad!(nlp::LMModel, x::AbstractVector{T}, g::AbstractVector{T}) where{T}
47+
@lencheck nlp.meta.nvar x
48+
@lencheck nlp.meta.nvar g
49+
increment!(nlp, :neval_grad)
50+
nlp.v .= nlp.F
51+
@. nlp.g = nlp.σ .* x
52+
mul!(nlp.v, nlp.J, x, one(T), one(T))
53+
mul!(g, nlp.J', nlp.v, one(T), one(T))
54+
return g
55+
end

‎src/LM_alg.jl‎

Lines changed: 193 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,195 @@
1-
export LM
1+
export LM, LMSolver, solve!
2+
3+
import SolverCore.solve!
4+
5+
mutable struct TRDHSolver # FIXME
6+
end
7+
8+
mutable struct LMSolver{
9+
T <: Real,
10+
G <: ShiftedProximableFunction,
11+
V <: AbstractVector{T},
12+
M <: AbstractLinearOperator{T},
13+
ST <: AbstractOptimizationSolver,
14+
PB <: AbstractRegularizedNLPModel,
15+
} <: AbstractOptimizationSolver
16+
xk::V
17+
∇fk::V
18+
mν∇fk::V
19+
Fk::V
20+
Fkn::V
21+
Jk::M
22+
ψ::G
23+
xkn::V
24+
s::V
25+
has_bnds::Bool
26+
l_bound::V
27+
u_bound::V
28+
l_bound_m_x::V
29+
u_bound_m_x::V
30+
subsolver::ST
31+
subpb::PB
32+
substats::GenericExecutionStats{T, V, V, T}
33+
end
34+
35+
function LMSolver(
36+
reg_nls::AbstractRegularizedNLPModel{T, V};
37+
subsolver = R2Solver,
38+
) where{T, V}
39+
x0 = reg_nls.model.meta.x0
40+
l_bound = reg_nls.model.meta.lvar
41+
u_bound = reg_nls.model.meta.uvar
42+
43+
xk = similar(x0)
44+
∇fk = similar(x0)
45+
mν∇fk = similar(x0)
46+
Fk = similar(x0, reg_nls.model.nls_meta.nequ)
47+
Fkn = similar(Fk)
48+
Jk = jac_op_residual(reg_nls.model, xk)
49+
xkn = similar(x0)
50+
s = similar(x0)
51+
has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf)) || subsolver == TRDHSolver
52+
if has_bnds
53+
l_bound_m_x = similar(xk)
54+
u_bound_m_x = similar(xk)
55+
@. l_bound_m_x = l_bound - x0
56+
@. u_bound_m_x = u_bound - x0
57+
else
58+
l_bound_m_x = similar(xk, 0)
59+
u_bound_m_x = similar(xk, 0)
60+
end
61+
62+
ψ =
63+
has_bnds ? shifted(reg_nls.h, xk, l_bound_m_x, u_bound_m_x, reg_nls.selected) :
64+
shifted(reg_nls.h, xk)
65+
66+
sub_nlp = LMModel(Jk, Fk, T(1), x0)
67+
subpb = RegularizedNLPModel(sub_nlp, ψ)
68+
substats = RegularizedExecutionStats(subpb)
69+
subsolver = subsolver(subpb)
70+
71+
return LMSolver(
72+
xk,
73+
∇fk,
74+
mν∇fk,
75+
Fk,
76+
Fkn,
77+
Jk,
78+
ψ,
79+
xkn,
80+
s,
81+
has_bnds,
82+
l_bound,
83+
u_bound,
84+
l_bound_m_x,
85+
u_bound_m_x,
86+
subsolver,
87+
subpb,
88+
substats
89+
)
90+
end
91+
92+
function SolverCore.solve!(
93+
solver::LMSolver{T, G, V},
94+
reg_nls::AbstractRegularizedNLPModel{T, V},
95+
stats::GenericExecutionStats{T, V};
96+
callback = (args...) -> nothing,
97+
x::V = reg_nls.model.meta.x0,
98+
atol::T = √eps(T),
99+
rtol::T = √eps(T),
100+
verbose::Int = 0,
101+
max_iter::Int = 10000,
102+
max_time::Float64 = 30.0,
103+
max_eval::Int = -1,
104+
σk::T = eps(T)^(1 / 5),
105+
σmin::T = eps(T),
106+
η1::T = √√eps(T),
107+
η2::T = T(0.9),
108+
γ::T = T(3),
109+
θ::T = 1/(1 + eps(T)^(1 / 5)),
110+
) where {T, V, G}
111+
reset!(stats)
112+
113+
# Retrieve workspace
114+
selected = reg_nls.selected
115+
h = reg_nls.h
116+
nls = reg_nls.model
117+
118+
xk = solver.xk .= x
119+
120+
# Make sure ψ has the correct shift
121+
shift!(solver.ψ, xk)
122+
123+
Fk = solver.Fk
124+
Fkn = solver.Fkn
125+
Jk = solver.Jk
126+
∇fk = solver.∇fk
127+
JdFk = solver.JdFk
128+
Jt_Fk = solver.Jt_Fk
129+
ψ = solver.ψ
130+
xkn = solver.xkn
131+
s = solver.s
132+
133+
has_bnds = solver.has_bnds
134+
if has_bnds
135+
l_bound = solver.l_bound
136+
u_bound = solver.u_bound
137+
l_bound_m_x = solver.l_bound_m_x
138+
u_bound_m_x = solver.u_bound_m_x
139+
end
140+
141+
# initialize parameters
142+
improper = false
143+
hk = @views h(xk[selected])
144+
if hk == Inf
145+
verbose > 0 && @info "LM: finding initial guess where nonsmooth term is finite"
146+
prox!(xk, h, xk, one(eltype(x0)))
147+
hk = @views h(xk[selected])
148+
hk < Inf || error("prox computation must be erroneous")
149+
verbose > 0 && @debug "LM: found point where h has value" hk
150+
end
151+
improper = (hk == -Inf)
152+
improper == true && @warn "LM: Improper term detected"
153+
improper == true && return stats
154+
155+
if verbose > 0
156+
@info log_header(
157+
[:outer, :inner, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :normJ, :arrow],
158+
[Int, Int, T, T, T, T, T, T, T, T, Char],
159+
hdr_override = Dict{Symbol, String}(
160+
:fx => "f(x)",
161+
:hx => "h(x)",
162+
:xi => "√(ξ1/ν)",
163+
:normx => "‖x‖",
164+
:norms => "‖s‖",
165+
:normB => "‖J‖²",
166+
:arrow => "R2N",
167+
),
168+
colsep = 1,
169+
)
170+
end
171+
172+
local ξ1::T
173+
local ρk::T = zero(T)
174+
175+
residual!(nls, xk, Fk)
176+
Jk = jac_op_residual(nls, xk)
177+
mul!(∇fk, Jk', Fk)
178+
fk = dot(Fk, Fk) / 2
179+
180+
σmax, found_σ = opnorm(Jk)
181+
found_σ || error("operator norm computation failed")
182+
ν = θ / (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ
183+
sqrt_ξ1_νInv = one(T)
184+
185+
@. mν∇fk = -ν * ∇fk
186+
187+
φ1(d) = let Fk = Fk, Jk = Jk,
188+
d -> dot(Fk, Fk) / 2
189+
end
190+
191+
return
192+
end
2193

3194
"""
4195
LM(nls, h, options; kwargs...)
@@ -143,7 +334,7 @@ function LM(
143334
Resid_hist[k] = nls.counters.neval_residual
144335

145336
# model for first prox-gradient iteration
146-
φ1(d) = begin
337+
φ1(d) = begin # || Fk ||^2/2 + d*Jk'*Fk
147338
jtprod_residual!(nls, xk, Fk, Jt_Fk)
148339
dot(Fk, Fk) / 2 + dot(Jt_Fk, d)
149340
end

‎src/RegularizedOptimization.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ include("splitting.jl")
1919
include("TR_alg.jl")
2020
include("TRDH_alg.jl")
2121
include("R2_alg.jl")
22+
include("LMModel.jl")
2223
include("LM_alg.jl")
2324
include("LMTR_alg.jl")
2425
include("R2DH.jl")

0 commit comments

Comments
 (0)