Jukka Aho (@ahojukka5) | Jäähdytysmentelmä

Jäähdytysmenetelmä1 (engl. simulated annealing, SA) on satunnaisuuteen ja todennäköisyyksiin perustuva globaali optimointistrategia. Menetelmän perusajatus on ottaa satunnaisia askelia parametriavaruudessa ja hyväksyä uusi luku minimiksi todennäköisyysfunktion avulla. Todennäköisyysfunktio antaa arvon 1 mikäli uusi arvo on pienempi kuin entinen, mutta antaa myös todennäköisyyden valita huonompi arvo jollakin todennäköisyydellä, joka pienenee iteraatioiden funktiona. Tällä tavalla algoritmi pystyy pakenemaan lokaalista minimistä.

using Plots

Tarkastellaan funktiota \(f(x) = x^2 - 10\cos(5x) + 10\), jolla on lukuisia lokaaleja optimipisteitä sekä tasan yksi globaali optimipiste origossa.

f = x -> x^2 - 10*cos(5*x) + 10
p = plot(-5:0.01:5, f)
savefig("figs/function_1d.svg")
p

svg

Kaikki gradienttipohjaiset optimointimenetelmät epäonnistuvat tällaisen funktion optimoinnissa, mikäli alkuarvaus ei ole todella hyvä.

Lähdetään naputtelemaan SA koodi. Ensin tarvitaan transitiomalli, eli mitenkä tässä parametriavaruudessa olisi tarkoitus tepastella. Valitaan esimerkiksi normaalijakaumasta satunnainen luku.

using Distributions

function transition_model(x)
    return x + rand(Normal())
end

transition_model (generic function with 1 method)

Sen lisäksi tarvitaan todennäköisyysfunktio. Yleensä \(P(e, e', T)\) määritellään niin että \(P = 1\), jos \(e' < e\) ja muutoin \(\exp(-(e'-e)/T)\), joten menemme sillä. Tässä \(e\) on “energia” joka minimoituu.

function acceptance_probability(f, f_new, T)
    if f_new < f
        return 1.0
    end
    return exp(-(f_new-f)/T)
end

acceptance_probability (generic function with 1 method)

Lisäksi tarvitaan jäähtymismalli, joka pienentää lämpötilaa iteraatioiden funktiona. Otetaan tässä yksinkertaisin mahdollinen, eli lämpötila pienenee lineaarisesti.

function temperature(fraction)
    return 1 - fraction
end

temperature (generic function with 1 method)

Eli: transitiomallin ehdottama uusi \(x\) valitaan uudeksi minimiksi, mikäli \(f(x_{new}) < f(x)\), mutta myös siinä tapauksessa, että funktion arvo pisteessä \(x_{new}\) on isompi, jollakin todennäköisyydellä joka pienenee iteraatioiden funktiona.

Varsinainen pääalgoritmi tekee satunnaisia siirtymiä parametriavaruudessa ja hyväksyy tai hylkää uuden minimin perustuen todennäköisyysmalliin. Tallennetaan trace jotta voidaan tutkia konvergenssia:

function optimize_sa(f, x0::T, acceptance_probability,
        transition_model, temperature; k_max=10) where T

    x = x0
    accepted = Int[]
    rejected = Int[]
    trace = T[]
    for k=1:k_max
        x_new = transition_model(x)
        push!(trace, x_new)
        temp = temperature(k/k_max)
        if acceptance_probability(f(x), f(x_new), temp) > rand()
            x = x_new
            push!(accepted, k)
        else
            push!(rejected, k)
        end
    end
    return trace, accepted, rejected
end

optimize_sa (generic function with 1 method)

Ajetaan optimointialgoritmia esimerkiksi 300 iteraatiota esimerkkifunktiolle:

using Random
Random.seed!(0)
trace, accepted, rejected = optimize_sa(f, 5.0,
    acceptance_probability, transition_model,
    temperature; k_max=300)

iters = 1:length(trace)
ia = iters[accepted]
ir = iters[rejected]
xa = trace[accepted]
xr = trace[rejected]
xa[end]

0.012187171992082857

Vaikuttaisi että optimi on löytynyt. Algoritmille ominaista on että se ei anna tarkkaa ratkaisua vaan jotakin sen läheltä. Loput voidaan sitten hoitaa esimerkiksi Newtonin menetelmällä.

scatter(ia, xa,
    marker=:o, label="Accepted",
    title="Trace of x")
p = scatter!(ir, xr,
    marker=:x, label="Rejected")
savefig("figs/trace_x.svg")
p

svg

scatter(ia, f.(xa),
    marker=:o, label="Accepted",
    title="Function f")
p = scatter!(ir, f.(xr),
    marker=:x, label="Rejected")
savefig("figs/function_f.svg")
p

svg

Jäljestä näkee että optimi on löytynyt aika lyhyen iteroinnin jälkeen, alle 50 iteraatiossa. Optimoinnin puolessa välissä näkee kuinka algoritmi on ottanut pienen satunnaishypyn joka ei ole itse asiassa parantanut tulosta.

Case study: Challengerin katastrofi.

https://tieku.fi/maailmankaikkeus/avaruuslennot/challengerin-katastrofi-astronautit-tiesivat-kuolevansa 28. tammikuuta 1986 avaruussukkula Challenger muuttui taivaalla jättimäiseksi tulipalloksi.

Räjähdyksen syyksi selvisi yksinkertainen virhe – ja pian myös tiedettiin, että mukana olleet seitsemän astronauttia olivat vielä hengissä sekunteja ensimmäisten räjähdysten vavisuttaessa sukkulaa.

Vian syyksi paljastui O-renkaat, jotka eivät kestäneet Floridan hyyviä pakkaskelejä. Analysoidaan nyt tuosta saatua dataa logistisella mallilla.

data = """
6 0 66  50  1
6 1 70  50  2
6 0 69  50  3
6 0 68  50  4
6 0 67  50  5
6 0 72  50  6
6 0 73 100  7
6 0 70 100  8
6 1 57 200  9
6 1 63 200 10
6 1 70 200 11
6 0 78 200 12
6 0 67 200 13
6 2 53 200 14
6 0 67 200 15
6 0 75 200 16
6 0 70 200 17
6 0 81 200 18
6 0 76 200 19
6 0 79 200 20
6 2 75 200 21
6 0 76 200 22
6 1 58 200 23
"""

function read_data(data)
    lines = split(strip(data), "\n")
    N = length(lines)
    x = zeros(N)
    y = zeros(N)
    for i in 1:N
        vals = split(lines[i], " ")
        x[i] = parse(Int, vals[3])
        y[i] = parse(Int, vals[2])
    end
    return x, y
end

x, y = read_data(data)
for i in 1:length(y)
    if y[i] != 0.0
        y[i] = 1.0
    end
end

Todennäköisyys, että O-rengas on vaurioitunut, on

\[\begin{equation} p\left(x\right)=\frac{1}{1+\exp\left(-k\left(x-x\_{0}\right)\right)} \end{equation}\]

Logistinen kustannus on \[\begin{equation} L\left(\theta\right)=-\frac{1}{m}\sum\_{i=1}^{m}\left[\hat{y}^{\left(i\right)}\log\left(p\left(\theta\right)\right)+\left(1-\hat{y}^{\left(i\right)}\right)\log\left(1-p\left(\theta\right)\right)\right] \end{equation}\]

function loss(θ; x=x, y=y)
    k, x0 = θ
    @assert length(x) == length(y)
    J = 0.0
    N = length(x)
    for i in 1:N
        z = k*(x[i] - x0)
        a = 1.0 / (1.0 + exp(-z))
        L = -1.0 * (y[i] * log(a) + (1.0 - y[i]) * log(1.0 - a))
        J += L
    end
    return J / N
end

loss (generic function with 1 method)

Kevyt visualisointi jotta nähdään miltä funktio näyttää:

function do_plot()
    ks = -3.0:0.01:0.0
    ws = collect(50.0 : 0.1 : 100.0)

    px = []
    py = []
    for ki in ks
        for wi in ws
            l = loss((ki, wi))
            if (isnan(l) || isinf(l))
                push!(px, ki)
                push!(py, wi)
                break
            end
        end
    end

    p = contour(ks, ws, (x, y) -> loss((x, y)))
    plot!(px, py, color=:black, lw=1.5, label=nothing)
    xlims!(-3.0, 0.0)
    ylims!(50.0, 90.0)
    xlabel!("Parameter k")
    ylabel!("Parameter x0")
    title!("Logistic cost")
    return p
end

p = do_plot()
savefig("figs/logistic_loss.svg")
p

svg

Nyt sitten optimoidaan. Eli minimoidaan logistista kustannusfunktiota. Jäähdytysmenetelmällä.

function transition_model(θ)
    k, x0 = θ
    k_new = k + 1/π * rand(Normal())
    x0_new = x0 + π * rand(Normal())
    return (k=k_new, x0=x0_new)
end

transition_model (generic function with 1 method)

using Random
Random.seed!(0)
trace, accepted, rejected = optimize_sa(loss,
    (k=-3.0, x0=50.0),
    acceptance_probability, transition_model,
    temperature; k_max=1000)

trace_k = [o.k for o in trace]
trace_x0 = [o.x0 for o in trace]
iters = 1:length(trace)
ia = iters[accepted]
ir = iters[rejected]
ka = trace_k[accepted]
x0a = trace_x0[accepted]
kr = trace_k[rejected]
x0r = trace_x0[rejected]
trace[accepted][end]

(k = -0.23451914148538833, x0 = 65.48773395385525)

scatter(ia, x0a,
    marker=:o, label="Accepted",
    title="Trace of x0")
p = scatter!(ir, x0r,
    marker=:x, label="Rejected")
savefig("figs/trace_x0.svg")
p

svg

scatter(ia, ka,
    marker=:o, label="Accepted",
    title="Trace of k",
    legend=:topleft)
p = scatter!(ir, kr,
    marker=:x, label="Rejected")
savefig("figs/trace_k.svg")
p

svg

do_plot()
plot!(ka, x0a, color=:blue, markercolor=:blue,
    marker=:o, alpha=0.3, label="Accepted")
p = plot!(kr, x0r, color=:red, markercolor=:red,
    marker=:x, alpha=0.3, label="Rejected", legend=:topleft)
savefig("figs/trace_x0_and_k.svg")
p

svg

Lyhyen konvergenssin jälkeen jää algoritmi pyörimään minimin tuntumaan. Valitulla menetelmällä ratkaisuksi saadaan

opt = trace[accepted][end]

(k = -0.23451914148538833, x0 = 65.48773395385525)

opt2 = (k = -0.23, x0 = 64.75)
xc = collect(minimum(x) : 0.1 : maximum(x))
val(x) = 1.0 / (1.0 + exp(-opt.k*(x - opt.x0)))
yc = val.(xc)
plot(xc, yc, lw=1.5, color=:black, label="fitted model")
p = scatter!(x, y, color=:white, label="data")
title!("Challenger space shuttle O-ring failures")
xlabel!("Ambient temperature (Fahrenheit)")
ylabel!("Probability to damage")
ylims!(-0.1, 1.1)
savefig("figs/challenger_model.svg")
p

svg

p = do_plot()
ks = -3.0:0.01:0.0
ws = collect(50.0 : 0.1 : 100.0)
p = contour!(p, ks, ws, (x,y) -> loss((x, y)), levels=[0.45, 0.48, 0.55, 0.70, 1.0])
p = scatter!(p, [opt.k], [opt.x0], color=:red, label="minimum", legend=:topleft)
savefig("figs/challenger_logistic_cost_optimized.svg")
p

┌ Warning: Multiple series with different levels share a colorbar. Colorbar may not reflect all series correctly.
└ @ Plots /home/jukka/.julia/packages/Plots/7U0ob/src/backends/gr.jl:514

svg

┌ Warning: Multiple series with different levels share a colorbar. Colorbar may not reflect all series correctly.
└ @ Plots /home/jukka/.julia/packages/Plots/7U0ob/src/backends/gr.jl:514

Jäähdytysmentelmä

Case study: Challengerin katastrofi.

Tags

Recent posts

Jäähdytysmentelmä

Neuroverkkojen matematiikkaa

Konversioita Julian datatyyppien välillä

Elementtimenetelmän kantafunktioiden laskeminen

Toisen kertaluvun differentiaaliyhtälön ratkaisu elementtimenetelmällä

Archives