Laplacian eigenmaps

This example illustrates Laplacian eigenmaps using the Julia language. These eigenmaps provide nonlinear dimensionality reduction for data lying near manifolds (rather than near subspaces).

See Belkin & Niyogi, 2003.

This page comes from a single Julia file: eigmap.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: eigmap.ipynb, or open it in binder here: eigmap.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([
        "ImagePhantoms"
        "InteractiveUtils"
        "LaTeXStrings"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
    ])
end

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

using ImagePhantoms: rect, phantom
using InteractiveUtils: versioninfo
using LaTeXStrings
using LinearAlgebra: norm, Diagonal, eigen, svd, svdvals
using MIRTjim: jim, prompt
using Plots: default, gui, plot, savefig, scatter
using Plots.PlotMeasures: px
using Random: seed!
default(); default(label = "", markerstrokecolor = :auto,
 tickfontsize = 12, labelfontsize = 16,
)
seed!(0)

Random.TaskLocalRNG()

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

isinteractive() && prompt(:prompt);

Generate data

The example considered here uses synthetic data consisting of images of rectangles of various widths and rotations.

The only latent parameters here are one of the rectangle widths and the rotation angle. So all of the images lie in a manifold of dimension 2 in a $16^2 = 256$ dimensional ambient space.

function rect_phantom( ;
    T = Float32,
    nx = 40,
    ny = 40,
    sigma1::Real = 0f0,
    sigma2::Real = 0f0,
    level1::Function = () -> 0 + sigma1 * randn(T),
    level2::Function = () -> 1 + sigma2 * randn(T),
    angle::Function = () -> rand(T) * π/4f0,
    width_min::Real = min(nx,ny) / 4f0 + 3,
    width_max::Real = min(nx,ny) / 1.6f0, # √2
    width::Function = () -> [width_min/2,
        width_min + rand(T) * (width_max - width_min)],
    center::Function = () -> 0 * (rand(T, 2) .- 1//2),
    oversample::Int = 6,
)
    ob = rect(Tuple(center()), Tuple(width()), angle())
    l1 = level1()
    l2 = level2()
    x = (1:nx) .- (nx-1)/2
    y = (1:ny) .- (ny-1)/2
    return l1 .+ (l2 - l1) * phantom(x, y, [ob], oversample),
        ob.angle[1], ob.width[2]
end

function make_phantoms(nx, ny, nrep)
    data = zeros(Float32, nx, ny, nrep)
    angles = zeros(nrep)
    widths = zeros(nrep)
    for iz in 1:nrep
        data[:,:,iz], angles[iz], widths[iz] = rect_phantom(; nx, ny)
    end
    return data, angles, widths
end


if !@isdefined(data)
    nx,ny = 40,40
    nrep = 500
    @time data, angles, widths = make_phantoms(nx, ny, nrep)
    pj = jim(data[:,:,1:72]; title = "72 of $nrep images",
     xaxis = false, yaxis = false, colorbar = :none, # book
     size = (400, 480), nrow = 9,
    )
end
# savefig(pj, "eigmap-data.pdf")

These data do not lie in a 2-dimensional subspace. To see that, we plot the first 40 singular values. Even in the absence of noise here, there are many singular values that are are far from zero.

tmp = svdvals(reshape(data, nx*ny, nrep))
ps = scatter(tmp; title = "Data singular values", widen=true,
 xaxis = (L"k", (1, 40), [1, 40, nx*ny]),
 yaxis = (L"σ_k", (0, 216), [0, 32, 48, 72, 215]),
)
# savefig(ps, "eigmap-svd.pdf")

prompt()

Eigenmaps

Apply one of the Laplacian eigenmap methods. First compute the distances between all pairs of data points (images).

There is little if any humanly visible structure in the distance map.

distance = [norm(d1 - d2) for
    d1 in eachslice(data; dims=3),
    d2 in eachslice(data; dims=3)
]
pd = jim(distance; title = L"‖ X_j - X_i \; ‖", xlabel = L"i", ylabel = L"j")
# savefig(pd, "eigmap-dis.pdf")

Compute the weight matrix that describes an affinity between data points $i$ and $j$. There are many ways to do this; here we follow the approach given in Sanders et al. 2016.

α = Float64(sum(abs2, distance) / nrep^2) # per Sanders eqn. (4)
W = @. exp(-distance^2 / α)
pw = jim(W; title = L"W_{ij}", xlabel = L"i", ylabel = L"j")
# savefig(pw, "eigmap-w.pdf")

Graph Laplacian

Compute the graph Laplacian from the weight matrix.

d = vec(sum(W; dims=2))
D = Diagonal(d)
L = D - W # Laplacian
pl = jim(L; title = L"L_{ij}", xlabel = L"i", ylabel = L"j", color=:cividis)
# savefig(pl, "eigmap-l.pdf")

Compute the generalized eigendecomposition to find solutions of $L v = λ D v$.

Because the vector $\mathbf{1}$ is in the null space of $L$, the first (smallest) generalized eigenvalue is 0.

F = eigen(L, Matrix(D))

pe = scatter(F.values;
 title = "Generalized eigenvalues",
 xlabel=L"k", ylabel=L"λ_k",
 xticks = [1,nrep],
 yticks = [0, 0.7, 0.9, 1],
)
# savefig(pe, "eigmap-eigval.pdf")

prompt()

Here we are interested in the first couple of generalized eigenvectors, excluding the one corresponding to $λ=0$. Those vectors serve as features for the reduced dimension.

features = 1000 * F.vectors[:,2:3]
pf = scatter(features[:,1], features[:,2];
  title = "Eigenmap features",
  xlabel = "feature 1", ylabel = "feature 2")
# savefig(pf, "eigmap-f.pdf")

prompt()

To examine whether those features encode the relevant manifold here, we plot the features with color-coding corresponding to the rectangle width and angle.

The first feature is correlated with rectangle angle. The second feature is correlated with rectangle width.

pc = plot(
 scatter(features[:,1], features[:,2], marker_z = widths, color=:cividis,
  xaxis = ("feature 1", (-7,8), -6:6:6),
  yaxis = ("feature 2", (-6,5), -4:2:4),
# colorbar_title = "width", colorbar_fontsize = 15,
  annotate = (7, 4, "width"),
 ),
 scatter(features[:,1], features[:,2], marker_z = angles, color=:cividis,
  xaxis = ("feature 1", (-7,8), -6:6:6),
  yaxis = ("feature 2", (-6,5), -4:2:4),
  xlabel="feature 1", ylabel="feature 2",
# colorbar_title = "angle",
  annotate = (7, 4, "angle"),
 );
 layout = (2,1), size = (600, 400), right_margin = 10px,
)
# savefig(pc, "eigmap-c.pdf")

prompt()

As another way to illustrate this correlation, suppose we examine all images where feature 1 is in some interval like $(-1.5,-1)$. These images are all of rectangles of similar orientation, but various widths.

tmp = data[:,:, -1.5 .< features[:,1] .< -1]
pf1 = jim(tmp; title = "Feature 1 set")
# savefig(pf1, "eigmap-pf1.pdf")

Conversely, the images where feature 2 is in some interval like $(0.2,0.5)$ are all rectangles of similar widths, but various rotations.

tmp = data[:,:, 0.2 .< features[:,2] .< 0.5]
pf2 = jim(tmp; nrow=2, title = "Feature 2 set", size=(600,300))
# savefig(pf2, "eigmap-pf2.pdf")

Examine PCA for comparison

tmp = reshape(data, nx*ny, nrep)
U2 = svd(tmp).U[:,1:2]
pup = jim(reshape(U2, nx, ny, 2); title="First 2 PCA components")
# savefig(pup, "eigmap-pca-u.pdf")

Projection onto a 2-dimensional subspace works poorly:

lr2 = reshape(U2 * (U2' * tmp), nx, ny, :)
plr = jim(lr2[:,:,1:88], "Projection onto 2-dimensional subspace")
# savefig(plr, "eigmap-pca-lr2.pdf")

Nevertheless, the corresponding features are reasonably correlated with width and angle.

feat_pca = (U2' * tmp)' # leading coefficients

pp = plot(
 scatter(feat_pca[:,1], feat_pca[:,2], marker_z = widths, color=:cividis,
  xlabel="feature 1", ylabel="feature 2",
  colorbar_title="width",
 ),
 scatter(feat_pca[:,1], feat_pca[:,2], marker_z = angles, color=:cividis,
  xlabel="feature 1", ylabel="feature 2",
  colorbar_title="angle",
 ),
 plot_title = "First two PCA coefficients",
)
# savefig(pp, "eigmap-pca-f.pdf")

prompt()

Reproducibility

This page was generated with the following version of Julia:

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

11-element Vector{SubString{String}}:
 "Julia Version 1.11.1"
 "Commit 8f5b7ca12ad (2024-10-16 10:53 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-16.0.6 (ORCJIT, znver3)"
 "Threads: 1 default, 0 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`
  [6e4b80f9] BenchmarkTools v1.5.0
  [aaaa29a8] Clustering v0.15.7
  [35d6a980] ColorSchemes v3.27.1
⌅ [3da002f7] ColorTypes v0.11.5
⌃ [c3611d14] ColorVectorSpace v0.10.0
  [717857b8] DSP v0.7.10
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.7.0
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.8.0
  [587475ba] Flux v0.14.25
  [a09fc81d] ImageCore v0.10.4
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.0.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.20.1
  [7035ae7a] MIRT v0.18.2
  [170b2178] MIRTjim v0.25.0
  [eb30cadb] MLDatasets v0.7.18
  [efe261a4] NFFT v0.13.5
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.5
  [429524aa] Optim v1.10.0
  [91a5bcdd] Plots v1.40.9
  [f27b6e38] Polynomials v4.0.11
  [2913bbd2] StatsBase v0.34.3
  [d6d074c3] VideoIO v1.1.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [44cfe95a] Pkg v1.11.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.