Interpolations and Barycenters

Some distances distance define the existence of a shortest path, a geodesic. An interpolation is essentially a datapoint on that shortest path. We provide some functionality to interpolate between different spectra and models under transport-based metrics.

Below is an example usage of interpolations. We initially create two random systems, we then define the distance under which to interpolate and then calculate the frequency response for some different values of the interpolation parameter $t \in (0,1)$

using SpectralDistances, ControlSystemsBase, Distances, Plots, Random
plotly()
Random.seed!(0)

n = 4
r1 = complex.(-0.01 .+ 0.001randn(3), 2randn(3))
r1 = ContinuousRoots([r1; conj.(r1)])

r2 = complex.(-0.01 .+ 0.001randn(3), 2randn(3))
r2 = ContinuousRoots([r2; conj.(r2)])

r1,r2 = normalize_energy.((r1, r2))

A1 = AR(r1)
A2 = AR(r2)

##
fig1   = plot()
t      = 0.1
dist   = RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(0., exp10(1.01)))
interp = SpectralDistances.interpolator(dist, A1, A2)
w      = exp10.(LinRange(-1.5, 1, 300))
for t = LinRange(0, 1, 7)
    Φ = clamp.(interp(w,t), 1e-10, 100)
    plot!(w, sqrt.(Φ), xscale=:log10, yscale=:log10, line_z = t, lab="", xlabel="", title="W_2", ylims=(1e-3, 1e1), colorbar=false, l=(1,), c=:viridis)
end

rdist  = EuclideanRootDistance(domain = Continuous(), p = 2)
interp = SpectralDistances.interpolator(rdist, A1, A2, normalize=false)
fig2   = plot()
for t = LinRange(0, 1, 7)
    Φ = interp(w,t)
    plot!(w, sqrt.(Φ), xscale=:log10, yscale=:log10, line_z = t, lab="", xlabel="", title="RD", ylims=(1e-3, 1e1), colorbar=false, l=(1,), c=:viridis)
end

fig3 = plot()
Φ1   = bode(tf(A1), w)[1][:]
Φ2   = bode(tf(A2), w)[1][:]
for t = LinRange(0, 1, 7)
    plot!(w, (1-t).*Φ1 .+ t.*Φ2, xscale=:log10, yscale=:log10, line_z = t, lab="", xlabel="Frequency", title="L_2", ylims=(1e-3, 1e1), colorbar=false, l=(1,), c=:viridis)
end

fig = plot(fig1, fig2, fig3, layout=(3,1))

Barycenters

A barycenter is a generalization the the arithmetic mean to metrics other than the Euclidean. A barycenter between models is calculated like this

bc = barycenter(distance, models)

It can be useful to provide some options to the solvers:

options = (solver=sinkhorn_log!, tol=1e-8, iters=1_000_000, γ=0.0, uniform=true, inneriters=500_000, innertol=1e-6)
distance = OptimalTransportRootDistance(domain=Continuous(), p=2, β=0.01, weight=simplex_residueweight)
bc = barycenter(distance, models; options...)

We can plot the barycenters:

using SpectralDistances, ControlSystemsBase, Plots
models   = examplemodels(3)
distance = OptimalTransportRootDistance(domain=Continuous())
bc       = barycenter(distance, models)
w        = exp10.(LinRange(-0.5, 0.5, 350)) # Frequency vector
G        = tf.(models) # Convert models to transfer functions from ControlSystemsBase.jl
plot()
bodeplot!.(G, Ref(w), plotphase=false, lab="Input models", linestyle=:auto)
bodeplot!(tf(bc), w, plotphase=false, lab="Barycenter", xscale=:identity, c=:green)

Barycenters between spectrograms

We can also calculate a barycenter between spectrograms (or arbitrary matrices) using an efficient convolutional method. The most important parameter to tune in order to get a good result, apart from the regularization parameter β, is the dynamic_floor. This parameter determines where (in log space) the floor of the PSD is. This serves as a denoising, why the barycenter appears with a very dark background in the image below.

using SpectralDistances, DSP, Plots
N     = 24_000
t     = 1:N
f     = range(0.8, stop=1.2, length=N)
y1    = sin.(t .* f) .+ 0.1 .* randn.()
y2    = sin.(t .* reverse(f .+ 0.5)) .+ 0.1 .* randn.()
S1,S2 = spectrogram.((y1,y2), 1024)

A = [S1,S2]
β = 0.0001     # Regularization parameter (higher implies more smoothing and a faster, more stable solution)
λ = [0.5, 0.5] # Barycentric coordinates (must sum to 1)
B = barycenter_convolutional(A, β=β, tol=1e-6, iters=2000, ϵ=1e-100, dynamic_floor=-2)
plot(
    plot(S1, title="S1"),
    plot(B, title="Barycenter"),
    plot(S2, title="S2"),
    layout=(1,3),
    colorbar=false
)

Note that in order to calculate the barycenter, the sum of each input spectrogram is normalized.

This function works for any vector of matrices as long as all entries are positive and each matrix has an equal sum.

For a more thourogh example, see whistle.jl.

Trade off between frequency and time

There is currently no way of having different costs between transport in time and transport along the frequency axis other than to change the resolution of the spectrogram.

Barycentric coordiantes

The inverse problem to that of finding a barycenter is that of finding the barycentric coordinates λ of a query point $Q$, such that the resulting barycenter is as close as possible to the query point. Given a set of rational spectra $\left\{ G_i \right\}$, a nonlinear projection of a spectrum $Q$ onto this set can be obtained by solving the following nested optimization problem

\[\begin{aligned} λ &= \argmin_{\bar{λ}} \, W\big(Q, Q^*(\bar{λ})\big)\\ Q^*(\bar{λ}) &= \argmin_{\bar{Q}} \sum_i \bar{λ}_i W(G_i, \bar{Q}) \end{aligned}\]

where $λ$ are the barycentric coordinates belonging to the probability simplex. Problems of this type are sometimes referred to as histogram regression.

A nonlinear projection onto a basis consisting of spectra can be useful for, e.g., spectral dictionary learning, basis pursuit, topic modelling, denoising and detection. The function barycentric_coordinates is available for select distances:

SpectralDistances.barycentric_coordinates — Function

λ = barycentric_coordinates(pl, ql, p, q; options, kwargs...)

Compute the barycentric coordinates λ such that sum(λᵢ W(pᵢ,q) for i in eachindex(p)) is minimized.

This function works best with the sinkhorn_log! solver, a large β (around 1) and small tolerance. These are set using kwargs....

Arguments:

pl: Atoms in measures p, vector, length n_measures, of matrices of size n_dims × n_atoms
ql: Atoms in measure q
p: Measures p, a matrix of weight vectors, size n_atoms × n_measures that sums to 1
q: the veight vector for measure q, length is n_atoms
options: For the Optim solver. Defaults are options = Optim.Options(store_trace=false, show_trace=false, show_every=0, iterations=20, allow_f_increases=true, time_limit=100, x_tol=1e-5, f_tol=1e-6, g_tol=1e-6, f_calls_limit=0, g_calls_limit=0)
solver: = sinkhorn_log! solver
tol: = 1e-7 tolerance
β: = 0.1 entropy regularization. This function works best with rather large regularization, hence the large default value.
kwargs: these are sent to the solver algorithm.

Example:

using SpectralDistances, ControlSystemsBase, Optim
models = examplemodels(10)

d = OptimalTransportRootDistance(
    domain = SpectralDistances.Continuous(),
    p      = 2,
    weight = residueweight,
    β      = 0.01,
)
Xe = barycenter(d, models, solver=sinkhorn_log!)

G = tf.(models)
plot()
pzmap!.(G)
pzmap!(tf(Xe), m=:c, title="Barycenter OptimalTransportRootDistance", lab="BC")

options = Optim.Options(store_trace       = true,
                        show_trace        = false,
                        show_every        = 1,
                        iterations        = 50,
                        allow_f_increases = true,
                        time_limit        = 100,
                        x_tol             = 1e-7,
                        f_tol             = 1e-7,
                        g_tol             = 1e-7,
                        f_calls_limit     = 0,
                        g_calls_limit     = 0)


method = LBFGS()
λ = barycentric_coordinates(d, models, Xe, method,
    options = options,
    solver  = sinkhorn_log!,
    robust  = true,
    uniform = true,
    tol     = 1e-6,
)
bar(λ, title="Barycentric coorinates")

G = tf.(models)
plot()
pzmap!.(G, lab="")
pzmap!(tf(Xe), m = :c, title = "Barycenter OptimalTransportRootDistance", lab = "BC")
# It's okay if the last system dot does not match the barycenter exactly, there are limited models to choose from.
pzmap!(G[argmax(λ)], m = :c, lab = "Largest bc coord", legend = true)

Interpolations and Barycenters

Barycenters

Barycenters between spectrograms

Trade off between frequency and time

Barycentric coordiantes

K-Barycenters