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, ControlSystems, 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, ControlSystems, 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 ControlSystems.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_coordinatesFunction
λ = 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, ControlSystems, 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)
source
barycentric_coordinates(d::ConvOptimalTransportDistance, X::Vector{<:AbstractMatrix}, q::AbstractMatrix; method = LBFGS(), kwargs...)

Calculate the barycentric coordinates of a vector of matrices X using the convolutional method.

Arguments:

Optim options

The default options are

options = Optim.Options(
        store_trace       = true,
        show_trace        = false,
        show_every        = 1,
        iterations        = 10,
        allow_f_increases = false,
        time_limit        = 150,
        x_tol             = 1e-3,
        f_tol             = 1e-3,
        g_tol             = 1e-4,
    )
source

Below is an example that assumes that you have access to a vector with a signal called y and the corresponding sample rate fs (in this example, we use the y constructed in the previous example. In the example, we simply use one of the input spectra as query point. This way, we know what barycentric coordinates we expect

using SpectralDistances, DSP, Plots
models = map(DSP.arraysplit(y1, 8*512, 0)) do y
    spectrogram(y, 512, fs = 24_000, window=hanning)
end

matrices = s1.(normalize_spectrogram.(models))

β      = 0.01
dist   = ConvOptimalTransportDistance(β=β)
Q      = matrices[2] # Use the second spectrogram as query
λ, res = barycentric_coordinates(dist, matrices, Q)
plot(
    bar(λ, title="Barycentric coordinates", lab=""),
    plot(spectrogram(y1, window=hanning), title="Signal"),
    heatmap(Q, title="Query point")
)

K-Barycenters

Below, we show an example of how one can run the K-barycenter algorithm on a collection of sound signals. sounds is expected to be of type Vector{Vector{T}}. The example further assumes that there is a vector of labels::Vector{Int} that contain the true classes of the datapoints, which you do not have in an unsupervised setting.

using SpectralDistances, ControlSystems
fitmethod = TLS(na=12)
models = SpectralDistances.fitmodel.(fitmethod, sounds)
G = tf.(models) # Convert to transfer functions for visualization etc.

##
using Clustering
dist = OptimalTransportRootDistance(domain=Continuous(), β=0.01, weight=simplex_residueweight)
@time clusterresult = SpectralDistances.kbarycenters(
    dist,
    models,
    n_classes, # number of clusters
    seed       = :rand,
    solver     = sinkhorn_log!,
    tol        = 2e-6,
    innertol   = 2e-6,
    iters      = 100000,
    inneriters = 100000,
    verbose    = true,
    output     = :best,
    uniform    = true,
    kiters     = 10
)

bc,ass = clusterresult.barycenters, clusterresult.assignments

# Visualize results
using MLBase, Plots.PlotMeasures, AudioClustering
newass,perm = AudioClustering.associate_clusters(labels,ass)
classinds   = 1:n_classes
yt          = (classinds, [label_strings[findfirst(labels .== i)] for i in classinds])

@show mean(labels .== newass)
cm = confusmat(n_classes,labels,newass)
heatmap(cm./sum(cm,dims=2), xlabel="Cluster assignment",ylabel="Best matching class", color=:viridis)
anns = [(reverse(ci.I)..., text(val,12,:gray)) for (ci,val) in zip(CartesianIndices(cm)[:], vec(cm))]
annotate!(anns)
yticks!(yt)
xticks!(yt, xrotation=45)
current()

The figure should look like the last figure in the paper.

SpectralDistances.ISAFunction
ISA(X, w = nothing; iters = 100, printerval = typemax(Int))

Iterative swapping algorithm from "On the Computation of Wasserstein barycenters", Giovanni Puccetti et al.

Arguments:

  • X: vector of d×k matrices where d is dimension and k number of atoms
  • w: weights. See the files test_barycenter.jl for different uses.
  • iters: maximum number of iterations
  • printerval: print this often
source
SpectralDistances.barycenterFunction
barycenter(d::OptimalTransportRootDistance, models; normalize = true, kwargs...)
barycenter(d, models)
barycenter(d, models, λ; normalize, uniform, solver, kwargs...)

Approximately calculate the barycenter supported on the same number of atoms as the number of poles in the models.

The solver can be selected by providing a keword argument, example: solver=IPOT.

Uses the algorithms from "Fast Computation of Wasserstein Barycenters"

Example:

models = examplemodels(10)

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

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

Arguments:

  • models: vector of AR models
  • normalize: make sure weights sum to 1
  • kwargs: are sent to the solver
source
SpectralDistances.barycenterMethod
barycenter(X::Vector{<:AbstractArray}, λ)

Calculate the weighted barycenter for point clouds in X. Each X[i] has the shame n_dims × n_atoms λ is the weight vector that should sum to 1.

source
SpectralDistances.barycenterMethod
barycenter(d::EuclideanRootDistance, models::AbstractVector, [λ])

Example:

models = examplemodels(10)

Xe = barycenter(EuclideanRootDistance(domain=SpectralDistances.Continuous(),p=2), models)

G = tf.(models)
plot()
pzmap!.(G)
pzmap!(tf(Xe), m=:c, title="Barycenter EuclideanRootDistance")
current()
source
SpectralDistances.barycenter_convolutionalFunction
barycenter_convolutional(models::Vector{<:DSP.Periodograms.TFR}, λ = Fill(1 / length(models), length(models)); dynamic_floor = default_dynamic_floor(models), kwargs...)

Covenience function for the calculation of spectrograms. This function transforms the spectrograms to log-power and adjusts the floor to dynamic_floor, followed by a normalization to sum to 1.

This function will be called if barycenter is called with ConvOptimalTransportDistance as first argument.

Arguments:

  • dynamic_floor: Sets the floor of the spectrogram in log-domain, i.e., all values below this will be truncated. The default value is based on a quantile of the spectrogram powers. If your spectrograms are mostly low entropy, you can try to increase this number to get sharper results.
  • kwargs: Same as for the base method

Example:

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=200, ϵ=1e-100, dynamic_floor=-2)
plot(plot(S1, title="S1"), plot(B, title="Barycenter"), plot(S2, title="S2"), layout=(1,3), colorbar=false)
source
SpectralDistances.barycenter_convolutionalMethod
barycenter_convolutional(A, [λ]; β = 0.01, iters = 1000, tol = 1e-9, ϵ = 1e-90, verbose = false)

Convolutional barycenters.

´βis the regularization andλ(optional) are the weights (barycentric coordinates). To reuse allocated space between successive calls, use the "workspace" method from the example below.ϵ` is a truncation parameter for numerical stability.

a1      = zeros(10, 10)
a1[2,2] = 1
a2      = zeros(10, 10)
a2[6,6] = 1
A       = [a1,a2]
β       = 0.01
λ       = [0.5, 0.5] # Barycentric coordinates, must sum to 1
w       = BCWorkspace(A, β)
b       = barycenter_convolutional(w,A,λ)
plot(heatmap(a1), heatmap(a2), heatmap(b))

Ref: J. Solomon, F. de Goes, G. Peyré, M. Cuturi, A. Butscher, A. Nguyen, T. Du, L. Guibas. Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains. 2015 https://people.csail.mit.edu/jsolomon/assets/convolutional_w2.compressed.pdf

source
SpectralDistances.barycentric_coordinatesMethod
barycentric_coordinates(d::ConvOptimalTransportDistance, X::Vector{<:AbstractMatrix}, q::AbstractMatrix; method = LBFGS(), kwargs...)

Calculate the barycentric coordinates of a vector of matrices X using the convolutional method.

Arguments:

Optim options

The default options are

options = Optim.Options(
        store_trace       = true,
        show_trace        = false,
        show_every        = 1,
        iterations        = 10,
        allow_f_increases = false,
        time_limit        = 150,
        x_tol             = 1e-3,
        f_tol             = 1e-3,
        g_tol             = 1e-4,
    )
source
SpectralDistances.barycentric_coordinatesMethod
λ = 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, ControlSystems, 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)
source
SpectralDistances.embeddingMethod
embedding([::Type{Vector}], m, [full=true])

Returns a Vector/Matrix containing the roots of m. full indicates whether or not to use all poles or only one half-plane.

source
SpectralDistances.kbarycentersMethod
kbarycenters(X, p, k; seed = :rand, kiters = 10, verbose = false, output = :best, kwargs...)

Clustering using K-barycenters. If you want to cluster spectra, consider the method that accepts models instead.

Arguments:

  • X: Support of input measures, Vector{Matrix} where each matrix is ndims×natoms
  • p: Weights of input measures Vector{Vector} where each matrix is of length n_atoms and should sum to 1.
  • k: number of clusters
  • seed: :rand or :eq
  • kiters: number of iterations
  • verbose: print stuff?
  • output: output lowest cost clustering or :latest?
  • kwargs: are sent to the inner solvers.
source
SpectralDistances.kbarycentersMethod
kbarycenters(d::OptimalTransportRootDistance, models::Vector{<:AbstractModel}, k; normalize = true, kwargs...)

This function is only available if using Clustering.

Example:

clusterresult = kbarycenters(
    dist,
    models,
    k,         # number of clusters
    seed       = :rand,
    solver     = sinkhorn_log!,
    tol        = 2e-6,
    innertol   = 2e-6,
    iters      = 100000,
    inneriters = 100000,
    verbose    = true,
    output     = :best,
    uniform    = true,
    kiters     = 10
)

The docs contain a more detailed example

Arguments:

  • models: A vector of models
  • k: number of clusters
  • normalize: Whether or not to normalize the weight vectors (recommended)
  • kwargs: are sent to inner solvers, (solver,tol,iters)
source