Skip to content
Snippets Groups Projects
Commit 72924e1a authored by Johan Grönqvist's avatar Johan Grönqvist
Browse files

Test more validation metrics for gradients

parent 03c5ec9f
Branches
No related tags found
No related merge requests found
...@@ -75,6 +75,7 @@ function bisect(; Δ=1e-8, kwargs...) ...@@ -75,6 +75,7 @@ function bisect(; Δ=1e-8, kwargs...)
λu = λc λu = λc
λc *= 2 λc *= 2
end end
λℓ = λc
end end
while λu-λℓ > Δ while λu-λℓ > Δ
...@@ -88,9 +89,18 @@ function bisect(; Δ=1e-8, kwargs...) ...@@ -88,9 +89,18 @@ function bisect(; Δ=1e-8, kwargs...)
λℓ λℓ
end end
# %%
ϵ = 1e-4
# %% # %%
bisect(ζ=.1) bisect(ζ=.1)
# %%
bisect(ζ=.1+ϵ)
# %%
bisect(ζ=.2-ϵ)
# %% # %%
bisect(ζ=.2) bisect(ζ=.2)
...@@ -104,16 +114,16 @@ bisect(ζ=.2) ...@@ -104,16 +114,16 @@ bisect(ζ=.2)
valid, P = validate(λ=λ, ζ=.1, ret_P=true) valid, P = validate(λ=λ, ζ=.1, ret_P=true)
# %% # %%
validate(λ=λ, ζ=.1, P=P) @show validate(λ=λ, ζ=.1, P=P) validate(λ=λ, ζ=.1+ϵ, P=P)
# %% # %%
validate(λ=λ, ζ=.1, P=P) bisect(ζ=.1, P=P)
# %% # %%
validate(λ=λ, ζ=.2, P=P) bisect(ζ=.1+ϵ, P=P)
# %% # %%
bisect(ζ=.1, P=P) bisect(ζ=.2-ϵ, P=P)
# %% # %%
bisect(ζ=.2, P=P) bisect(ζ=.2, P=P)
...@@ -129,9 +139,110 @@ P ...@@ -129,9 +139,110 @@ P
# %% # %%
valid, P = validate(λ=λ, ζ=.2, ret_P=true) valid, P = validate(λ=λ, ζ=.2, ret_P=true)
# %%
@show validate(λ=λ, ζ=.2-ϵ, P=P) validate(λ=λ, ζ=.2, P=P)
# %%
bisect(ζ=.1, P=P)
# %%
bisect(ζ=.1+ϵ, P=P)
# %%
bisect(ζ=.2-ϵ, P=P)
# %%
bisect(ζ=.2, P=P)
# %% # %%
P P
# %% # %%
# %% # %%
function validate_all(;λ, ζs, ω=1.0, ϵ=1e-8, P = nothing, ret_P = false)
model = JuMP.Model(MosekTools.Optimizer)
set_silent(model)
M = [ 1+ϵ 0; 0 ϵ ]
A(ζ) = [ λ 1; -ω^2 λ-2*ζ*ω ]
if isnothing(P)
lhss = []
@variable model P_var[1:2, 1:2] PSD
for ζ in ζs
lhs = A(ζ)'*P_var+P_var*A(ζ)+M
@constraint model lhs<=0 PSDCone()
push!(lhss, lhs)
end
optimize!(model)
# @show termination_status(model) value.(lhs) eigvals(value.(lhs))
valid = maximum([ maximum(eigvals(value.(lhs))) for lhs in lhss ]) < 0
if ret_P
(valid, value.(P_var))
else
valid
end
else
lhss = []
for ζ in ζs
lhs = A(ζ)'*P+P*A(ζ)+M
push!(lhss, lhs)
end
valid = maximum([ maximum(eigvals(lhs)) for lhs in lhss]) < 0
valid
end
end
# %%
function bisect_all(; Δ=1e-8, kwargs...)
if validate_all(λ=0.0; kwargs...)
λℓ = 0.0
λc = 1.0
while validate_all(λ = λc; kwargs...)
λℓ = λc
λc *= 2
end
λu = λc
else
λu = 0.0
λc = -1.0
while !validate_all(λ=λc; kwargs...)
λu = λc
λc *= 2
end
λℓ = λc
end
while λu-λℓ > Δ
λc = (λu+λℓ)/2
if validate_all(λ=λc; kwargs...)
λℓ = λc
else
λu = λc
end
end
λℓ
end
# %%
λ = bisect(ζ = .1)
# %%
valid, P = validate(λ=λ, ζ=.1, ret_P=true)
# %%
validate_all(λ=λ, ζs = [.1, .1+ϵ])
# %%
validate_all(λ=λ, ζs = [.2-ϵ, .2])
# %%
bisect_all(ζs=[.1, .2])
# %%
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment