Logistic regression

Binary classification via logistic regression in Julia.

This page comes from a single Julia file: logistic1.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: logistic1.ipynb, or open it in binder here: logistic1.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "ADTypes"
        "InteractiveUtils"
        "LaTeXStrings"
        "LinearAlgebra"
        "MIRTjim"
        "Optim"
        "Plots"
        "Random"
        "Statistics"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using ADTypes: AutoForwardDiff
using InteractiveUtils: versioninfo
using LaTeXStrings
using LinearAlgebra: dot, eigvals
using MIRTjim: prompt
using Optim: optimize
import Optim # Options
using Plots: default, gui, savefig
using Plots: histogram!, plot, plot!, scatter, scatter!
using Random: seed!
using Statistics: mean
default(); default(markersize=6, linewidth=2, markerstrokecolor=:auto, label="",
 tickfontsize=12, labelfontsize=18, legendfontsize=18, titlefontsize=18)

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? prompt(:prompt) : prompt(:draw);

Data

Generate synthetic data from two classes

if !@isdefined(yy)
    seed!(0)
    n0 = 60
    n1 = 50
    mu0 = [-1, 1]
    mu1 = [1, -1]
    v0 = mu0 .+ randn(2,n0) # class -1
    v1 = mu1 .+ randn(2,n1) # class 1
    nex = 0
    if true # 2017-01-18
        nex = 4 # extra dim (beyond the 2 shown) to make "larger scale"
        v0 = [v0; rand(nex,n0)] # (2+nex, n0)
        v1 = [v1; rand(nex,n1)] # (2+nex, n1)
    end
    M = n0 + n1 # how many samples
    yy = [-ones(Int, n0); ones(Int, n1)] # (M) labels
    vv = [[v0 v1]; ones(1,M)] # (npar, M) training data - with intercept
    npar = 3 + nex # unknown parameters
end;

Scatter plot and initial decision boundary

if !@isdefined(ps)
    x0 = [-1; 3; rand(nex); 5]
    v1p = range(-1,1,101) * 4
    v2p_fun(x) = @. (-x[end] - x[1] * v1p) / x[2]

    ps = plot(aspect_ratio = 1, size = (550, 500), legend = :topright,
     xaxis = (L"v_1", (-4, 4), [-4 -1 0 1 4]),
     yaxis = (L"v_2", (-4, 4), [-4 -1 0 1 4]),
    )
    plot!(v1p, v2p_fun(x0), color=:red, label="initial")
    plot!(v1p, v1p, color=:yellow, label="ideal")
    alpha = 0.7
    scatter!(v0[1,:], v0[2,:], color = :green; alpha)
    scatter!(v1[1,:], v1[2,:], color = :blue, marker = :square; alpha)
    # savefig(ps, "demo_fgm1_ogm1_s0.pdf")
end
ps

prompt()

Cost function

Logistic regression with Tikhonov regularization involves minimizing the following cost function:

\[f(x) = 1_M' h.(A x) + (β/2) ‖ x ‖_2^2\]

where $h(z) = \log(1 + e^{-z})$ is the logistic loss function.

Here $A$ is $M × N$ matrix with $M$ samples of $N$ features along each row (typically including the intercept $1$). The $m$th row of $A$ has already been multiplied by the $m$th binary class label that is ±1.

The cost function gradient is $∇ f(x) = A' \dot{h}.(A x) + β x$, and its Lipschitz constant is $‖A‖_2^2 / 4 + β$.

After optimizing $x$, the classifier is simply $\text{sign}(⟨v,x⟩)$ where the feature vector $v$ typically includes the intercept $1$.

if !@isdefined(cost)
    pot(t) = log(1 + exp(-t)) # logistic
    dpot(t) = -1 / (exp(t) + 1)
    tmp = vv * vv' # (npar, npar) covariance
    tmp = eigvals(tmp)
    @show maximum(tmp) / minimum(tmp)
    pLip = maximum(tmp) / 4 # 1/4 comes from logistic curvature

    reg = 2^0 # todo: use cross validation to select
    Lip = pLip + reg # Lipschitz constant

    A = yy .* vv' # M × N matrix of features times labels
    gfun(x) = A' * dpot.(A * x) + reg * x # gradient
    if false
        tmp = gfun(x0)
        @show size(tmp)
    end

    cost(x::AbstractVector) = sum(pot, A * x) + reg/2 * sum(abs2, x)
    cost(x::AbstractMatrix) = cost.(eachcol(x)) ## to handle arrays
end;

maximum(tmp) / minimum(tmp) = 79.5564516910791

GD

Iterate GD

tol = 1e-6
tol_gd = tol
tol_n1 = tol
tol_o1 = tol

function gd(x0, gfun::Function, niter::Int)
    xg = copy(x0)
    xgs = copy(xg)
    for n in 1:niter
        xold = xg
        xg -= 1/Lip * gfun(xg)
        if false && (norm(xg - xold, Inf) / norm(x0, Inf) < tol_gd)
            @show n
            break
        end
        any(isnan, xg) && throw("nan")
        # projected GD ??
        # xg = xg / norm(xg) # decision boundary unaffected by norm!
        xgs = [xgs xg] # archive
    end
    return xgs
end

if !@isdefined(xgs)
    niter_gd = 300
    xgs = gd(x0, gfun, niter_gd)
    pgs = plot(xgs', xlabel="Iteration", title = "GD")
end
pgs

prompt()

Nesterov FGM

do_restart = true;

function fgm(x0, grad, Lip::Real, niter::Int)
    re_nest = Int[]
    told = 1
    x = copy(x0)
    zold = copy(x0)
    xns = copy(x)
    for n in 1:niter
        xold = copy(x)
        grad = gfun(x)
        znew = x - 1/Lip * grad
        zdiff = znew - zold # dk - for restart

        tnew = 1/2 * (1 + sqrt(1 + 4 * told^2))

        x = znew + (told - 1) / tnew * (znew - zold)
        zold = copy(znew)
        told = tnew
        if false && (norm(x - xold, Inf) / norm(x0, Inf) < tol_n1)
            @show n
            break
        end
        any(isnan, x) && throw("nan")

        if do_restart && (dot(grad, zdiff) > 0) # dk new version
            told = 1
            x = copy(znew) # check
            zold = copy(x)
            @show "nest. restart", n
            push!(re_nest, n)
        end
        xns = [xns x] # archive
    end
    return xns, re_nest
end

if !@isdefined(xns)
    niter_n1 = 300
    niter_n1 = 500 # nex
    xns, re_nest = fgm(x0, gfun, Lip, niter_n1)
    pns = plot(xns', xlabel = "Iteration", title = "FGM")
end
pns

prompt()

OGM

Optimized gradient method

function ogm(
    x0,
    gfun,
    Lip::Real,
    niter::Int,
)
    re_ogm1 = Int[]
    told = 1
    x = copy(x0)
    zold = copy(x0)
    x1s = copy(x)
    for n in 1:niter
        xold = x
        grad = gfun(x)
        znew = x - 1/Lip * grad

        tnew = 1/2 * (1 + sqrt(1 + 4 * told^2))

        x = znew + (told - 1) / tnew * (znew - zold) + told/tnew * (znew - xold)
        zdiff = znew - zold # dk - for restart
        zold = znew
        told = tnew
        if false && (norm(x - xold, Inf) / norm(x0, Inf) < tol_n1)
            @show n
            break
        end
        any(isnan, x) && throw("nan")

        if do_restart && (dot(grad, zdiff) > 0) # dk new version
            push!(re_ogm1, n)
            told = 1
            x = znew # check
            zold = x # dk fixed from x0
            @info "ogm1 restart $n"
        end
        x1s = [x1s x] # archive
    end
    return x1s, re_ogm1
end


if !@isdefined(x1s)
    niter_o1 = 300
    niter_o1 = 500 # nex
    x1s, re_ogm1 = ogm(x0, gfun, Lip, niter_o1)
    po1 = plot(x1s', xlabel = "Iteration", title = "OGM")
end
po1

prompt()

L-BFGS Quasi-Newton optimizer

opt = Optim.Options(
 store_trace = true,
 show_warnings = false,
 extended_trace = true, # for trace of x
)
outq = optimize(cost, gfun, x0, opt;
    inplace = false, autodiff = AutoForwardDiff())
xqs = hcat(Optim.x_trace(outq)...)
xq = outq.minimizer
xh = xqs[:,end] # final estimate

7-element Vector{Float64}:
  1.5208531001394034
 -1.8351989096919836
 -0.21781865185010876
 -0.5977703450489538
 -0.5768183189101287
  0.33008172864648383
 -0.017115054005187102

Impartial version of x ͚

xh_tmp = [xgs[:,end] xns[:,end] x1s[:,end] xqs[:,end]]
xh = vec(mean(xh_tmp[:,2:end], dims=2)); # GD too slow to include

Plot cost

ifun(xs) = 0:(size(xs,2)-1)
extra = do_restart ? " (restart)" : ""
pc = plot(xaxis = ("iteration", (0,10)), yaxis = ("Cost function",))
plot!(0:niter_gd, cost(xgs) .- cost(xh), label = "GD" * extra)
plot!(0:niter_n1, cost(xns) .- cost(xh), label = "FGM" * extra)
plot!(0:niter_o1, cost(x1s) .- cost(xh), label = "OGM1" * extra)
plot!(ifun(xqs), cost(xqs) .- cost(xh), label = "QN")

prompt()

Plot decision boundaries

if true
    psh = deepcopy(ps)
    v2p = @. (-xh[end] - xh[1] * v1p) / xh[2]
    plot!(psh, v1p, v2p, color = :magenta, label = "final")
# savefig(psh, "demo-fgm1-fgm1a.pdf")
end
psh

prompt()

Plot iterate convergence

efun1(x) = vec(sqrt.(sum(abs2, x .- xh, dims=1)))
efun(x) = do_restart ? log10.(efun1(x)) : efun1(x)

pic = plot(
 xaxis = ("Iteration", (0, 40+10*nex), 0:20:80),
 yaxis = do_restart ?
  (L"\log_{10}(‖ \mathbf{x}_k - \mathbf{x}_* ‖)", (-3, 1), -3:1) :
  (L"‖ \mathbf{x}_k - \mathbf{x}_* ‖", (0, 8), [-1, 0, 8]),
 legend = :topright,
)
plot!(ifun(xgs), efun(xgs), color = :green, label = "GD")
plot!(ifun(xns), efun(xns), color = :blue, label = "Nesterov FGM" * extra)
plot!(ifun(x1s), efun(x1s), color = :red, label = "OGM1" * extra)
plot!(ifun(xqs), efun(xqs), label = "QN", marker = :o)
if do_restart
    scatter!(re_nest, efun(xns[:, re_nest .+ 1]), color = :blue)
    scatter!(re_ogm1, efun(x1s[:, re_ogm1 .+ 1]), color = :red)
end
pic

prompt()

# savefig demo_fgm1_ogm1c # restart
# savefig demo_fgm1_ogm1b # no restart

Plot 1D separation

inprod0 = [v0; ones(1,n0)]' * xh
inprod1 = [v1; ones(1,n1)]' * xh

accuracy0 = round(count(<(0), inprod0) / n0 * 100, digits=1)
accuracy1 = round(count(>(0), inprod1) / n1 * 100, digits=1)

plot(xaxis=("⟨x,v⟩",))
bins = -12:12
alpha = 0.5
histogram!(inprod0; alpha, bins, color = :green, linecolor = :green,
 label = "class 0: $accuracy0%")
histogram!(inprod1; alpha, bins, color = :blue, linecolor = :blue,
 label = "class 1: $accuracy1%")

prompt()

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

12-element Vector{SubString{String}}:
 "Julia Version 1.12.5"
 "Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)"
 "Build Info:"
 "  Official https://julialang.org release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)"
 "  GC: Built with stock GC"
 "Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()

Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [47edcb42] ADTypes v1.21.0
  [6e4b80f9] BenchmarkTools v1.6.3
  [aaaa29a8] Clustering v0.15.8
  [35d6a980] ColorSchemes v3.31.0
  [3da002f7] ColorTypes v0.12.1
  [c3611d14] ColorVectorSpace v0.11.0
  [717857b8] DSP v0.8.4
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.17.0
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.10.0
  [587475ba] Flux v0.16.9
  [a09fc81d] ImageCore v0.10.5
  [9ee76f2b] ImageGeoms v0.11.2
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.1.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.21.0
  [7035ae7a] MIRT v0.18.3
  [170b2178] MIRTjim v0.26.0
  [eb30cadb] MLDatasets v0.7.21
  [efe261a4] NFFT v0.14.3
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.10
  [429524aa] Optim v2.0.1
  [91a5bcdd] Plots v1.41.6
  [f27b6e38] Polynomials v4.1.1
  [2913bbd2] StatsBase v0.34.10
  [1986cc42] Unitful v1.28.0
  [d6d074c3] VideoIO v1.6.1
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0

This page was generated using Literate.jl.