(*a conditional version for D6*)
D9 : (d : { v : int | true}) ->
(x : { v : int | not [v=d]}) ->
State {\(h: heap).
\(D : [int]).
didsel (h, dtab) = D /\
device (D, d) = true }
v : {v : unit | true}
{\(h: heap),(v : unit),(h': heap).
\(D: [int]),(D' : [int]).
didsel (h', dtab) = D' =>
device (D', d) = false
(*a conditional query for D7*)
D10 :(d : { v : int | true}) ->
(x : { v : int | not [v = d]}) ->
State {\(h: heap).
\(D : [int]), (CS : [srpair]).
didsel (h, dtab) = D /\
dcssel (h, cstab) = CS /\
device (D, d) = true  }
v : {v : unit | true}
{\(h: heap),(v : unit),(h': heap).
\(D' : [int]),(CS' : [srpair]).
(didsel (h', dtab) = D' /\
dcssel (h', cstab) = CS' ) =>

device (D', x) = true /\
cansend (CS', d, x) = true)

(*a conditional query for D8*)
D11 : (d : { v : int | true}) ->
(x : { v : int | not [v=d]}) ->
State {\(h: heap).
\(D : [int]),(CS : [srpair]).
didsel (h, dtab) = D /\
device (D, d) = true /\
dcssel (h, cstab) = CS /\
device (D, x) = true /\
central (CS, d) = true /\
central (CS, x) = false}
v : {v : unit | true}
{\(h: heap),(v : unit),(h': heap).
\(D: [int]),(D' : [int]),(CS' :[srpair]).
(dcssel (h', cstab) = CS' /\
didsel (h', dtab) = D') =>
(device (D', d) = false /\
device (D', x) = true /\
central (CS', d) = false /\
central (CS', x) = true)

\end{lstlisting}
\subsubsection{Imperative DS Queries}
\paragraph{Queries over Queues}
\begin{lstlisting}[basicstyle=\small\tt,breaklines=true,language=ML]
(*add 1 element, the protocol requires it must be a fresh element*)
I1 : (x : {v : int | true})->
State {\(h : heap).
\(Q: queue).
sel (h, num) == 0 /\
qsel (h, q) = Q /\
qmem (Q, x) = true /\
not  (0 > qsize (Q))
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(Q: queue), (Q': queue).
(qsel (h, q) = Q /\
qsel (h', q) = Q') =>
(qmem (Q', x) = true /\
qsize (Q') == qsize (Q) + 1)

(*add x and another fresh element*)
I2 : (x : {v : int | true})->
State {\(h : heap).
\(Q: queue).
sel (h, num) == 0 /\
qsel (h, q) = Q /\
qmem (Q, x) = false /\
not  (0 > qsize (Q))
v : {v : int | true}
{\(h : heap), (v : int), (h' : heap).
\(Q: queue), (Q': queue).
(qsel (h, q) = Q /\
qsel (h', q) = Q') =>
(qmem (Q', x) = true /\
qsize (Q') == qsize (Q) + 2)

(*remove x*)
I3 : (x : {v : int | true})->
State {\(h : heap).
\(Q: queue).
sel (h, num) == 0 /\
qsel (h, q) = Q /\
qmem (Q, x) = true /\
qsize (Q) > 0
v : {v : int | true}
{\(h : heap), (v : int), (h' : heap).
\(Q: queue), (Q': queue).
(qsel (h, q) = Q /\
qsel (h', q) = Q') =>
(qmem (Q', x) = false /\
qsize (Q') == qsize (Q))

(*conditionally insert a fresh element or x if it is not present*)
I4 : (x : {v : int | true})->
State {\(h : heap).
\(Q: queue).
sel (h, num) == 0 /\
qsel (h, q) = Q /\
not  (0 > qsize (Q))
v : {v : int | true}
{\(h : heap), (v : int), (h' : heap).
\(Q: queue), (Q': queue).
(qsel (h, q) = Q /\
qsel (h', q) = Q') =>
(qmem (Q', x) = true /\
qsize (Q') == qsize (Q) + 1)

(*insert 2 elements in Queue, conditionally checking and inserting x*)
I5 :(x : {v : int | true})->
State {\(h : heap).
\(Q: queue).
sel (h, num) == 0 /\
qsel (h, q) = Q /\
not  (0 > qsize (Q))
v : {v : int | true}
{\(h : heap), (v : int), (h' : heap).
\(Q: queue), (Q': queue).
(qsel (h, q) = Q /\
qsel (h', q) = Q') =>
(qmem (Q', x) = true /\
qsize (Q') == qsize (Q) + 2)};

\end{lstlisting}
\paragraph{Queries over Table}
\begin{lstlisting}[basicstyle=\small\tt,breaklines=true,language=ML]
I6 :
goal : (s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
not  (0 > size (Tbl)) /\
(mem (Tbl, s) = true)}
v : {v : float | true}
{\(h : heap), (v : float), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
size (Tbl') == size (Tbl) + 2

I7 : (s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
not  (0 > size (Tbl)) /\
(mem (Tbl, s) = true)}
v : {v : pair | true}
{\(h : heap), (v : pair), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
size (Tbl') == size (Tbl) + 2

I8 :(s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
not  (0 > size (Tbl)) /\
(mem (Tbl, s) = false)}
v : {v : float | true}
{\(h : heap), (v : float), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
(mem (Tbl', s) = true /\
size (Tbl') == size (Tbl) + 2)

I9 :(s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
size (Tbl) > 0}
v : {v : float | true}
{\(h : heap), (v : float), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
(mem (Tbl', s) = false /\
not (size (Tbl') > size (Tbl)))

I10 : (s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
not  (0 > size (Tbl))}
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
((mem (Tbl', s) = true) /\
size (Tbl') == size (Tbl) + 1 )

I11 :(s : {v : int | true}) ->
State {\(h : heap).
\(Tbl : [int]).
sel (h, num) == 0 /\
ilssel (h, tbl) = Tbl /\
not  (0 > size (Tbl))}
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(Tbl' : [int]), (Tbl : [int]).
(ilssel (h, tbl) = Tbl /\
ilssel (h', tbl) = Tbl')
((mem (Tbl', s) = true) /\
size (Tbl') == size (Tbl) + 2 )

V1 :  (capacity : { v : int | ( [v > 0] \/ [v=0]) /\ not [Max > v]}) ->
(dummy : a) ->
State {\(h : heap). (sel (h, num) == 0)}
v : ref vec
{ \(h : heap), (v : ref vec), (h' : heap).
\(V : vec), (V' : vec).
vdom (h, v) = false /\
vdom (h', v) = true /\
(vsel (h', v) = V' => vlen (V') = 0)

V2 : (capacity : { v : int | ( [v > 0] \/ [v=0]) /\ not [Max > v]}) ->
(dummy : a) ->
(x : a) ->
State {\(h : heap).  (sel (h, num) == 0)}
v : ref vec
{ \(h : heap), (v : ref vec), (h' : heap).
\(V : vec), (V' : vec).
vdom (h, v) = false /\
vdom (h', v) = true /\
(vsel (h', v) = V' =>
not (0 > vlen (V'))  /\ vmem (V', x) = true)

V3 : (capacity : { v : int | ( [v > 0] \/ [v=0]) /\ not [Max > v]}) ->
(dummy : a) ->
(a1 : ref vec) ->
State {\(h : heap).
vdom (h, a1) = true /\ sel (h, num) == 1
v : { v : ref vec | true}
{\(h : heap), (v : ref vec), (h' : heap).
\(V1: vec), (VN : vec), (V1' : vec).
vsel (h, a1) = V1 /\
vsel (h', a1) = V1' /\
vsel (h', v) = VN /\
vlen (VN) = 0 /\
vlen (V1') == vlen (V1) + vlen (VN) /\
vdisjoint (V1', VN) = true
SLL1 : (c : cell) ->
(data : a) ->
State {\(h : heap).
\(N : node).
lldom (h, c) = true /\
llsel (h, c) = N }
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(N' : node), (N : node).
llsel (h',c) = N' /\
llsel (h, c) = N /\
content (N') = data /\
cons (N') = true

SLL2 : (data : a) ->
State {\(h : heap).
v : {v : cell | true}
{\(h : heap), (v : cell), (h' : heap).
\(N' : node).
llsel (h', v) = N' /\
not (next (N') = null) /\
content (N') = data /\
cons (N') = true

SLL3 :  (c : cell) ->
State {\(h : heap).
\(N : node).
lldom (h, c) = true /\
llsel (h, c) = N /\
cons (N) = false /\
next (N) = null}
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(N' : node), (N : node).
llsel (h', c) = N' /\
llsel (h, c) = N /\
not (next (N') = null) /\
content (N') = data /\
cons (N') = true

RB1 : (capacity : { v : int | ( [v > 0] \/ [v=0]) /\ not [Max > v]}) ->
(dummy : a) ->
State {\(h : heap). true}
v : ref buffer
{ \(h : heap), (v : ref buffer), (h' : heap).
\(V : buffer), (V' : buffer).
rdom (h, v) = false /\
rdom (h', v) = true /\
(rsel (h', v) = V' => rlen (V') = 0)
RB2 : (a1 :  ref buffer) ->
State {\(h : heap). rdom (h, a1) = true}
v : ref buffer
{ \(h : heap), (v : ref buffer), (h' : heap).
\(V1 : buffer), (VN : buffer).
rdom (h', v) = true /\
(rsel (h, a1) = V1 /\
rsel (h, v) = VN ) =>

(rsel (h', a1) = rsel (h, a1) /\
[VN = V1] /\
rlen (VN) = rlen (V1) /\
rdisjoint (V1, VN) = true)

Q1 : (x : a) ->
State {\(h : heap). true}
v : ref queue
{ \(h : heap), (v : ref queue), (h' : heap).
\(Q' : queue).
qdom (h', v) = false  /\
(qsel (h', v) = Q' =>
(qlen (Q') = 1 /\
vqmem (Q', x) = true))

Q2 : (x : a) ->
(q1 : ref queue) ->
State {\(h : heap). qdom (h, q1) = true}
v : ref queue
{\(h : heap), (v : ref queue), (h' : heap).
\(Q1: queue), (Q1' : queue), (QN' : queue).
qsel (h, q1) = Q1 /\
qsel (h', q1) = Q1' /\
qsel (h', v) = QN' ) =>
(qlen (QN') == (qlen (Q1) + 1) /\
qlen (Q1') = 0 /\
not [q1 = v])

Q3 : (x : a) ->
(q : ref queue) ->
State {\(h : heap).
\(Q1 : queue).
qdom (h, q1) = true /\
qsel (h, q1) = Q1 /\
qlen (Q1, x) > 0}
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(Q1: queue), (Q1' : queue).
(qsel (h, q1) = Q1 /\
qsel (h', q1) = Q1' )
(vqmem (Q1', x) = true /\
qlen (Q1') = qlen (Q1))

ZL1 : (x : a) ->
State {\(h : heap). true}
v : {v : ref ziplist | true}
{\(h : heap), (v : ref ziplist), (h' : heap).
\(Z' : ziplist).
zsel (h, v) = Z /\
zsel (h', v) = Z' /\
zlen (Z') = 1 /\
zllen (Z') = 0/\
zrlen (Z') =  1

ZL2 : (x : a) ->
State {\(h : heap). true}
v : {v : ref ziplist | true}
{\(h : heap), (v : ref ziplist), (h' : heap).
\(Z' : ziplist).
zsel (h, v) = Z /\
zsel (h', v) = Z' /\
zlen (Z') = 1 /\
zllen (Z') = 0/\
zrlen (Z') =  1
ZL3 : (x : a) ->
(z : ref ziplist) ->
State {\(h : heap).
\(Z : ziplist). dom (h, z) = true /\
zsel (h , z) = Z /\
zlen (Z) = 1}
v : {v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(Z' : ziplist).
zsel (h, v) = Z /\
zsel (h', v) = Z' /\
zlen (Z') = 2 /\
zllen (Z') = 0/\
zrlen (Z') =  2
PQ1 :(x : int)  ->
(pqueue : ref pq) ->
State {\(h : heap).
\(P : pq).
pqdom (h, pqueue) = true}
v : { v : int | true}
{\(h : heap), (v : a), (h' : heap).
\(P : pq), (P' : pq).
pqsel (h, pqueue) = P /\
pqsel (h', pqueue) = P')
(pqmem (P', x) = true /\
(minimum (P) > x  => [v = x]) /\
(not minimum (P) > x  =>
not [v = x])
PQ2 : (pqueue : ref pq) ->
{x : int | Top > x } ->
State {\(h : heap).
\(P : pq).
pqdom (h, pqueue) = true}
v : { v : ref pq | true}
{\(h : heap), (v : ref pq), (h' : heap).
\(P' : pq).
(pqdom (h, v) = true /\
pqsel (h', v) = P')
(pqmem (P', x) = true /\
minimum (P') = x)

HT1 : (ht : ref table)  ->
(k : key) ->
(val : a) ->
State {\(h : heap).
\(H : table).
hdom (h, ht) = true /\
hsel (h, ht) = H /\
hmem (H,k) = true
v : { v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(H: table), (H' : table).
hsel (h, ht) = H /\
hsel (h', ht) = H' /\
hvmem (H', val) = true /\
hmem (H', k) = true

HT2 : (ht : ref table)  ->
(k : key) ->
(val : a) ->
State {\(h : heap).
\(H : table).
hdom (h, ht) = true /\
hsel (h, ht) = H
v : { v : unit | true}
{\(h : heap), (v : unit), (h' : heap).
\(H: table), (H' : table).
hsel (h, ht) = H /\
hsel (h', ht) = H' /\
hvmem (H', val) = true /\
hmem (H', k) = true

HT3 : (ht : ref table)
(k : key) ->
(val : a) ->
State {\(h : heap).
v : { v : ref table | true}
{\(h : heap), (v : ref table), (h' : heap).
\(H: table), (H' : table).
hsel (h, ht) = H/\
hsel (h', ht) = H' /\
hvmem (H', val) = true /\
hmem (H', k) = true

\end{lstlisting}
\clearpage
\subsubsection{Parsers Queries}
\paragraph{Queries over PNG and CDCL}
\begin{lstlisting}[basicstyle=\small\tt,breaklines=true,language=ML]
P1 :State {\(h : heap).
sel (h, fuel) == 0}
v : { v : pngtriple | true}
{\(h: heap),(v : pngtriple),(h': heap).
\(Len : int), (Tys:quad), (Con:[char]).
(pnglen (v) = Len /\
pngts (v) = Tys /\
pngdata (v) = Con) =>
sel (h', fuel) == 0 /\
Len == (length (Con) + 4)
P2: State {\(h : heap).
sel (h, fuel) == 0}
v : { v : pngtriple | true}
{\(h: heap),(v : pair),(h': heap).
\(Len : int), (Tys:quad), (Con:[char]).
(ppfst (v) = Len /\
ppsnd (v) = Con) =>
sel (h', fuel) == 0 /\
Len == (length (Con))
P3 : (x : { v : unit | true}) -> State {\(h : heap).
sel (h, fuel) == 0}
v : { v : pngtriple | true}
{\(h: heap),(v : pngtriple),(h': heap).
\(Len : int), (Tys:quad), (Con:[char]).
(pnglen (v) = Len /\
pngts (v) = Tys /\
pngdata (v) = Con) =>

sel (h', fuel) == 0 /\
(   Len == length (Con) + 4 \/
Len == length (Con) + quadlength (Tys)
P4 : State {\(h : heap).
\(Id : [string]), (Ty : [string]).
(idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
ldisjoint (Id, Ty) = true)
v : {v : tdecl | true}
{   \(h : heap), (v : typeexpr), (h' : heap).
\(Id : [string]), (Ty : [string]), (Id' : [string]), (Ty' : [string]).
(idsel (h', ids) = Id' /\
idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
tysel (h', types) = Ty') =>
ldisjoint (Id', Ty') = true)
P5 :State {\(h : heap).
\(Id : [string]), (Ty : [string]).
(idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
ldisjoint (Id, Ty) = true)
v : {v : externfun | true}
{   \(h : heap), (v : externfun), (h' : heap).
\(Id : [string]), (Ty : [string]), (Id' : [string]), (Ty' : [string]).
(idsel (h', ids) = Id' /\
idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
tysel (h', types) = Ty') =>
(ldisjoint (Id', Ty') = true)
P6 :State {\(h : heap).
\(Id : [string]), (Ty : [string]).
(idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
ldisjoint (Id, Ty) = true)
v : {v : externvar | true}
{   \(h : heap), (v : externvar), (h' : heap).
\(Id : [string]), (Ty : [string]), (Id' : [string]), (Ty' : [string]).
(idsel (h', ids) = Id' /\
idsel (h, ids) = Id /\
tysel (h, types) = Ty /\
tysel (h', types) = Ty') =>
(ldisjoint (Id', Ty') = true)

\end{lstlisting}
\end{document}


# SpectroMap: Peak detection algorithm for audio fingerprinting

Aarón López-García
Department of Computer Science
University of Valencia
<EMAIL_ADDRESS>
(September 2022)

###### Abstract

Audio fingerprinting is a technique used to identify and match audio
recordings based on their unique characteristics. It involves creating a
condensed representation of an audio signal that can be used to quickly
compare and match against other audio recordings. The fingerprinting process
involves analyzing the audio signal to extract certain features, such as
spectral content, tempo, and rhythm, among other things. In this paper, we
present SpectroMap, an open-source GitHub repository for audio fingerprinting
written in Python programming language. It is composed of a peak search
algorithm that extracts topological prominences from a spectrogram via time-
frequency bands. In this paper, we introduce the algorithm functioning with
two experimental applications in a high-quality urban sound dataset and
environmental audio recordings to describe how it works and how effective it
is in handling the input data. Finally, we have posed two Python scripts that
would reproduce the proposed case studies in order to ease the reproducibility
of our audio fingerprinting system.

##  1 Introduction

In computer science, fingerprinting is a procedure that summarizes the input
data by mapping it to a much shorter item Broder (1993). Similarly to human
fingerprints, such transformation contains the essential information and
properties of the original data, so it can be used to identify it among other
samples Wagner (1983).

Regarding the acoustic field, audio fingerprinting is understood as an
algorithm that extracts the main components taking into account the perceptual
characteristics of the audio Cano et al. (2005). Most of the time, these
techniques are applied over the spectrogram representation of the signal.
Then, the pattern extraction is conducted by means of time domain Ramalingam
and Krishnan (2006), frequency domain Seo et al. (2006), or a combination of
both called time-frequency domain Lu (2002). As far as implementation is
concerned, there are some techniques created for this purpose, although they
have their own advantages and limitations. The phase-based Arnold et al.
(2014); Wang et al. (2015) and chroma-based Kim et al. (2008); Kim and
Narayanan (2008) fingerprinting techniques are widely used. In regard to data
transformations, wavelets have been very effective in this field Baluja and
Covell (2008); Kamaladas and Dialin (2013); Jiang et al. (2019). Nonetheless,
this is not the only feature utilized for this purpose Burges et al. (2002);
Seo et al. (2005); Miller et al. (2005); Serrano et al. (2022).

For applicability purposes, Wang (2003) developed the idea of a constellation
map for Shazam Entertainment in order to implement an audio search algorithm.
Over the years, many different techniques have been developed Cano and Batlle
(2005). However, it is worth mentioning that their implementation in machine
learning tasks is very useful to reduce training costs, and so enable faster
implementations. For example, we can also find recognition of activities of
daily living via audio fingerprinting Pires et al. (2018).

This paper presents the SpectroMap algorithm for creating audio fingerprints
from a given audio signal. The method has been designed in order to deal with
both raw audio excerpts and pre-processed spectrograms. The main objective is
to cover the audio matching task because it can be considerably time-consuming
Kurth and Muller (2008); Garcia-Hernandez and Gomez-Ricardez (2019). In
essence, this paper is motivated by our previous work López-García et al.
(2022), where an in-depth example of the application of audio fingerprinting
was applied to music plagiarism.

##  2 Methodology

The algorithm presented in this paper has been designed in order to carry out
the entire process required to get the fingerprint from a given audio signal.
In this manner, we facilitate open-source software able to produce large-scale
signal processing. Depending on the objective of the user, we can use a raw
signal or an already computed spectrogram as input when initializing the
SpectroMap object. In case we decide to use raw signals, we can also include
the required parameters for the signal processing step. Thereupon, the
algorithm computes a local search to extract the topological prominences of
the given spectrogram. The architecture of SpectroMap is depicted in Figure 1.
In this section, we have detailed the two steps that perform the fingerprint
extraction of our algorithm.

Figure 1: Architecture of the SpectroMap algorithm with the two possible input
audio samples.

###  2.1 Signal processing

With the aim of implementing a fingerprint extraction for a given musical
signal $X_{t}$, we have designed an algorithm that computes a global peak
detection over the spectrogram associated to give us its constellation map.
Let $N_{FFT}$ and $N_{O}$ be the length of the Fast Fourier Transform (FFT)
window and the number of elements to overlap between segments respectively, we
first compute the spectrogram of the signal ($S_{tfa}$), by using the Hamming
window, in order to get the (time, frequency, amplitude) vectors by
considering these two parameters. Such representation contains the amplitude
spatial information to analyze. Our engine search determines whether a time-
frequency point can be considered locally relevant according to its
neighborhood. Then, the detection is processed regarding a required band. Let
$\\{T_{i}\\}_{i=1}^{n}$ and $\\{F_{j}\\}_{j=1}^{m}$ be the time and frequency
bands of the spectrogram with the amplitude of the event, we can reformulate
the spectrogram $S_{tfa}=(T_{i})_{i=1}^{n}=(F_{j})_{j=1}^{m}$ as its rows and
columns representations.

As part of the engine search, we define two windows $\phi_{T}^{d_{T}}$ and
$\phi_{F}^{d_{F}}$ to process the local pairwise comparisons with a respective
length of $d_{T}$ and $d_{F}$, whose functionality is to extract a number of
elements of the band and return the local maximum. Without limiting the
generality of the foregoing, we can mathematically describe the time-band
window mechanism with a length of $0<d_{T}\leq n$ and structure
$T_{i}=(T_{i}^{1},...,T_{i}^{n})$ as:

$\displaystyle\phi_{T}^{d_{T}}(T_{i})=\left(\max{\\{T_{i}^{k},\dots,T_{i}^{k+d_{T}}\\}}\right)_{1\leq
k\leq n-d_{T}-1},$ (1)

per each band $i\in\\{1,\dots,n\\}.$

When we group all the values we drop those elements that have an equal index
to avoid duplicates. Hence, we can group the window of each band to create the
set:

$\displaystyle\Phi_{T}^{d_{T}}=\\{\phi_{T}^{d_{T}}(T_{i})\\}_{i=1}^{n}.$ (2)

This way, we get the topologically prominent elements per each feature vector.
Owing to the equation (1), it is easy to note that even though there are
$n-d_{T}-1$ matches, the window $\phi_{T}^{d_{T}}(T_{i})$ may contain a
smaller number of elements whenever $d_{T}>2$. Depending on how restrictive we
need to be, we can proceed with just one of the bands or combine them to
create a more stringent search and distortion resistance since it is returned
only the peaks that are prominent in both directions. Finally, the algorithm
merges all the band-dependent peaks, as shown in equation (2), to give us the
total number of spatial points that determine the so-called audio fingerprint.
In Figure 2, we can see a graphical example of an audio fingerprint.

Figure 2: Example of the spectrogram of an acoustic signal with its
fingerprint stacked. The magnitudes are presented as seconds on the X-axis,
Hertz on the Y-axis, and Decibels depicted as a color map.

###  2.2 Algorithm

Our engine search, which boosts SpectroMap, processes audio signals in order
to return an output file with the (time, frequency, amplitude) peaks detected
in its spectrogram representation. Thus, it can be combined with the Mercury
software to complete an in-depth comparison between music excerpts. Figure 3
has a cursory description of the performance of SpectroMap. The algorithm
basically batches the files by means of the following steps:

  1. Step 1

Decide the window to use and set the parameters $N_{FFT}$ and $N_{O}$.

  2. Step 2

Read the audio file to get its amplitude vector and its sample rate.

  3. Step 3

Compute the spectrogram through the associated Fourier transformations.

  4. Step 4

Set a fixed window length ($d_{T}$, $d_{F}$ or both) for the pairwise
comparisons.

  5. Step 5

Choose the settings to proceed with the peak detection over a selected band or
a combination of both.

  6. Step 6

Create an identification matrix consisting of a binary matrix with the same
shape as the spectrogram with the position of the highlighted prominences.

  7. Step 7

Extract such elements and create a file with the (time, frequency, amplitude)
vectors.

Regarding step Step 5, the authors highly recommend selecting both bands to
perform the peak detection since the output is more filtered and spatially
consistent. For the remainder steps, the choice is a personal decision that
depends on the scope of the research. It is worth mentioning that the
limitations of the method depend on the functionality of the Signal module of
the SciPy library. Both installation and usage are described in our GitHub
repository López-García (2022).

Figure 3: Flowchart with the inner architecture of the algorithm implemented
that detects the topological peaks of a spectrogram.

##  3 Case study: Processing of environmental and urban sound events

The aim of this section is to present an experiment in which the performance
of SpectroMap is analyzed in terms of computational cost. To this end, we have
evaluated the speed of our algorithm over two datasets. On the one hand, Urban
Sound 8K Salamon et al. (2014) is an audio dataset that contains 8732 labeled
sound excerpts. The files are pre-sorted into ten folds in order to help in
the reproduction and comparison of machine-learning experiments. The samples
have a duration of $\approx 4$s and they are classified as urban sounds from
10 classes: air conditioner, car horn, children playing, dog bark, drilling,
engine idling, gunshot, jackhammer, siren, and street music. On the other
hand, ESC-10 Piczak (2015) is a labeled collection of 400 environmental audio
recordings suitable for benchmarking methods of environmental sound
classification. In particular, ESC-10 is a subset of the major dataset ESC-50,
which contains 2000 audio excerpts with a total size of $\approx 600$MB
publicly available (https://github.com/karolpiczak/ESC-50#download).

For both datasets, the main use commonly attached is the classification task
via supervised AI models. We can find robust performance ($94.6\%$ accuracy)
utilizing CNN architectures Hassan et al. (2019) for the Urban Sound 8K set
and other applications in low-cost monitoring devices Mydlarz et al. (2017).
For the ESC-50 dataset, and so the ESC-10, it has been shown that deep
architectures such as Transformers Chen et al. (2022) and CNNs Elizalde et al.
(2022) can learn with high precision ratios from this kind of audio sources
with $97.00\%$ and $96.70\%$ of accuracy respectively.

###  3.1 Computational costs

Table 1 is presented the computational cost associated with the audio
fingerprinting extraction task. All the experiments were produced by using
seconds as time magnitude. For both datasets, it is computed the peak
detection per folder and per audio sample. The Python script utilized to
obtain the Table 1 is displayed in A. The computer that conducted the
experiments was equipped with an AMD Ryzen 7 3700u with 16GB RAM running in
Ubuntu 20.04.3 LTS OS.

Table 1: Summary of the computational costs in seconds, with four significant figures, produced during the audio fingerprinting task for both Urban Sound 8K and ESC-10 experiments. The mean column contains the average plus-minus its standard deviation. |  |  | Processing times (seconds) |
---|---|---|---|---
| Set | Files | Min | Mean | Max | Total | it/s
Urban Sound 8K | 1 | 874 | 0.0046 | 0.0885 $\pm$ 0.0257 | 0.2414 | 77.4187 | 11.2892
2 | 889 | 0.0047 | 0.0921 $\pm$ 0.0382 | 0.2686 | 81.9609 | 10.8466
3 | 926 | 0.0055 | 0.1154 $\pm$ 0.0469 | 0.2682 | 106.9068 | 8.6617
4 | 991 | 0.0098 | 0.0957 $\pm$ 0.0385 | 0.2546 | 94.8629 | 10.4466
5 | 937 | 0.0075 | 0.0905 $\pm$ 0.0367 | 0.2606 | 84.8130 | 11.0478
6 | 824 | 0.0079 | 0.1135 $\pm$ 0.0486 | 0.2799 | 93.5494 | 8.8081
7 | 839 | 0.0064 | 0.1040 $\pm$ 0.0407 | 0.2265 | 87.3159 | 9.6087
8 | 807 | 0.0094 | 0.1114 $\pm$ 0.0502 | 0.3355 | 89.9648 | 8.9701
9 | 817 | 0.0067 | 0.1277 $\pm$ 0.0500 | 0.2494 | 104.3812 | 7.8270
10 | 838 | 0.0087 | 0.0929 $\pm$ 0.0350 | 0.2497 | 77.8634 | 10.7624
ESC-10 | 1 | 80 | 0.1256 | 0.1803 $\pm$ 0.0244 | 0.2859 | 14.4258 | 5.5456
2 | 80 | 0.1434 | 0.1775 $\pm$ 0.0127 | 0.2104 | 14.2069 | 5.6310
3 | 80 | 0.1400 | 0.1799 $\pm$ 0.0153 | 0.2335 | 14.3952 | 5.5573
4 | 80 | 0.1169 | 0.1777 $\pm$ 0.0193 | 0.2210 | 14.2231 | 5.6246
5 | 80 | 0.1082 | 0.1781 $\pm$ 0.0180 | 0.2142 | 14.2512 | 5.6135

###  3.2 Graphical representation of the outputs for the environmental sound
dataset

Figure 4: Ten samples of the ESC-10 dataset representing each of the different
acoustic classes.

When we are conducting the signal processing stage for extracting the audio
fingerprint of the audio samples, the representation of the fingerprints per
class gives us significant information about the events. In order to give an
overview of such a depiction, we have analyzed the classes of the ESC-10
dataset. The main point is to present the coordinates (time-frequency)
relevant in terms of membership. On the one hand, figure 4 contains a random
sample per each of the 10 available classes. On the other hand, once we have
carried out our algorithm, we have stored the coordinates that represent a
peak within the fingerprint per each sample as a sequence
$\\{(t_{i},f_{i})\\}_{i=1}^{I_{n}}$ so that each fingerprint contains $I_{n}$
topological prominences. With that information, we have generated a global
class fingerprint consisting of natural entries that determine the number of
times a coordinate has been selected as a peak per each sample of the same
class. With the same notation as (2), we can define the global class
fingerprint per each class $k$ as

$FP_{k}=\sum_{i=1}^{N_{k}}\mathbf{1}(\Phi_{i}^{k}),$ (3)

where the summation stands for the matrix sum operator and $\mathbf{1}$ for
the matrix characteristic function of each fingerprint, $\Phi_{i}^{k}$ for the
$i$ fingerprint of the class $k$, and $N_{k}$ for the number of elements in
the class $k$.

Considering all the mathematical notation aforementioned, Figure 5 shows each
of the $FP_{k}$ in a viridis color palette indicating that brighter colors
have a major impact in the representation of $FP_{k}$.

Figure 5: Constellation map of the different classes of the ESC-10 datasets.
Each subfigure contains the stacked aggregation per each fingerprint of each
acoustic class indicating that the brighter the color, the more relevant the
coordinate.

##  4 Discussion and future work

The SpectroMap algorithm has shown great performance when dealing with real-
world acoustic scenarios. In the case studies conducted, our algorithm took a
total of 899.03” (14’ and 59.03”) for the Urban Sound 8K and 71.50” (1’ and
1.50”) for the ESC-10 datasets. This can be summarized in a number of 9.82
$\pm$ 1.14 and 5.59 $\pm$ 0.03 iterations per second on average respectively.
Therefore, SpectroMap can be considered an effective publicly available
technique for audio fingerprinting.

One of the major advantages that arise from our experiments is that we can
efficiently process audio signals (or even many others) for further analysis.
Additionally, all the stages from which the acoustic sample is transformed are
clearly defined, thus removing any kind of black boxes. Then, from these
contributions, a potential future work would be to approach machine learning
tasks by means of distance measures or similarity functions between audio
samples. On the one hand, we could attach classification problems with a
similar strategy to the KNN algorithm Fix and Hodges (1951); Cover and Hart
(1967); Jiang et al. (2007). Basically, we would predict the class of some
audio regarding its distance to the already known fingerprints. Another
alternative would be the use of AutoEncoders Zhang et al. (2020), with a semi-
supervised approach, that reconstructs some audio Mnasri et al. (2020) from
the information of a given set of fingerprints. On the other hand, we could
perform an unsupervised strategy to determine the different sound sources
based on the distribution that they present, such as K-means Lloyd (1982) or
DBSCAN Ester et al. (1996).

Finally, it is important to remark that the use of Hertz as a frequency scale
has been used for simplicity. Our main purpose has been to introduce the
SpectroMap algorithm and show its applicability and performance. Then, we
conducted a basic signal process to convert a signal into a spectrogram.
However, there exist many choices to get different scales or units. For
instance, the Mel-scale Stevens et al. (1937); O’Shaughnessy (1987) would be a
great alternative in order to get the perceptual scale of pitches of the
events studied. A further application can be found in Shen et al. (2018) and
Jang et al. (2019)).

##  5 Conclusions

We have introduced SpectroMap, a peak detection algorithm whose main
application is the extraction of audio fingerprints. The algorithm not only
processes raw signals but also preprocessed spectrograms, which means a major
advantage in this field. Apart from a detailed explanation of the procedure
and structure of the algorithm, we have also evaluated its performance in
state-of-the-art datasets for audio analysis. It has been shown that
SpectroMap is an effective and fast algorithm with an average of 1.340 and
3.336 iterations per second for the datasets presented in the case study
(Urban Sound 8K and ESC-10). Further interpretations and representations have
been shown in order to give a better understanding of the outputs of our
algorithm. The code and Python implementation of the package has been
presented in a straightforward manner in order to ease applicability and
reproducibility. Even though we have not emphasized the underlying application
on audio signals with comparison purposes, an instance of such an application
can be found in our last paper López-García et al. (2022).

## References

  * Broder [1993] A. Broder.  Some applications of rabin’s fingerprinting method.  _Sequences II_ , pages 143–152, 1993.
  * Wagner [1983] Neal R. Wagner.  Fingerprinting.  In _1983 IEEE Symposium on Security and Privacy_ , pages 18–18, 1983\.  doi: 10.1109/SP.1983.10018.
  * Cano et al. [2005] P. Cano, E. Batlle, E. Gómez, L. C. T. Gomes, and M. Audio Fingerprinting: Concepts and Bonnet.  Applications.  _Computational Intelligence For Modelling And Prediction_ , pages 233–245, 2005.  URL https://link.springer.com/chapter/10.1007/10966518_17.
  * Ramalingam and Krishnan [2006] Arunan Ramalingam and Sridhar Krishnan.  Gaussian mixture modeling of short-time fourier transform features for audio fingerprinting.  _IEEE Transactions on Information Forensics and Security_ , 1(4):457–463, 2006.  doi: 10.1109/TIFS.2006.885036.
  * Seo et al. [2006] J.S. Seo, Minho Jin, Sunil Lee, Dalwon Jang, Seungjae Lee, and C.D. Yoo.  Audio fingerprinting based on normalized spectral subband moments.  _IEEE Signal Processing Letters_ , 13(4):209–212, 2006.  doi: 10.1109/LSP.2005.863678.
  * Lu [2002] Chun-Shien Lu.  Audio fingerprinting based on analyzing time-frequency localization of signals.  In _2002 IEEE Workshop on Multimedia Signal Processing._ , pages 174–177, 2002.  doi: 10.1109/MMSP.2002.1203275.
  * Arnold et al. [2014] M. Arnold, X. Chen, P. Baum, U. Gries, and G. A Doërr.  Phase-based audio watermarking system robust to acoustic path propagation.  _IEEE Transactions On Information Forensics And Security_ , 9:411–425, 2014.
  * Wang et al. [2015] X. Wang, L. Gao, and S. PhaseFi Mao.  Phase fingerprinting for indoor localization with a deep learning approach.  _2015 IEEE Global Communications Conference (GLOBECOM)_ , pages 1–6, 2015.
  * Kim et al. [2008] S. Kim, E. Unal, and S. Narayanan.  Music fingerprint extraction for classical music cover song identification.  _2008 IEEE International Conference On Multimedia And Expo_ , pages 1261–1264, 2008.
  * Kim and Narayanan [2008] S. Kim and S. Narayanan.  Dynamic chroma feature vectors with applications to cover song identification.  _2008 IEEE 10th Workshop On Multimedia Signal Processing_ , pages 984–987, 2008.
  * Baluja and Covell [2008] Shumeet Baluja and Michele Covell.  Waveprint: Efficient wavelet-based audio fingerprinting.  _Pattern Recognition_ , 41(11):3467–3480, 2008\.  ISSN 0031-3203.  doi: https://doi.org/10.1016/j.patcog.2008.05.006.  URL https://www.sciencedirect.com/science/article/pii/S0031320308001702.
  * Kamaladas and Dialin [2013] M.Davidson Kamaladas and M.Maxina Dialin.  Fingerprint extraction of audio signal using wavelet transform.  In _2013 International Conference on Signal Processing , Image Processing & Pattern Recognition_, pages 308–312, 2013.  doi: 10.1109/ICSIPR.2013.6497946.
  * Jiang et al. [2019] Yuantao Jiang, Chunxue Wu, Kaifa Deng, and Yan Wu.  An audio fingerprinting extraction algorithm based on lifting wavelet packet and improved optimal-basis selection.  _Multimedia Tools and Applications_ , 78(21):30011–30025, Nov 2019.  ISSN 1573-7721.  doi: 10.1007/s11042-018-6802-y.  URL https://doi.org/10.1007/s11042-018-6802-y.
  * Burges et al. [2002] Christopher J. C. Burges, John C. Platt, and Soumya Jana.  Extracting noise-robust features from audio data.  In _2002 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , volume 1, pages I–1021–I–1024, 2002.  doi: 10.1109/ICASSP.2002.5743968.
  * Seo et al. [2005] J.S. Seo, Minho Jin, Sunil Lee, Dalwon Jang, Seungjae Lee, and C.D. Yoo.  Audio fingerprinting based on normalized spectral subband centroids.  In _Proceedings. (ICASSP ’05). IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005._ , volume 3, pages iii/213–iii/216 Vol. 3, 2005.  doi: 10.1109/ICASSP.2005.1415684.
  * Miller et al. [2005] Matthew L. Miller, Manuel Acevedo Rodriguez, and Ingemar J. Cox.  Audio fingerprinting: Nearest neighbor search in high dimensional binary spaces.  _Journal of VLSI signal processing systems for signal, image and video technology_ , 41(3):285–291, Nov 2005.  ISSN 0922-5773.  doi: 10.1007/s11265-005-4152-2.  URL https://doi.org/10.1007/s11265-005-4152-2.
  * Serrano et al. [2022] Salvatore Serrano, Murtadha Arif Bin Sahbudin, Chakib Chaouch, and Marco Scarpa.  A new fingerprint definition for effective song recognition.  _Pattern Recognition Letters_ , 160:135–141, 2022.  ISSN 0167-8655.  doi: https://doi.org/10.1016/j.patrec.2022.06.009.  URL https://www.sciencedirect.com/science/article/pii/S0167865522002033.
  * Wang [2003] A. Wang.  An industrial strength audio search algorithm.  _Ismir_ , 2003:7–13, 2003.
  * Cano and Batlle [2005] P. Cano and E. A Batlle.  Review of audio fingerprinting.  _Journal Of VLSI Signal Processing_ , 41:271–284, 2005\.
  * Pires et al. [2018] I. Pires, R. Santos, N. Pombo, N. Garcia, F. Flórez-Revuelta, S. Spinsante, R. Goleva, and E. Zdravevski.  Recognition of activities of daily living based on environmental analyses using audio fingerprinting techniques: A systematic review.  _Sensors_ , 18, 2018.  URL https://www.mdpi.com/1424-8220/18/1/160.
  * Kurth and Muller [2008] F. Kurth and M. Muller.  Efficient index-based audio matching.  _IEEE Transactions On Audio, Speech, And Language Processing_ , 16:382–395, 2008.
  * Garcia-Hernandez and Gomez-Ricardez [2019] Jose Juan Garcia-Hernandez and Juan Jose Gomez-Ricardez.  Hardware architecture for an audio fingerprinting system.  _Computers & Electrical Engineering_, 74:210–222, 2019\.  ISSN 0045-7906.  doi: https://doi.org/10.1016/j.compeleceng.2019.01.017.  URL https://www.sciencedirect.com/science/article/pii/S0045790618303537.
  * López-García et al. [2022] A. López-García, B. Martínez-Rodríguez, and V. A Liern.  Proposal to compare the similarity between musical products.  _One More Step for Automated Plagiarism Detection?._ , pages 192–204, 2022.
  * López-García [2022] A. López-García.  Spectromap is a peak detection algorithm that computes the constellation map for a given signal.  _GitHub Repository_ , 1:1, 2022.  URL https://pypi.org/project/spectromap/.
  * Salamon et al. [2014] J. Salamon, C. Jacoby, and J. A Bello.  Dataset and taxonomy for urban sound research.  _Proceedings Of The 22nd ACM International Conference On Multimedia_ , pages 1041–1044, 2014.
  * Piczak [2015] K. Esc Piczak.  Dataset for environmental sound classification.  _Proceedings Of The 23rd Annual ACM Conference On Multimedia_ , 2015:1015–1018, 2015.  URL http://dl.acm.org/citation.cfm?doid=2733373.2806390.
  * Hassan et al. [2019] S. Hassan, Zeeshan Khan, and Ghani Khan M.  M. & saleem, s. robust sound classification for surveillance using time frequency audio features.  _2019 International Conference On Communication Technologies (ComTech)_ , pages 13–18, 2019.
  * Mydlarz et al. [2017] C. Mydlarz, J. Salamon, and J. Bello.  The implementation of low-cost urban acoustic monitoring devices.  _Applied Acoustics_ , 117:207–218, 2017.  URL https://www.sciencedirect.com/science/article/pii/S0003682X1630158X.
  * Chen et al. [2022] K. Chen, X. Du, B. Zhu, Z. Ma, T. Berg-Kirkpatrick, and S. Hts-at: A Dubnov.  Hierarchical token-semantic audio transformer for sound classification and detection.  _CoRR_ , abs/2202.00874, 2022.  URL https://arxiv.org/abs/2202.00874.
  * Elizalde et al. [2022] B. Elizalde, S. Deshmukh, M. Ismail, and H. Clap Wang.  Learning audio concepts from natural language supervision. (arxiv.  _arXiv_ , 2022, 2022.  URL https://arxiv.org/abs/2206.04769.
  * Fix and Hodges [1951] E. Fix and J. Discriminatory Analysis Hodges.  Nonparametric discrimination: Consistency properties.  _International Statistical Review / Revue Internationale De Statistique_ , 57:238–247, 1951.  URL http://www.jstor.org/stable/1403797.
  * Cover and Hart [1967] T. Cover and P. Hart.  Nearest neighbor pattern classification.  _IEEE Transactions On Information Theory_ , 13:21–27, 1967\.
  * Jiang et al. [2007] L. Jiang, Z. Cai, D. Wang, and S. Jiang.  Survey of improving k-nearest-neighbor for classification.  _Fourth International Conference On Fuzzy Systems And Knowledge Discovery (FSKD 2007)_ , 1:679–683, 2007.
  * Zhang et al. [2020] G. Zhang, Y. Liu, and X. A Jin.  survey of autoencoder-based recommender systems.  _Frontiers Of Computer Science_ , 14:430–450, 2020.  URL https://doi.org/10.1007/s11704-018-8052-6.
  * Mnasri et al. [2020] Z. Mnasri, S. Rovetta, and F. Masulli.  Audio surveillance of roads using deep learning and autoencoder-based sample weight initialization.  _2020 IEEE 20th Mediterranean Electrotechnical Conference ( MELECON)_ , pages 99–103, 2020.
  * Lloyd [1982] S. Lloyd.  Least squares quantization in pcm.  _IEEE Transactions On Information Theory_ , 28:129–137, 1982.
  * Ester et al. [1996] M. Ester, H. Kriegel, J. Sander, and X. A Xu.  Density-based algorithm for discovering clusters in large spatial databases with noise.  _Proceedings Of The Second International Conference On Knowledge Discovery And Data Mining_ , pages 226–231, 1996.
  * Stevens et al. [1937] S. S. Stevens, J. Volkmann, and E. B. Newman.  A scale for the measurement of the psychological magnitude pitch.  _The Journal of the Acoustical Society of America_ , 8(3):185–190, 1937.  doi: 10.1121/1.1915893.
  * O’Shaughnessy [1987] Douglas O’Shaughnessy.  _Speech communication: human and machine_ , volume 2.  Wiley-IEEE Press, 1987.  ISBN 0780334493.  URL https://ieeexplore.ieee.org/book/5263814.
  * Shen et al. [2018] Jonathan Shen, Ruoming Pang, Ron J. Weiss, Mike Schuster, Navdeep Jaitly, Zongheng Yang, Zhifeng Chen, Yu Zhang, Yuxuan Wang, Rj Skerrv-Ryan, Rif A. Saurous, Yannis Agiomvrgiannakis, and Yonghui Wu.  Natural tts synthesis by conditioning wavenet on mel spectrogram predictions.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4779–4783, 2018.  doi: 10.1109/ICASSP.2018.8461368.
  * Jang et al. [2019] Byeong-Yong Jang, Woon-Haeng Heo, Jung-Hyun Kim, and Oh-Wook Kwon.  Music detection from broadcast contents using convolutional neural networks with a mel-scale kernel.  _EURASIP Journal on Audio, Speech, and Music Processing_ , 2019(1):11, Jun 2019.  ISSN 1687-4722.  doi: 10.1186/s13636-019-0155-y.  URL https://doi.org/10.1186/s13636-019-0155-y.

##  Appendix A Python implementation

This section is dedicated to the application of the SpectroMap algorithm to
some kind of example. In particular, the module is designed to process either
a raw signal or a spectrogram. For the first case, we make use of the
`spectromap` object. For the second case, we apply the `peak_search` function.
In addition, the script that reproduces the results shown is displayed at the
end in order to easy reproducibility purposes.

The library was written with the Python 3.8 version and its usage depends just
on NumPy 1.19 and SciPy 1.6.3. packages. The repository is under the GNU
General Public License v3.0.

###  A.1 Application over a raw signal

import numpy as np from spectromap.functions.spectromap import spectromap

# Generate a random signal y = np.random.rand(44100) kwargs = ’fs’: 22050,
’nfft’: 512, ’noverlap’:64

# Instantiate the SpectroMap object SMap = spectromap(y, **kwargs)

# Get the spectrogram representation plus its time and frequency bands f, t, S
= SMap.get_spectrogram()

# Extract the topological prominent elements from the spectrogram. #
Coordinates matrix as (time, freq) # Peak matrix. fraction = 0.15 # Fraction
of spectrogram to compute local comparisons condition = 2 # Axis to analyze
(0: Time, 1: Frequency, 2: Time+Frequency) id_peaks, peaks =
SMap.peak_matrix(fraction, condition)

# Get the peaks coordinates as as (s, Hz, dB)-array. extraction_t_f_dB =
SMap.from_peaks_to_array()

###  A.2 Application over a given spectrogram

from spectromap.functions.spectromap import peak_search

fraction = 0.05 # Fraction of spectrogram to compute local comparisons
condition = 2 # Axis to analyze (0: Time, 1: Frequency, 2: Time+Frequency)
id_peaks, peaks = peak_search(spectrogram, fraction, condition)

###  A.3 Application over a dataset

Here is presented the script that reproduces the experimental cases for the
Urban Sound 8K described in section 3

# OS packages import os import time # Parallel computing from pqdm.processes
import pqdm

# Sinal processing modules import numpy as np import librosa from
spectromap.functions.spectromap import spectromap

NFFT = 1024 FRACTION = 1/3 # Fraction of spectrogram to compute local
comparisons CONDITION = 2 # Axis to analyze (0: Time, 1: Frequency, 2:
Time+Frequency)

def spectromap_for_urbansound8K(file): # Transform file file =
os.path.join(folder_path, file) # Read file signal, sample_rate =
librosa.load(file, mono=True) kwargs = ’fs’: sample_rate, ’nfft’: NFFT #
Initialize SpectroMap t0 = time.time() SMap = spectromap(signal, **kwargs) #
Get the spectrogram representation plus its time and frequency bands f, t, S =
SMap.get_spectrogram() # Extract the topological prominent elements from the
spectrogram id_peaks, peaks = SMap.peak_matrix(FRACTION, CONDITION) # Get the
peaks coordinates as as (s, Hz, dB)-array. extraction_t_f_dB =
SMap.from_peaks_to_array() cost_time = time.time() - t0 return cost_time

# Path to Urban Sound 8K audio excerpts path_UrbanSound8K = os.path.join(’./’,
’UrbanSound8K/audio’)

computational_costs = [] FOLDERS = os.listdir(path_UrbanSound8K) for folder in
FOLDERS: print(’: ’.format(folder)) folder_path =
os.path.join(path_UrbanSound8K, folder) # Start process t0 = time.time()
folder_time = pqdm(os.listdir(folder_path), spectromap_for_urbansound8K,
n_jobs = os.cpu_count()) folder_time = np.array(folder_time) # End process
computational_costs.append(folder_time)

np.save(’UrbanSound8K_computational_costs.npy’, computational_costs)


# Conflict-Averse Gradient Optimization of Ensembles for Effective Offline
Model-Based Optimization

Sathvik Kolli
UC Berkeley Electrical Engineering & Computer Science
<EMAIL_ADDRESS>
###### Abstract

Data-driven offline model-based optimization (MBO) is an established practical
approach to black-box computational design problems for which the true
objective function is unknown and expensive to query. However, the standard
approach which optimizes designs against a learned proxy model of the ground
truth objective can suffer from distributional shift. Specifically, in high-
dimensional design spaces where valid designs lie on a narrow manifold, the
standard approach is susceptible to producing out-of-distribution, invalid
designs that “fool” the learned proxy model into outputting a high value.
Using an ensemble rather than a single model as the learned proxy can help
mitigate distribution shift, but naive formulations for combining gradient
information from the ensemble, such as minimum or mean gradient, are still
suboptimal and often hampered by non-convergent behavior.

In this work, we explore alternate approaches for combining gradient
information from the ensemble that are robust to distribution shift without
compromising optimality of the produced designs. More specifically, we explore
two functions, formulated as convex optimization problems, for combining
gradient information: multiple gradient descent algorithm (MGDA) [?] and
conflict-averse gradient descent (CAGrad) [?]. We evaluate these algorithms on
a diverse set of five computational design tasks [?]. We compare performance
of ensemble MBO with MGDA and ensemble MBO with CAGrad with three naive
baseline algorithms: (a) standard single-model MBO, (b) ensemble MBO with mean
gradient, and (c) ensemble MBO with minimum gradient.

Each algorithm produces 128 optimized designs, and we report performance of
these designs under three metrics: (a) max ground truth score, (b) average
ground truth score, (c) 50th percentile ground truth score.

  * •

For the max ground truth score, we find that MGDA is in the top 2 best-
performing algorithms on 4 tasks and is the best-performing algorithm on 2
tasks. CAGrad is the in the top 2 best-performing algorithms on 2 tasks.

  * •

For the 50th percentile ground truth score, MGDA is the best-performing
algorithm on 1 task. CAGrad is in the top 2 best-performing algorithms on 3
tasks and is the best-performing algorithm on 2 tasks.

  * •

Finally, for the average ground truth score, MGDA is the best-performing
algorithm on 2 tasks. CAGrad is in the top 2 best-performing algorithms on 3
tasks and is the best-performing algorithm on 2 tasks.

Our results suggest that MGDA and CAGrad strike a desirable balance between
conservatism and optimality. In general, we noticed that MGDA and CAGrad
performed equally well, if not better, than other algorithms on the max ground
truth score. However, both algorithms lead to significant improvement on the
average and 50th percentile ground truth scores, suggesting they may be more
conservative and less susceptible to being “fooled” by invalid designs. Our
results demonstrate that MGDA and CAGrad can help robustify data-driven
offline MBO without compromising optimality of designs.

##  1 Introduction

We study the problem of computational design, which arises in settings ranging
from synthetic biology to robot design. Specifically, we focus on the setting
of black-box optimization, which attempts to generate optimal designs where
the objective function and constraints are unknown. Put simply, we want to
find the optimal design, $x$, that maximizes some unknown objective function,
$f(x)$:

$\operatorname*{arg\,max}_{x}\enspace f(x)$

Examples of black-box optimization problems include optimizing robot
morphologies, biological sequences (proteins, genes), computer chips, neural
network architectures, or superconducting materials.

###  1.1 Offline Model Based Optimization (MBO)

One promising approach to solving black-box optimization problems is data-
driven MBO, where a proxy model of the unknown objective function is learned
from empirically collected data and used to guide the design procedure.

In order to model the true objective function with high fidelity, it is often
critical to actively collect additional data during the training procedure
[?]. However, in many design problems, active real-world data collection is
expensive (e.g. requires synthesizing protein structures for protein
optimization or building and testing a robot for robot design) or dangerous
(e.g. when optimizing over aircraft designs). Thus, we focus instead on the
more practical setting of offline MBO, where we are given a static dataset of
designs and cannot make any queries to the ground truth.

In essence, when we use offline MBO to solve black-box optimization problems,
we are trying to solve the problem

$\operatorname*{arg\,max}_{x}\enspace f(x)$

with two key assumptions:

  1. 1.

Black-box assumption: $f(x)$ is an unknown function

  2. 2.

Offline assumption: $f(x)$ is expensive to query

The general offline MBO workflow is illustrated in Figure 1.

Figure 1: A typical offline MBO workflow [?]. We are given a static dataset of
designs, which we use to learn a proxy model, $\hat{f}_{\theta}$, of the true
objective. Then, our design procedure is guided by the learned proxy model.

###  1.2 Distribution Shift in MBO

The most basic approach to offline, data-driven model-based optimization
involves the following steps [?]:

  1. 1.

We have a static dataset $D$ of input designs and their corresponding
objective values:

$\\{(x_{1},y_{1}),\ldots,(x_{N},y_{N})\\}$

We assume this paired data was generated from a true, unknown objective
function, $y=f(x)$.

  2. 2.

Using the dataset $D$, we learn a proxy model $\hat{f}_{\theta}(x)$ of the
unknown objective function $f(x)$, via supervised regression on the training
dataset.

  3. 3.

Finally, we find an optimal generated design $x^{*}$, by optimizing data point
$x\in D$ against the learned model $\hat{f}_{\theta}(x)$ (e.g. using $T$
gradient ascent/descent steps on the learned function):

$x_{k+1}\leftarrow x_{k}+\alpha\nabla_{x}f_{\theta}(x)|_{x=x_{k}},\text{ for
}k\in[1,T]$

As described, the above approach does not perform well in high-dimensional
input spaces, where the space of valid input designs lie on a narrow manifold,
because overestimation errors in the proxy model $\hat{f}_{\theta}(x)$ would
cause to the optimization procedure in step (3) to produce out-of-
distribution, invalid, and low-valued designs. Consequently, for the above
method to work, it is critical that we ensure that the proxy model,
$f_{\theta}(x)$, does not overestimate the objective value of out-of-
distribution points.

Some existing approaches to prevent overestimation of out-of-distribution
inputs include generative modeling, explicit density estimation, or
regularization techniques that incentivize conservatism in regions with
limited data (known as Conservative Objective Models [?]).

###  1.3 Ensemble MBO

One simple approach to addressing the issue of distribution shift is to use an
ensemble rather than a single model as a proxy for the true objective
function. The rationale behind this is that it is less likely for multiple
different models to be “fooled” by the same out-of-distribution input than it
is for a single model to suffer from overestimation errors.

The ensemble can have individual models with varied architectures and
regularization techniques. Although we don’t study this in this project,
future work can try including density models and/or Conservative Objective
Models with varying levels of conservatism in the ensemble.

When we used a single proxy model, we used gradient ascent on the model in
order to generate new designs:

$x_{k+1}\leftarrow x_{k}+\alpha\nabla_{x}\hat{f}_{\theta}(x_{k})$

Now, that we use an ensemble $\\{\hat{f}_{1}(x),\ldots,\hat{f}_{m}(x)\\}$, we
need to update our design using gradient information from all the models in
the ensemble. Thus, we get the following update:

$x_{k+1}\leftarrow x_{k}+\alpha
g(\nabla_{x}\hat{f}_{1}(x_{k}),\ldots,\nabla_{x}\hat{f}_{m}(x_{k}))$

where $g$ is some function of all the gradients.

Two naive approaches for the function $g$ are:

  * •

Mean Gradient:

$\nabla_{x}\left(\sum_{i=1}^{m}f_{i}(x)\right)$

The benefit of this approach is that it captures gradient information from all
the models in the ensemble in each gradient step.

However, while this approach may lead to optimization that is more robust to
distribution shift, it is still possible for the optimization to be “fooled”
by an out-of-distribution input, particularly when a model or group of models
within the ensemble dominate the update. Furthermore, it is also possible for
the optimization to get stuck and fail to optimize further (sometimes becoming
stuck in a oscillating manner) due to conflicting gradients.

  * •

Minimum Gradient:

$\nabla_{x}\min\left(f_{1}(x),f_{2}(x),\ldots,f_{m}(x)\right)$

The benefit of this approach is that it is conservative and is less likely to
be “fooled” by an out-of-distribution input. In fact, we can interpret the
minimum gradient update as searching for a design for which all the models in
the ensemble assign a high score. By definition, this algorithm would produce
an out-of-distribution point if and only if every model in the ensemble
overestimated the value of that point.

The drawback of this method is that it has poor convergence guarantees and is
susceptible to oscillatory behavior.

###  1.4 Multiple Gradient Descent Algorithm (MGDA) and Conflict-Averse
Gradient Descent (CAGrad)

We consider two alternative functions for $g$ that combine the gradient
information from the models in the ensemble: multiple gradient descent
algorithm (MGDA) [?] and conflict-averse gradient descent (CAGrad) [?]. Both
functions are formulated as convex optimization problems. MGDA was developed
in the setting of multi-objective optimization, while CAGrad was developed for
multi-task learning.

For MGDA, the gradient is defined in terms of the following convex
optimization problem:

$\max_{d}\enspace\min_{i}\enspace\langle d,g_{i}\rangle-\frac{1}{2}\left\lVert
d\right\rVert^{2}$

where

$g_{i}=\nabla_{x}\hat{f}_{i}(x),\enspace i=1,\ldots,m$

For CAGrad, the gradient is defined in terms of the following convex
optimization problem:

$\max_{d}\enspace\min_{i}\enspace\langle
d,g_{i}\rangle\enspace\text{s.t.}\enspace\left\lVert d-g_{0}\right\rVert\leq
c\left\lVert g_{0}\right\rVert$

where

$g_{i}=\nabla_{x}\hat{f}_{i}(x),\enspace i=1,\ldots,m$

$c\in[0,1)$ is a hyper-parameter, and
$g_{0}=\frac{1}{m}\nabla_{x}\hat{f}_{i}(x)$ is the average gradient.

The high-level intuition behind these methods is that they search for a design
with a high model-predicted objective value, while leveraging the worst local
improvement of the models in the ensemble to regularize the optimization
trajectory.

More specifically, assume we update our design $x$ by $x\leftarrow x+\alpha
d$, where $\alpha$ is a step size and $d$ is the update vector. Then, the
minimum improvement rate across the models in the ensemble is given by:

$\displaystyle R(x,d)$
$\displaystyle=\min_{i\in\\{1,\ldots,m\\}}\left(\frac{1}{\alpha}(f_{i}(x+\alpha
d)-f_{i}(x))\right)$ $\displaystyle\approx\min_{i\in\\{1,\ldots,m\\}}\langle
g_{i},d\rangle$

where we use the first-order Taylor approximation, assuming $\alpha$ is small.

We can view both MGDA and CAGrad as looking for the “best” update vector $d$
within a local ball, where we define the “best” update vector as the one that
maximizes the worst improvement rate. The difference between the two is the
local ball within which we search for the optimal update vector, where MGDA is
centered at zero, while CAGrad is centered at the average gradient $g_{0}$.

Figure 2: Comparison of methods from “Conflict-Averse Gradient Descent for
Multi-task Learning” [?]. We compare naive gradient descent (GD, top) to MGDA
(middle) and CAGrad (bottom).

In theory, MGDA provably converges to an arbitrary point on the Pareto set
[?], while CAGrad provably converges to a stationary point of the average
proxy objective, when $0\leq c<1$ [?].

The differences between mean gradient, MGDA, and CAGrad are illustrated in
Figure 2, taken from [?].

##  2 Methods

We evaluate the performance of five algorithms

  * •

Single Proxy Model, Naive Gradient Ascent

  * •

Ensemble Proxy Model, Mean Gradient Ascent

  * •

Ensemble Proxy Model, Min Gradient Ascent

  * •

Ensemble Proxy Model, MGDA

  * •

Ensemble Proxy Model, CAGrad

on five diverse benchmark tasks selected from [?]. The tasks and relevant
details are listed in Table 1.

|  | Total
---
Dataset
Size
| MBO
---
Dataset
Size
Dimensions | Type |  | Oracle (Spear.
---
Correlation)
| Primal or
---
Dual? (# of
decision var.)
TF Bind 8 | 65536 | 32898 | (8, 4) | Discrete | Lookup Table | Primal (32)
TF Bind 10 | 1048576 | 50000 | (10, 4) | Discrete | Lookup Table | Primal (40)
ChEMBL | 1093 | 546 | (31, 591) | Discrete | Random Forest (0.792) | Dual (6)
Hopper Controller | 3200 | 3200 | 5126 | Continuous | Exact | Dual (6)
Ant Morphology | 25009 | 15005 | 60 | Continuous | Exact | Primal (60)

|  | Grad. Asc.
---
Parameters
(# of rounds, $\alpha$)
| Design Model
---
Architecture
(# of parameters)
| Design Model
---
Results (Val. Spear.,
Val. Loss)
| CAGrad
---
Parameter
TF Bind 8 | (200, 10) | FullyConnected (1140801) | (0.43, 0.15) | c=0.5
TF Bind 10 | (200, 50) | FullyConnected (2122753) | (0.97, 0.00) | c=0.5
ChEMBL | (200, 100) | FullyConnected (6582401) | (0.77, 0.13) | c=0.5
Hopper Controller | (200, 1) | FullyConnected (5252097) | (0.87, 0.28) | c=0.3
Ant Morphology | (200, 0.03) | FullyConnected (8460289) | (0.57, 0.33) | c=0.2

Table 1: Comparison of the TF Bind 8, TF Bind 10, ChEMBL, Hopper, and Ant
Morphology tasks.

###  2.1 General Procedure

Each task has a corresponding dataset, which we refer to as the “total
dataset” for the task. In order to evaluate our MBO algorithms, we take a
subset of the total dataset, which we call the “MBO dataset.”

The general procedure we use for evaluation is as follows:

  1. 1.

The design/proxy model(s) are trained on the MBO dataset, which is a subset
(bottom K% of target values) of the total dataset for the task.

  2. 2.

We take the top 128 (i.e. highest target values) inputs in the MBO dataset and
optimize each one using our algorithm to produce 128 designs.

  3. 3.

Next, we calculate the (a) max ground-truth score, (b) the 50th percentile
ground-truth score, and (c) the average ground-truth score of the 128 designs.

For all the ensemble proxy model methods, we use six ensemble models, which
all have the same architecture, but are trained and validated on different
subsets of the MBO dataset.

The motivation for the above approach is that, in real design scenarios, it is
usually impractical to empirically test and use every design produced by a
given algorithm. Instead, we would only select the top-scored designs
according to our algorithm and experimentally validate them and use them.
Thus, ideally, we want an algorithm whose best set of designs is actually the
best under the true objective function.

###  2.2 Oracles

For the third step in the above procedure, some tasks have an exact oracle,
meaning that the ground truth for every possible permutation of inputs is
provided in a lookup table or that individual designs are cheap to evaluate,
while other tasks use a learned oracle (e.g. neural network, random forest
model) as a proxy for the ground truth. In this case, we train a separate
oracle model on the total dataset for the task. Details on the oracles for
each task are listed in Table 1.

###  2.3 Dual Formulations

When we use the primal formulation of MGDA and CAGrad, which are presented
above, the number of decision variables is equal to the dimensionality of the
input designs. In some cases, this is computationally feasible. However, for
some of our tasks, the inputs are high-dimensional, and we need to use the
dual formulation of MGDA and CAGrad, where the number of decision variables is
equal to the number of models in the ensemble (i.e. 6). For both MGDA and
CAGrad, since the primal problem is convex and Slater’s condition holds, we
have strong duality.

As before, we define:

$g_{i}=\nabla_{x}\hat{f}_{i}(x),\enspace i=1,\ldots,m$

$c\in[0,1)$ is a hyper-parameter for CAGrad, and
$g_{0}=\frac{1}{m}\nabla_{x}\hat{f}_{i}(x)$ is the average gradient.

The dual formulation of MGDA is

$\min_{w}\frac{1}{2}\left\lVert\sum_{i=1}^{K}w_{i}g_{i}\right\rVert^{2}\enspace\text{s.t.}\enspace\sum_{i=1}^{K}w_{i}=1\enspace\text{and}\enspace\forall
i,w_{i}\geq 0$

[?].

The dual formulation of CAGrad is

$\min_{w}g_{w}^{T}g_{0}+\sqrt{\phi}||g_{w}||\enspace\text{s.t.}\enspace\sum_{i=1}^{K}w_{i}=1\enspace\text{and}\enspace\forall
i,w_{i}\geq 0$

where $g_{w}=\sum_{i}w_{i}g_{i}$ and $\phi=c^{2}||g_{0}||^{2}$. The optimal
update vector is $d^{*}=g_{0}+g_{w^{*}}/\lambda^{*}$, where
$\lambda^{*}=\left\lVert g_{w^{*}}\right\rVert/\sqrt{\phi}$, [?].

We can interpret the dual formulations as finding the weights for a weighted
average of the gradients from each model in the ensemble.

Details regarding which formulation we use for each benchmark task are listed
in Table 1.

###  2.4 Gradient Ascent Procedure

For discrete tasks, we perform gradient updates in one-hot space. After each
update, we map the sequence back to discrete space by taking an argmax.

For continuous tasks, we normalize the inputs, so that each position has zero
mean and unit variance. We then perform gradient updates in this space.

###  2.5 Hyperparameter Selection

Each task has three relevant hyperparameters:

  * •

Design Model Architecture and Training Parameters

  * •

Gradient Ascent Parameters: (Number of Gradient Update Steps, Learning Rate)

  * •

CAGrad Hyperparameter $c$

In order to select and train the design model architecture, we randomly split
the MBO Dataset into a training set and validation set, and we selected the
model with the best validation spearman correlation and validation mean-
squared loss.

For the gradient ascent parameters, we fixed the number of gradient update
steps to 200, and we selected the learning rate by visually analyzing
optimization trajectories (i.e. plot of proxy model(s) prediction vs. number
of rounds of mutation). We selected the CAGrad hyperparameter $c$ in a similar
fashion.

It is important that hyperparameter tuning should be done purely offline,
without any access to the ground truth objective or oracle.

Details regarding the hyperparameters we use for each task are listed in Table
1.

###  2.6 Discrete Benchmark Tasks

We detail the three discrete tasks on which we perform evaluation.

The TF Bind 8 task [?] is based on an empirical dataset of measurements of
binding activity between a variety of human transcription factors and every
possible length-8 DNA sequence. The optimization goal is to identify DNA
sequences that maximizing the binding activity score for each TF. The design
space for sequences is comprised of four categorical variables, one
representing each nucleotide (A, T, C, or G). The oracle for TF Bind 8 is
exact (a lookup table containing ground-truth values for every possible
permutation of inputs).

The TF Bind 10 task [?] is a neural network produced dataset of predicted
estimates of the relative binding affinities between all unique length-10 DNA
sequences and each of two protein targets. The optimization goal is to
identify DNA sequences that maximize the predicted binding affinity to
targets. The design space for sequences is A, T, C, and G as before. The
oracle for TF Bind 10 is exact (a lookup table containing ground-truth values
for every possible permutation of inputs).

The ChEMBL task [?] is a dataset of pairs of molecules and assays which test
for specific functional properties of those molecules. The optimization goal
is to design a molecule that achieves a high functional property score on a
specific assay. The design space for molecules is based on SMILES encodings
(rather than amino acids), resulting in a design space of categorical
variables that take one of 591 values, for sequences of length 31. The oracle
for ChEMBL is a random forest, which achieves a spearman correlation of 0.792
on the dataset.

###  2.7 Continuous Benchmark Tasks

We detail the two continuous tasks on which we perform evaluation.

The Hopper Controller task is an OpenAI gym locomotion task [?]. The
optimization goal is to design a set of weights for a controller neural
network (representing a policy) that will optimize for expected return on the
locomotion task. While Hopper is typically a reinforcement learning task, we
formulate it as offline MBO by utilizing a supervised dataset of neural
network controlled weights matched with expected return values. There are 5126
continuous variables corresponding to the flattened weights of this neural
network. In order to evaluate the ground truth score for a design, we simply
load in the weights of the neural network controller and run 1000 steps of
simulation in the MuJoCo simulator [?] used with this environment.

The Ant Morphology task is an OpenAI gym task [?]. The goal is to optimize the
morphology (e.g. size, orientation, location of limbs) of Ant, a simulated
robot whose goal is to run fast (i.e. a locomotion task) in its environment.
There are 60 continuous variables corresponding to these morphological
parameters. We obtain a design’s ground truth score by running robotic
simulation in the MuJoCo simulator [?] for 100 time steps, averaging 16
independent trials.

| TF Bind 8 | TF Bind 10 |  | ChEMBL
---
(Random Forest Oracle)
| Hopper
---
Controller
| Ant
---
Morphology
dataset | 0.439 | 0.240 | 0.635 | 1.0 | 0.747
single model | 0.976 | 0.682 | 0.808 | 1.544 | 0.807
ensemble, mean | 0.973 | 0.754 | 0.777 | 2.829 | 0.944
ensemble, min | 0.976 | 0.726 | 0.788 | 3.040 | 0.977
ensemble, MGDA | 0.979 | 0.734 | 0.800 | 3.579 | 0.949
ensemble, CAGrad | 0.976 | 0.735 | 0.774 | 2.815 | 0.924
Table 2: Max (Normalized) ground-truth score of the top 128 generated designs. For each task, the best score is green, and the second best score is blue. The “dataset” row in yellow is the normalized score of the best design in the starting offline MBO dataset. | TF Bind 8 | TF Bind 10 |  | ChEMBL
---
(Random Forest Oracle)
| Hopper
---
Controller
| Ant
---
Morphology
single model | 0.758 | 0.568 | 0.740 | 0.658 | 0.415
ensemble, mean | 0.750 | 0.685 | 0.770 | 0.655 | 0.689
ensemble, min | 0.737 | 0.688 | 0.685 | 0.646 | 0.707
ensemble, MGDA | 0.683 | 0.686 | 0.760 | 0.650 | 0.735
ensemble, CAGrad | 0.811 | 0.692 | 0.768 | 0.629 | 0.704
Table 3: 50th Percentile (Normalized) ground-truth score of the top 128 generated designs. For each task, the best score is green, and the second best score is blue. | TF Bind 8 | TF Bind 10 |  | ChEMBL
---
(Random Forest Oracle)
| Hopper
---
Controller
| Ant
---
Morphology
single model | 1.524 | 0.404 | 0.547 | 563.45 | 19.426
ensemble, mean | 1.494 | 1.112 | 0.706 | 578.11 | 251.980
ensemble, min | 1.410 | 1.138 | 0.285 | 598.55 | 281.266
ensemble, MGDA | 1.055 | 1.117 | 0.615 | 602.12 | 306.254
ensemble, CAGrad | 1.895 | 1.150 | 0.695 | 564.72 | 236.802
Table 4: Average (Unnormalized) ground-truth score of the top 128 generated
designs. For each task, the best score is green, and the second best score is
blue.

##  3 Results

We report three metrics from the top 128 designs of each algorithm: (a) max
ground truth score (Table 2), (b) 50th percentile ground truth score (Table
3), (c) average ground truth score (Table 4).

In order to report performance on the same order of magnitude across tasks, we
normalize the max ground truth score and the 50th percentile ground truth
scores using the formula

$y_{\text{normalized}}(y)=\frac{y-y_{\min}}{y_{\max}-y_{\min}}$

where $y_{\max}$ and $y_{\min}$ are the maximum and minimum objective values
in the total dataset for each task. By definition, a normalized score greater
than 1 means that we have designed an input that is better than any input in
the total dataset for the task.

We do not normalize the average ground truth scores.

Here are the results we observed:

  * •

For the max ground truth score, we find that MGDA is in the top 2 best-
performing algorithms on 4 tasks and is the best-performing algorithm on 2
tasks. CAGrad is the in the top 2 best-performing algorithms on 2 tasks.

  * •

For the 50th percentile ground truth score, MGDA is the best-performing
algorithm on 1 task. CAGrad is in the top 2 best-performing algorithms on 3
tasks and is the best-performing algorithm on 2 tasks.

  * •

Finally, for the average ground truth score, MGDA is the best-performing
algorithm on 2 tasks. CAGrad is in the top 2 best-performing algorithms on 3
tasks and is the best-performing algorithm on 2 tasks.

In general, we observe that MGDA and CAGrad performed roughly as well, if not
better, than other algorithms on the max ground truth score. However, when
considering the 50th percentile and average ground truth scores, we found that
MGDA performed much better than other algorithms on the continuous tasks and
CAGrad performed much better than other algorithms on the discrete tasks. This
suggests that that MGDA and CAGrad are more conservative and less susceptible
to being “fooled” by invalid designs.

##  4 Discussion and Future Work

###  4.1 Interpretation of Results

Based on our results, MGDA and CAGrad seem to robustify data-driven offline
MBO without compromising optimality of designs. In real-world design
scenarios, utilizing MGDA or CAGrad over mean or minimum gradient could be
well-motivated in contexts where we care about generating a diverse dataset of
good designs rather than one-off good designs. Often, computational design
projects involve repeated iteration between design generation and experimental
validation. Usually, it is more practical and efficient to experimentally
validate many proposed good designs at a time (e.g. a dataset), rather than
repeated iteration over a single design.

###  4.2 Future Work

We faced some key challenges in this work that present opportunities for
future research. First, hyperparameter selection in a purely offline manner is
difficult, and future work should explore better, more rigorous methods for
offline hyperparameter tuning.

Second, there are a lot of different approaches for using a gradient ascent
optimizer in discrete space: gradient normalization, alternating between
updates in soft-space and updates in hard-space, and more. In our
experimentation, we found that certain methods, such as gradient
normalization, improved the performance of MGDA on discrete tasks
significantly, but for consistency, we present results for a relatively simple
gradient ascent optimizer. Future work should study alternate methods for
gradient ascent in discrete space.

Finally, a key challenge we faced was with the ChEMBL task and other tasks we
tried out which don’t have an exact oracle. Using a learned oracle to evaluate
how robust our design algorithms are to out-of-distribution designs is
unreliable, because learned oracles usually suffer from the same distribution
shift problem as the proxy design models. Finding a way to reliably evaluate
design algorithms using learned oracles is an important area for future
research, because many real-world design tasks don’t have exact oracles.

Finally, although we use a diverse set of tasks, future work can study CAGrad
and MGDA on more tasks, especially with a focus on the unique characteristics
of tasks that may make MGDA more suitable than CAGrad, or vice versa.

## References

  * [Gau+11] Anna Gaulton et al.  “ChEMBL: a large-scale bioactivity database for drug discovery”  In _Nucleic Acids Research_ 40.D1, 2011, pp. D1100–D1107  DOI: 10.1093/nar/gkr777
  * [Dés12] Jean-Antoine Désidéri  “Multiple-gradient descent algorithm (MGDA) for multiobjective optimization”  In _Comptes Rendus Mathematique_ 350.5, 2012, pp. 313–318  DOI: https://doi.org/10.1016/j.crma.2012.03.014
  * [SLA12] Jasper Snoek, Hugo Larochelle and Ryan P Adams  “Practical Bayesian Optimization of Machine Learning Algorithms”  In _Advances in Neural Information Processing Systems_ 25 Curran Associates, Inc., 2012  URL: https://proceedings.neurips.cc/paper/2012/file/05311655a15b75fab86956663e1819cd-Paper.pdf
  * [TET12] Emanuel Todorov, Tom Erez and Yuval Tassa  “MuJoCo: A physics engine for model-based control”  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033  DOI: 10.1109/IROS.2012.6386109
  * [Bar+16] Luis A. Barrera et al.  “Survey of variation in human transcription factors reveals prevalent DNA binding changes”  In _Science_ 351.6280, 2016, pp. 1450–1454  DOI: 10.1126/science.aad2257
  * [Bro+16] Greg Brockman et al.  “OpenAI Gym”, 2016  eprint: arXiv:1606.01540
  * [Le+18] Daniel D. Le et al.  “Comprehensive, high-resolution binding energy landscapes reveal context dependencies of transcription factor binding”  In _Proceedings of the National Academy of Sciences_ 115.16, 2018, pp. E3702–E3711  DOI: 10.1073/pnas.1715888115
  * [Liu+21] Bo Liu et al.  “Conflict-Averse Gradient Descent for Multi-task Learning”  In _CoRR_ abs/2110.14048, 2021  arXiv: https://arxiv.org/abs/2110.14048
  * [Tra+21] Brandon Trabucco, Aviral Kumar, Xinyang Geng and Sergey Levine  “Conservative Objective Models for Effective Offline Model-Based Optimization”  In _CoRR_ abs/2107.06882, 2021  arXiv: https://arxiv.org/abs/2107.06882
  * [Tra+22] Brandon Trabucco, Xinyang Geng, Aviral Kumar and Sergey Levine  “Design-Bench: Benchmarks for Data-Driven Offline Model-Based Optimization”  In _CoRR_ abs/2202.08450, 2022  arXiv: https://arxiv.org/abs/2202.08450


# Tight Lower Bound for Average Number of Terms in Optimal Double-base Number
System

Vorapong Suppakitpaisarn11footnotemark: 1<EMAIL_ADDRESS>Graduate
School of Information Science and Technology, The University of Tokyo

###### Abstract

We show in this note that the average number of terms in the optimal double-
base number system is in $\Omega(n/\log n)$. The lower bound matches the upper
bound shown earlier by Dimitrov, Imbert, and Mishra (Math. of Comp. 2008).

###### keywords:

Double-base Number System , Asymptotic Analysis

##  1 Introduction and Notations

Given a non-negative integer $m$. A tuple $\left[k,X=\langle
x_{i}\rangle_{i=1}^{k},Y=\langle y_{i}\rangle_{i=1}^{k}\right]$ where
$k,x_{i},y_{i}$ are non-negative integers is a representation of $m$ in
double-base number system if $\sum\limits_{i=1}^{k}2^{x_{i}}3^{y_{i}}=m$. For
some $m$, there might be more than one representations of $m$ in the system.
For example, both $\left[2,\langle 0,2\rangle,\langle 1,1\rangle\right]$ and
$\left[4,\langle 0,1,2,3\rangle,\langle 0,0,0,0\rangle\right]$ are
representations of $15$.

For a representation $\left[k,X=\langle x_{i}\rangle_{i=1}^{k},Y=\langle
y_{i}\rangle_{i=1}^{k}\right]$, we call $k$ as the number of terms of the
representation. A representation with $k$ terms is a minimum representation of
$m$ if there is no representation of $m$ with $k^{\prime}$ terms such as
$k^{\prime}<k$. The number of terms of a minimum representation is denoted by
$k^{*}_{m}$ in this paper.

To speed up the calculation of scalar multiplication in elliptic curve
cryptography, many researchers devise algorithms to calculate a representation
with small number of terms for a given $m$. Those include the algorithm by
Dimitrov, Imbert, and Mishra [1]. Let $k^{\prime}_{m}$ be the number of terms
of representations obtained from the algorithm. The authors have shown that
$k^{\prime}_{m}\in O(\lg m/\lg\lg m)$ for all $m$. The result is surprising
as, in all variations of binary representations, there are infinite number of
$m$ of which all representations have $\Omega(\log m)$ terms [2]. It leads to
a smaller asymptotic complexity for calculating the scalar multiplication.

One may ask if the asymptotic complexity can be even smaller with the double-
base number system. However, in [3], we have shown that a smaller asymptotic
complexity cannot be obtained. There are infinite number of $m$ of which
$k_{m}^{*}\in\Omega(\lg m/\lg\lg m)$.

Until now, we have discuss the worst-case computation time of the scalar
multiplication. However, in literature, it is more common to analyze the
average-case computation time than to analyze the worst-case computation time
[4, 5, 6]. The average-case computation time usually depends on
$\mathcal{A}^{\prime}(n):=\sum\limits_{m=0}^{2^{n}-1}k_{m}^{\prime}/2^{n}$
when $n$ is a positive integer and $k_{m}^{\prime}$ is the number of terms in
the representation obtained from an algorithm. When the algorithm is that
proposed in [1], we have

$\mathcal{A}^{\prime}(n)=\sum\limits_{m=0}^{2^{n}-1}\frac{k_{m}^{\prime}}{2^{n}}=\sum\limits_{m=0}^{2^{n}-1}O(\frac{\lg
m}{\lg\lg m})/{2^{n}}=\sum\limits_{m=0}^{2^{n}-1}O(n/\lg n)/{2^{n}}=O(n/\lg
n).$

Let $\mathcal{A}^{*}(n):=\sum\limits_{m=0}^{2^{n}-1}k_{m}^{*}/2^{n}$. We know
that $\mathcal{A}^{*}(n)\leq\mathcal{A}^{\prime}(n)=O(n/\lg n)$.

Although it is known that the worst-case asymptotic complexity using double
base number system cannot be further improved, we may be able to further
improve the average-case asymptotic complexity. We have tried to improve our
algorithm and our analysis to have $\mathcal{A}^{\prime}(n)\in o(n/\lg n)$,
but have not been successful.

##  2 Our Result

We show that the asymptotic of $\mathcal{A}^{\prime}(n)$ cannot be further
improved, i.e. $\mathcal{A}^{*}(n)\in\Omega(n/\lg n)$. This implies that the
upper bound of $\mathcal{A}^{\prime}(n)$ is asymptotically tight, and the
algorithm in [1] is asymptotically optimal also on average case.

We will use the concept of prefix code and its properties to show the
tightness. The prefix code can be defined as follows:

######  Definition 1 (Prefix code).

Let $N$ be a positive integer and, for $i\in\\{0,\dots,N-1\\}$, let
$c_{i}\in\\{0,1\\}^{*}$. We say that $c_{0},\dots,c_{N-1}$ is a prefix code of
$0,\dots,N-1$ if $c_{i}$ is not a prefix of $c_{j}$ for $i\neq j$.

We can define an optimal prefix code as follows:

######  Definition 2 (Optimal prefix code).

Let $c_{0},\dots,c_{N-1}$ be a prefix code of $0,\dots,N-1$, let
$p_{0},\dots,p_{N-1}$ be a probability distribution on $0,\dots,N-1$, and let
$|c|$ be the length of $c\in\\{0,1\\}^{*}$

We say that $c_{0},\dots,c_{N-1}$ is an optimal prefix code under the
probability distribution $p_{0},\dots,p_{N-1}$ if there is no prefix code
$c^{\prime}_{0},\dots,c^{\prime}_{N-1}$ of $0,\dots,N-1$ such that
$\sum\limits_{i=0}^{N-1}|c^{\prime}_{i}|p_{i}<\sum\limits_{i=0}^{N-1}|c_{i}|p_{i}$.

We will use an algorithm in [7] to show our result. The algorithm takes a
probability distribution of $0,\dots,N-1$ and gives a prefix code
$c_{0},\dots,c_{N-1}$ as an output. The code obtained from the algorithm is
called Huffman code. The author of the paper has shown the following lemma.

######  Lemma 1 (Optimality of Huffman code [7]).

The output code obtained from the algorithm in [7] is an optimal prefix code
under the input probability distribution.

By Lemma 1, we can obtain the following lemma.

######  Lemma 2.

Suppose that $N$ can be written in the form of $2^{n}$ for some
$n\in\mathbb{Z}_{>0}$. Then, there is no prefix code $c_{0},\dots,c_{N-1}$
such that $\sum\limits_{i=0}^{N-1}\frac{|c_{i}|}{N}<n$.

###### Proof.

We omit the description of the algorithm of [7] in this paper. However, when
$N=2^{n}$ for some $n$, $p_{i}=1/N$ for all $i\in\\{0,\dots,N-1\\}$, and
$c_{0},\dots,c_{N-1}$ is a prefix code obtained from the algorithm, we know
that $|c_{i}|=n$ for all $i$. By the optimality of $c_{1},\dots,c_{N-1}$,
there is no prefix code $c_{1}^{\prime},\dots,c_{N-1}^{\prime}$ such that
$\sum\limits_{i=0}^{N-1}|c^{\prime}_{i}|p_{i}<\sum\limits_{i=0}^{N-1}|c_{i}|p_{i}$
or
$\sum\limits_{i=0}^{N-1}|c^{\prime}_{i}|/N<\sum\limits_{i=0}^{N-1}|c_{i}|/N=n$.
∎

Let a tuple $\left[k_{m}^{*},X^{*}=\langle
x^{*}_{i}\rangle_{i=1}^{k_{m}^{*}},Y^{*}=\langle
y^{*}_{i}\rangle_{i=1}^{k_{m}^{*}}\right]$ be a minimum representation of $m$
in double-base number system. Recall from Section 1 that
$\sum\limits_{i=1}^{k_{m}^{*}}2^{x^{*}_{i}}3^{y^{*}_{i}}=m$. We know that, for
all $i$, $x^{*}_{i},y^{*}_{i}\leq\lg m$, otherwise the term
$2^{x^{*}_{i}}3^{x^{*}_{i}}$ and the summation
$\sum\limits_{i=1}^{k}2^{x^{*}_{i}}3^{y^{*}_{i}}$ would be larger than $m$.
Therefore, we can represent $x^{*}_{i}$ and $y^{*}_{i}$ using a binary
representation length $\lg n$ when $m<2^{n}$. Also, we know that the number of
terms in the minimum representation would not be larger than $n$ when
$m<2^{n}$. We can represent $k_{m}^{*}$ using a $(\lg n)$-bit binary
representation.

Let $b[x]\in\\{0,1\\}^{\lg n}$ be a $(\lg n)$-bit binary representation of a
non-negative integer $x\in\\{0,\dots,n-1\\}$. Note that the length of $b[x]$
is fixed to $\lg n$ and the bit string $b[x]$ begins with $0$ when
$x<2^{n-1}$. We define the code
$c^{\prime\prime}_{0},\dots,c^{\prime\prime}_{2^{n}-1}$ of $0,\dots,2^{n}-1$
in the way that, for all $m$,

$c^{\prime\prime}_{m}=b[k_{m}^{*}]b[x_{1}^{*}]b[y_{1}^{*}]\dots
b[x_{k_{m}^{*}}^{*}]b[y_{k_{m}^{*}}^{*}].$

We can show the following lemma for the code.

######  Lemma 3.

The code $c^{\prime\prime}_{0},\dots,c^{\prime\prime}_{2^{n}-1}$ of
$0,\dots,2^{n}-1$ defined above is a prefix code.

###### Proof.

Assume a contradictory statement that there are $i\neq j$ such that
$c_{i}^{\prime\prime}$ is a prefix of $c_{j}^{\prime\prime}$. We must have
$b[k_{i}^{*}]=b[k_{j}^{*}]$ and $k_{i}^{*}=k_{j}^{*}$. We can then conclude
that the code $c_{i}^{\prime\prime}$ and $c_{j}^{\prime\prime}$ have the same
length. As $c_{i}^{\prime\prime}$ is a prefix of $c_{j}^{\prime\prime}$, we
have $c_{i}^{\prime\prime}=c_{j}^{\prime\prime}$. Thus, for $1\leq p\leq
k_{m}^{*}$, we have same $x_{p}^{*}$ and $y_{p}^{*}$ in the representations of
$i$ and $j$. The integer $i$ and $j$ have the same representation in double-
base number system. That is not possible when $i\neq j$. ∎

We are now ready to show our main result.

######  Theorem 4.

The average number of terms in the double-base number system, denoted by
$\mathcal{A}^{*}(n)$, is in $\Omega(n/\lg n)$.

###### Proof.

By the construction of $c_{m}^{\prime\prime}$, we have
$|c_{m}^{\prime\prime}|=(2k_{m}^{*}+1)\cdot\lg n$. Then,

$\sum\limits_{i=0}^{2^{n}-1}|c_{i}^{\prime\prime}|/2^{n}=\sum\limits_{i=0}^{2^{n}-1}(2k_{m}^{*}+1)\cdot\lg
n/2^{n}=(2\mathcal{A}^{*}(n)+1)\cdot\lg n.$

If $\mathcal{A}^{*}(n)\in o(n/\lg n)$, then

$\sum\limits_{i=0}^{2^{n}-1}|c_{i}^{\prime\prime}|/2^{n}=(2\cdot o(n/\lg
n)+1)\cdot\lg n=o(n).$

There is $n$ such that
$\sum\limits_{i=0}^{2^{n}-1}|c_{i}^{\prime\prime}|/2^{n}<n$. This contradicts
our result in Lemma 2. ∎

##  3 Concluding Remarks

In this paper, we have shown a tight lower bound for the double-base number
system in the previous section. Indeed, we can use the same argument to show
the tight lower bound for the multi-base number system. In other words, for
any constant $q$ and for any $b_{1},\dots,b_{q}$ such that $b_{i}$ and $b_{j}$
are co-prime for $i\neq j$, the average number of terms in the summation
$\sum\limits_{i=1}^{k}b_{1}^{\beta_{1}},\dots b_{q}^{\beta_{q}}=m$ would not
be smaller than $\Omega(n/\lg n)$ for $m\in\\{0,\dots,2^{n}-1\\}$.

The result also holds for the case that digit set is not $\\{0,1\\}$ as far as
the exponent is in $O(\log n)$. In other words, let the digit set be $D_{S}$,
the average number of terms in the summation
$\sum\limits_{i=1}^{k}d_{i}b_{1}^{\beta_{1}},\dots b_{q}^{\beta_{q}}=m$ for
$d_{i}\in D_{S}$ is not asymptotically smaller than that without $d_{i}$
unless we allow $\beta_{i}$ to be in $\Omega(\log n)$. This implies that our
algorithm for the multi-base number system in [3] is tight also on the average
case.

## Acknowledgement

The author would like to thank Prof. Kazuo Iwama (Kyoto University) for
suggesting us to use results from information theory in this research area.

## References

  * [1] V. Dimitrov, L. Imbert, P. Mishra, The double-base number system and its application to elliptic curve cryptography, Mathematics of Computation 77 (262) (2008) 1075–1104.
  * [2] V. Suppakitpaisarn, H. Imai, Worst case computation time for minimal joint hamming weight numeral system, in: ISITA 2014, 2014, pp. 138–142.
  * [3] D. Krenn, V. Suppakitpaisarn, S. Wagner, On the minimal hamming weight of a multi-base representation, Journal of Number Theory 208 (2020) 168–179.
  * [4] R. M. Avanzi, On multi-exponentiation in cryptography., IACR Cryptol. ePrint Arch. 2002 (2002) 154.
  * [5] E. Dahmen, K. Okeya, T. Takagi, An advanced method for joint scalar multiplications on memory constraint devices, in: ESAS 2005, Springer, 2005, pp. 189–204.
  * [6] H. Imai, V. Suppakitpaisarn, Improving width-3 joint sparse form to attain asymptotically optimal complexity on average case, IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences 98 (6) (2015) 1216–1222.
  * [7] D. A. Huffman, A method for the construction of minimum-redundancy codes, Proceedings of the IRE 40 (9) (1952) 1098–1101.