Distances

ModelDistance{D <: AbstractDistance} <: AbstractSignalDistance

A model distance operates on signals and works by fitting an LTI model to the signals before calculating the distance. The distance between the LTI models is defined by the field distance. This is essentially a wrapper around the inner distance that handles the fitting of a model to the signals. How the model is fit is determined by fitmethod.

Arguments:

• fitmethod::FitMethod: LS, TLS or PLR
• distance::D: The inner distance between the models

Example:

using SpectralDistances
innerdistance = OptimalTransportRootDistance(domain=Continuous(), β=0.005, p=2)
dist = ModelDistance(TLS(na=30), innerdistance)

The top level distance type

All subtypes of this type operates on rational transfer functions

All subtypes of this type operates on signals

BuresDistance <: AbstractDistance

Distance between pos.def. matrices

CoefficientDistance{D, ID} <: AbstractCoefficientDistance

Distance metric based on model coefficients

Arguments:

• domain::D: Discrete or Continuous
• distance::ID = SqEuclidean(): Inner distance between coeffs

DiscreteGridTransportDistance{DT} <: AbstractDistance

Optimal transport between two measures on a common discrete grid.

DiscreteGridTransportDistance(Cityblock(), Float32, n, n)

The constructor takes an inner distance, the numeric type and the length of the two measures.

DiscretizedRationalDistance{WT, DT} <: AbstractRationalDistance

This distance discretizes the spectrum before performing the calculations.

Arguments:

• w::WT = LinRange(0.01, 0.5, 300): Frequency set
• distmat::DT = distmat_euclidean(w, w): DESCRIPTION

EnergyDistance <: AbstractSignalDistance

std(x1) - std(x2) This distance can be added to a loss function to ensure that the energy in the two signals is the same. Some of the optimal transport-based distances are invariant to the energy in the signal, requiring this extra cost if that invariance is not desired. Combining distances is done by putting two or more in a tuple. Usage: combined_loss = (primary_distance, EnergyDistance())

EuclideanRootDistance{D, A, F1, F2} <: AbstractRootDistance

Simple euclidean distance between roots of transfer functions

Arguments:

• domain::D: Discrete or Continuous
• assignment::A = SortAssignement(imag): Determines how roots are assigned. An alternative is HungarianAssignement
• transform::F1 = identity: DESCRIPTION
• weight : A function used to calculate weights for the induvidual root distances. A good option is residueweight
• p::Int = 2 : Order of the distance

HungarianRootDistance{D, ID <: Distances.PreMetric, F} <: AbstractRootDistance

Similar to EuclideanRootDistance but does the pole assignment using the Hungarian method.

Arguments:

• domain::D = Continuous(): Discrete or Continuous
• distance::ID = SqEuclidean(): Inner distance
• transform::F = identity: If provided, this Function transforms all roots before the distance is calculated

KernelWassersteinRootDistance{D, F, DI} <: AbstractRootDistance

A kernel version of the root distance

Arguments:

• domain::D = Continuous(): Discrete or Continuous
• λ::Float64 = 1.0: Kernel precision, lower value means wider kernel.
• transform::F = identity: If provided, this Function transforms all roots before the distance is calculated

OptimalTransportHistogramDistance{DT} <: AbstractDistance

Optimal transport between two histograms. If you pass two vectors to this distance two histograms will be computed automatically (fit(Histogram,x)). Pass two histograms to get better manual control.

Arguments:

• p::Int = 1: order

OptimalTransportRootDistance{D, F1, F2, S} <: AbstractRootDistance

The Sinkhorn distance between roots. The weights are provided by weight, which defaults to residueweight.

Arguments:

• domain::D = Continuous(): Discrete or Continuous
• transform::F1 = identity: Probably not needed.
• weight::F2 =simplex_residueweight: A function used to calculate weights for the induvidual root distances.
• β::Float64 = 0.01: Amount of entropy regularization
• p::Int = 2 : Order of the distance
• divergence::S = nothing: A divergence that penalizes creation and destruction of mass.

RationalCramerDistance{DT} <: AbstractRationalDistance

Similar to RationalOptimalTransportDistance but does not use inverse functions.

Arguments:

• domain::DT = Continuous(): Discrete or Continuous
• p::Int = 2: order
• interval = (-(float(π)), float(π)): Integration interval

CationalOptimalTransportDistance{DT, MT} <: AbstractRationalDistanCe

calculates the Wasserstein distance using the closed-form sold*c for (d,c) in zip(D,C)s and inverse cumulative functions.

C ArgumentC:

• domain::DT = Continuous(): Discrete or [Continuous](@red*c for (d,c) in zip(D,C)r
• magnitude::MT = Identity():
• interval = (-(float(π)), float(π)): Integration interval

SymmetricDistance{D} <: AbstractDistance

Evaluates d(x,y) - 0.5(d(x,x) + d(y,y))

WelchLPDistance{AT <: Tuple, KWT <: NamedTuple, F} <: AbstractWelchDistance

Lᵖ distance between welch spectra, mean(abs(x1-x2)^p).

#Arguments:

• args::AT = (): These are sent to welch_pgram
• kwargs::KWT = NamedTuple(): These are sent to welch_pgram
• p::Int = 2: Order of the distance
• normalized::Bool = true: Normlize spectrum to sum to 1 (recommended)
• transform::F = identity: Optional function to apply to the spectrum, example log1p or sqrt. Must not produce negative values, so log is not a good idea. The function is applied like this: transform.(x1).

WelchOptimalTransportDistance{DT, AT <: Tuple, KWT <: NamedTuple} <: AbstractWelchDistance

Calculates the Wasserstein distance between two signals by estimating a Welch periodogram of each.

Arguments:

• distmat::DT: you may provide a matrix array for this
• args::AT = (): Options to the Welch function
• kwargs::KWT = NamedTuple(): Options to the Welch function
• p::Int = 2 : Order of the distance

discrete_grid_transportcost(x::AbstractVector{T}, y::AbstractVector{T}, p = 2, tol=sqrt(eps(T)); inplace=false) where T

Calculate the optimal-transport cost between two vectors that are assumed to have the same support, with sorted and equdistant support points.

The calculated cost corresponds to the following

D = [abs2((i-j)/(n-1)) for i in 1:n, j in 1:n] # note the 1/(n-1)
Γ = discrete_grid_transportplan(x, y)
dot(Γ, D)

• ´p´ is the power of the distance
• ´inplaceif true,x` is overwritten.

discrete_grid_transportplan(x::AbstractVector{T}, y::AbstractVector{T}, tol=sqrt(eps(T)); inplace=false) where T

Calculate the optimal-transport plan between two vectors that are assumed to have the same support, with sorted support points.

• ´inplaceif true,x` is overwritten.
• It's possible to supply keyword argument g to provide preallocated memory for the transport plan. Make sure it's all zeros if you do.

distmat(dist, e::AbstractVector; normalize, kwargs...) -> Any

Compute the symmetric, pairwise distance matrix using the specified distance.

• normalize: set to true to normalize distances such that the diagonal is zero. This is useful for distances that are not true distances due to d(x,y) ≠ 0 such as the OptimalTransportRootDistance

This function uses multiple threads if available and copies the distance object to each thread in case it stores an internal cache.

distmat_euclidean(e1::AbstractVector, e2::AbstractVector, p = 2)

The euclidean distance matrix between two vectors of complex numbers

distmat_euclidean!(D, e1::AbstractVector, e2::AbstractVector, p = 2) = begin

In-place version

domain(d::AbstractDistance)

domain(d) -> Discrete{Int64}

Return the domain of the distance

domain_transform(d::AbstractDistance, e)

Change domain of roots

precompute(d::AbstractDistance, As, threads=true)

Perform computations that only need to be donce once when several pairwise distances are to be computed

Arguments:

• As: A vector of models
• threads: Us multithreading? (true)

ConvOptimalTransportDistance <: AbstractDistance

Distance between matrices caluclated using sinkhorn_convolutional. This type automatically creates a workspace object that is reused between invocations. Functions that internally performs threaded computations copy this object to each thread, but if manual threading is performed, this must be handled manually, e.g.:

using ThreadTools
dists = [deepcopy(dist) for _ in 1:Threads.nthreads()]
D = tmap(eachindex(X)) do i
    evaluate(dists[Threads.threadid()], Q, X[i]; kwargs...)
end

It's important to tune the two parameters below, see the docstring for sinkhorn_convolutional for more help.

• To get sharp barycenters, a smaller β around 0.001 is recommended.

• To get smooth distance profiles (distance_profile), a slightly higher β than for barycenters is recommended. β around 0.01-0.05 should do fine.

• β = 0.001

• dynamic_floor = -10.0

• invariant_axis::Int = 0 If this is set to 1 or 2, the distance will be approximately invariant to translations along the invariant axis. As an example, to be invariant to a spectrogram being shifted slightly in time, set invariant_axis = 2.

See also SlidingConvOptimalTransportDistance

SCWorkspace(A, B, β)

Workspace object for sinkhorn_convolutional. Manually construct this in order to save allocations between consequtive calls to the solver.

Arguments:

• A: The first matrix
• B: The second matrix
• β: the regularization parameter

SlidingConvOptimalTransportDistance

Similar to ConvOptimalTransportDistance but lets the shorter signal slide across the longer signal and returns the minimum distance. Equivalent to calling

minimum(distance_profile(d::ConvOptimalTransportDistance, a, b; kwargs...))

Γ, u, v = IPOT(C, a, b; β=1, iters=1000)

The Inexact Proximal point method for exact Optimal Transport problem (IPOT) (Sinkhorn-like) algorithm. C is the cost matrix and a,b are vectors that sum to one. Returns the optimal plan and the dual potentials. See also sinkhorn. β does not have to go to 0 for this alg to return the optimal distance, in fact, if β is set too low, this alg will encounter numerical problems.

A Fast Proximal Point Method for Computing Exact Wasserstein Distance Yujia Xie, Xiangfeng Wang, Ruijia Wang, Hongyuan Zha https://arxiv.org/abs/1802.04307

mask_filter(X::Matrix, th=0.25)
mask_filter!(Y, X=Y, th=0.25)

Mask out parts of a spectrogram that only contains impulsive noise. ´th ∈ (0,1) determines the amount of masking, a higher value masks out more. This function is to be applied after something like normalize_spectrogram.

normalize_spectrogram(S, dynamic_floor = default_dynamic_floor(S))

Γ, u, v = sinkhorn(C, a, b; β=1e-1, iters=1000)

The Sinkhorn algorithm. C is the cost matrix and a,b are vectors that sum to one. Returns the optimal plan and the dual potentials. This function is relatively slow, see also sinkhorn_log! IPOT and sinkhorn_log for faster algorithms.

sinkhorn_convolutional(A::AbstractMatrix, B::AbstractMatrix; β = 0.001, kwargs...)

sinkhorn_convolutional(w::SCWorkspace{T}, A::AbstractMatrix, B::AbstractMatrix; β = 0.01, τ = 1 / eps(T), iters = 1000, tol = 1.0e-6, ϵ = eps(T) ^ 2, verbose = false, initUV = true) where T

Calculate the entropically regularizaed Sinkhorn distance between two matrices where the ground cost is squared euclidean. This function uses an efficient convolutional algorithm and is much more efficient than the corresponding sinkhorn_log! in this special case.

This function spends most of its time doing matrix multiplications. You may want to tune BLAS.set_num_threads in order to maximize performance. 1-2 threads is often best, especially if you thread outside of calls to this function.

Arguments:

• w: workspace object
• A: The first matrix
• B: The second matrix
• β: regularization parameter. To get smooth distance profiles (distance_profile), a slightly higher β than for barycenters is recommended. β around 0.01 should do fine.
• τ: stabilization parameter
• iters: maximum number of iterations
• tol: tolerance (change in oen of the dual variables)
• ϵ: other stabilization parameter
• verbose: print stuff?
• initUV: if true, initializes dual variables in w to one. It can be useful to set this to false if you have a good initial guess (warm start). Set the corresponding w.U, w.V to provide the initial guess.

Same as sinkhorn_log but operates in-place to save memory allocations. This function has higher performance than sinkhorn_log, but might not work as well with AD libraries.

This function can be made completely allocation free with the interface sinkhorn_log(w::SinkhornLogWorkspace{T}, C, a, b; kwargs...)

The sinkhorn_log! solver also accepts a keyword argument check_interval = 20 that determines how often the convergence criteria is checked. If β is large, the algorithm might converge very fast and you can save some iterations by reducing the check interval. If β is small and the algorithm requires many iterations, a larger number saves you from computing the check too often.

The workspace w is created linke this: w = SinkhornLogWorkspace(FloatType, length(a), length(b))

Γ, u, v = sinkhorn_log(C, a, b; β=1e-1, iters=1000, tol=1e-8)

The Sinkhorn algorithm (log-stabilized). C is the cost matrix and a,b are vectors that sum to one. Returns the optimal plan and the dual potentials. See also sinkhorn_log! for a faster implementation operating in-place, and IPOT for a potentially more exact solution.

When this function is being differentiated, warnings about inaccurate solutions are turned off. You may choose to manually asses the error in the constrains by ea, eb = SpectralDistances.ot_error(Γ, a, b).

The IPOT algorithm: https://arxiv.org/pdf/1610.06519.pdf

Γ, u, v = sinkhorn_unbalanced(C, a, b, divergence; β=1e-1, iters=1000, tol=1e-8)

The Unbalanced Sinkhorn algorithm (log-stabilized). C is the cost matrix and a,b are vectors that are not required to sum to one.

Ref: "Sinkhorn Divergences for Unbalanced Optimal Transport" https://arxiv.org/abs/1910.12958 Makes use of UnbalancedOptimalTransport.jl

ot_jump(D, P1, P2)

Solve the optimal transport problem using JuMP. This function is only available if using JuMP, GLPK.

Arguments:

• D: Distance matrix
• P1: Weight vector 1
• P2: Weight vector 2

ot_convex(D, P1, P2)

Solve the optimal transport problem using Convex.jl. This function is only available if using Convex.jl, GLPK.

Arguments:

• D: Distance matrix
• P1: Weight vector 1
• P2: Weight vector 2

D̃, S = complete_distmat(D, W, λ = 2)

Takes an incomplete squared Euclidean distance matrix D and fills in the missing entries indicated by the mask W. W is a BitArray or array of {0,1} with 0 denoting a missing value. Returns the completed matrix and an SVD object that allows reconstruction of the generating point set X.

NOTE This function is only available after using Convex, SCS.

Arguments:

• D: The incomplete matrix
• W: The mask
• λ: Regularization parameter. A higher value enforces better data fitting, which might be required if the number of entries in D is very small.

Example:

using Distances
P = randn(2,40)
D = pairwise(SqEuclidean(), P)
W = rand(size(D)...) .> 0.3 # Create a random mask
W = (W + W') .> 0           # It makes sense for the mask to be symmetric
W[diagind(W)] .= true
D0 = W .* D                 # Remove missing entries

D2, S = complete_distmat(D0, W)

@show (norm(D-D2)/norm(D))
@show (norm(W .* (D-D2))/norm(D))

The set of points that created D can be reconstructed up to an arbitrary rotation and translation, X contains the reconstruction in the d first rows, where d is the dimension of the point coordinates. To reconstruct X using S, do

X  = Diagonal(sqrt.(S.S))*S.Vt

# Verify that reconstructed `X` is correct up to rotation and translation
A = [X[1:2,:]' ones(size(D,1))]
P2 = (A*(A \ P'))'
norm(P-P2)/norm(P) # Should be small

Ref: Algorithm 5 from "Euclidean Distance Matrices: Essential Theory, Algorithms and Applications" Ivan Dokmanic, Reza Parhizkar, Juri Ranieri and Martin Vetterli https://arxiv.org/pdf/1502.07541.pdf