François Cayre

Marion Revolle

Isabelle Sivignon

Nicolas Le Bihan

SALZA: Practical Algorithmic Information Theory and Sharper Universal String Similarity

Keywords: Algorithmic Information Theory, Causal Inference, String Similarity, Universal Classification

Practical algorithmic information theory often relies on using off-the-shelf compressors, which both limits the range of target applications and introduces inaccuracies caused by the compressors implementation constraints and the fact that they are only able to handle multiple strings through concatenation.

We describe SALZA, a practical implementation of algorithmic information theory based on Lempel-Ziv routines, that is largely immune to the above issues. We focus on providing relationships that hold strictly in practice for strings of arbitrary length. We hope this work contributes to making clearer the links between Lempel-Ziv complexity, string similarity, and the use of compressors for computing information distances.

The capabilities of SALZA are highlighted by computing a universal semi-distance on strings, the NSD, and by performing causal inference using the PC algorithm, an application for which off-the-shelf compressors are hardly usable. SALZA was designed to take advantage of multi-core machines.

I. INTRODUCTION

I

N the discrete finite case, information theory in the probabilistic setting relies on the definition of a set A = {a k } and a discrete random variable X taking values in A equipped with a probability mass function (p.m.f.) p(a k ) = Pr {X = a k } to express the entropy H of X, the expected quantity of information that a string drawn from X conveys. Using base-2 logarithm and ignoring zero probability events, the Shannon entropy in bits is defined as:

H(X) = E p [-log p(X)] .
Defining another p.m.f. q on A, the relative entropy, or Kullback-Leibler (KL) divergence between p and q is defined as: D KL (p q) = E p log p(X) q(X) .

From there, one defines the mutual information I(X;Y ) between random variables X and Y on A as the KL-divergence between their joint p.m.f. and the product of X and Y p.m.f.'s (marginals).

And the conditional entropy H(X|Y ) is defined as H(X) -I(X;Y ).

Yet, as powerful as it is, probabilistic information theory says nothing about a particular realization x of X. This was the motivation for algorithmic information theory, independently developed by Solomonoff, Kolmogorov and Chaitin. Given a universal computer U accepting programs c of length |c|, the Kolmogorov (prefix1 ) complexity K(x) of a string x is defined as:

K(x) = min c:U(c)=x
|c|.

One similarly defines the conditional Kolmogorov complexity K(x|y) of a string x knowing another string y and eventually obtains the same kind of information calculus as in the probabilistic setting, albeit generally up to some additive constants.

With some simplification in the notations and assuming several i.i.d. x n ∼ ∏ n i=1 p(a i ), then (Theorem 7.3.1 in [START_REF] Cover | Elements of Information Theory[END_REF]):

E 1 n K(X n ) → H(X).
This relationship between Kolmogorov complexity and Shannon entropy shows the strong connection they share. Yet, the Kolmogorov complexity is uncomputable on a universal Turing computer and its use, following Solomonoff's advice, should be that of a golden standard against which practical solutions should be evaluated. In the same research line, the Normalized Information Distance (NID) has been proposed [START_REF] Bennett | Information Distance[END_REF], that serves as an ultimate distance between strings:

NID(x, y) = max {K(x|y), K(y|x)} max {K(x), K(y)} .

(1)

One of the most influential proposal for a computable complexity of an individual sequence was proposed by Lempel and Ziv [4]. This work soon after gave birth to the popular LZ family of source coding algorithms [START_REF] Ziv | A Universal Algorithm for Sequential Data Compression[END_REF], [START_REF]Compression of Individual Sequences via Variable-Rate Coding[END_REF], leading to optimal compressors for a range of source distributions. In turn, those compressors would eventually be used as practical proxies for the approximation of Kolmogorov complexity. In [START_REF] Ziv | A Measure of Relative Entropy Between Individual Sequences with Application to Universal Classification[END_REF], Ziv and Merhav used the cross-parsing

proposed in [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF] to derive an estimate of the KL-divergence between p.m.f.'s of which only one realization is known.

Using an off-the-shelf compressor C, the Normalized Compression Distance [START_REF] Li | The Similarity Metric[END_REF], [START_REF] Cilibrasi | Clustering by Compression[END_REF]

The expression of the NCD implicitly assumes |C(xy)| = |C(yx)|, which is an approximation that depends on the particular compressor used in practice. In effect, a compressor has to deal with a few constraints, and being aware of them may help selecting the right compressor for the data at hand [START_REF] Cebrián | Common Pitfalls Using the Normalized Compression Distance: What to Watch Out for in a Compressor[END_REF].

But another limitation in using compressors stems from the way they enforce the concatenationbased approximation to compute an estimate of conditional information C(y|x) ≈ C(xy) -C(x):

while this is both practical and theoretically satisfying when two strings are involved in the computation, approximations can only get worse when more strings are involved -for example, performing causal inference would require to compute something like C(y|x 1 , . . . , x n ) efficiently and reliably, for a number of varying y, x 1 , . . . , x n . It is not clear how one could guarantee to do that in practice. At the very least, this would certainly imply managing concatenations in a somewhat cumbersome and inefficient fashion.

We argue that using off-the-shelf compressors hinders what can be actually implemented to leverage the power of algorithmic information estimates for applications where conditioning must be computed with much greater flexibility. Although using algorithmic information has been justified for causal inference [START_REF] Janzing | Causal Inference Using the Algorithmic Markov Condition[END_REF], practical realizations on this side have been scarce.

Hence, we believe that an implementation of algorithmic information theory [START_REF] Steudel | Causal Markov Condition for Submodular Information Measures[END_REF] is both possible and needed, ensuring that (i) the required flexibility in conditioning is easily achieved and (ii)

the computed quantities verify usual information-theoretic relationships for strings of arbitrary length.

A. Notations

The empty string, set or multiset is denoted by / 0. The set of non-empty finite strings on alphabet A is denoted by A + and A = A + ∪ / 0. The length of alphabets, strings, sets and multisets are denoted with |.|. The symbol + is also used to denote addition of multisets. The first k strings in a set are denoted with x ≤k and by convention x ≤0 = / 0. The concatenation of two strings x and y is denoted with xy. Logarithms are in base 2. Terms related to implementation or simulations are typeset using monospaced font.

B. Lempel-Ziv symbols, between theory and practice

We shall make extensive use of Lempel-Ziv-type factorization [START_REF] Ziv | A Universal Algorithm for Sequential Data Compression[END_REF] and cross-factorization [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF].

Hence, we start by making clear some basic implementation choices.

We first remark that popular compressors from the Lempel-Ziv family (most notably gzip [START_REF] Deutsch | DEFLATE Compressed Data Format Specification, version 1.3[END_REF], [START_REF]GZIP File Format Specification[END_REF] and lzma [START_REF] Pavlov | 7z Format[END_REF]) nowadays are based on LZ77 [START_REF] Ziv | A Universal Algorithm for Sequential Data Compression[END_REF] instead of LZ78 [START_REF]Compression of Individual Sequences via Variable-Rate Coding[END_REF]. The reasons for this are that computer memory has become cheaper and it is easier to manage a sliding window (from 32KiB in [START_REF] Deutsch | DEFLATE Compressed Data Format Specification, version 1.3[END_REF] to 1MiB by default and maximum 4GiB in [START_REF] Pavlov | 7z Format[END_REF]) than an evolving dictionary of substrings not seen so far.

We will also use algorithms similar to those based on a sliding window, except that its size will be in practice limited by the host memory (pointers are 32-bit wide in SALZA).

Further, since the inception of the Lempel-Ziv source coding techniques [START_REF] Ziv | A Universal Algorithm for Sequential Data Compression[END_REF], [START_REF]Compression of Individual Sequences via Variable-Rate Coding[END_REF], [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF], where a symbol is composed of a reference to the next substring and the next following character, most expositions now use only references or characters -with characters being either described as length-1 references (most notably in the computer science literature, see [START_REF] Crochemore | A Simple Algorithm for Computing the Lempel-Ziv Factorization[END_REF] and related works)

or transmitted as such (as in off-the-shelf compressors, e.g. as in [START_REF] Deutsch | DEFLATE Compressed Data Format Specification, version 1.3[END_REF]). The same discrepancy also appears between the original exposition of [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF] and some of its later uses [START_REF] Ziv | A Measure of Relative Entropy Between Individual Sequences with Application to Universal Classification[END_REF], [START_REF] Hoobin | Relative Lempel-Ziv Factorization for Efficient Storage and Retrieval of Web Collections[END_REF].

In what follows, this distinction will be irrelevant, to the exception of Sec. II-D.

II. SALZA AS A MEASURE OF STRING SIMILARITY

The first step in this work is to explain how a generic factorization routine relates to string similarity and to previous works [START_REF] Lempel | On the Complexity of Finite Sequences[END_REF] on computable complexity. Also, we shall discuss why compressors may not always be the best tools for similarity assessment.

A. SALZA factorizations

SALZA is built around a routine for sequential factorization (parsing) of a string y when given as prior knowledge a set of strings x 1 , . . . , x n . Following the Lempel-Ziv basic idea, SALZA will look for the next longest substring that matches the beginning of the remaining part of y to be factorized. Each string is associated to its own search structure (see Appendix A-B for details) so that the next longest substring cannot overlap between strings.

Two kinds of factorization were implemented, depending on whether the part of y that was already factorized (the prefix of y) is included in the search space or not.

The first kind of factorization will be denoted with y|x 1 , . . . , x n : at each factorization step, the next longest substring is to be found either in one of the x 1 , . . . , x n , or in the prefix of y.

The second kind of factorization will be denoted with y| + x 1 , . . . , x n : at each factorization step, the next longest substring is to be found only in one of the x 1 , . . . , x n . This can be seen as a generalization of the procedure described in [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF] allowing for several "databases".

When the kind of factorization is left unspecified, it will be denoted with y x 1 , . . . , x n .

Hence, at each factorization step, we look for the longest substring in all prior knowledge strings x 1 , . . . , x n (also including the prefix of y for the factorization |) and we keep the longest of them as the next factor. Definition 2.1 (SALZA symbols and lengths): Given strings y, x 1 , . . . , x n in A , SALZA will factorize y into m symbols:

y x 1 , . . . , x n = (s 1 , l 1 , z 1 ) . . . (s k , l k , z k ) . . . (s m , l m , z m ),
where s k ∈ {y, x 1 , . . . , x n } is (a pointer to) the string in which the next longest substring was found, 1 ≤ l k ≤ |s k | is the length of the next longest substring, and 1 ≤ z k ≤ |s k | is a pointer in string s k where the next longest substring starts.

All SALZA symbol lengths are collected in the multiset L y x 1 ,...,

x n = {l k }, 1 ≤ k ≤ m.
Remark 1 (Invariance of y x 1 , . . . , x n ): y x 1 , . . . , x n is invariant with respect to the order of the prior knowledge strings x 1 , . . . , x n because the sequential factorization of y only evolves after all strings have been searched for the longest substring. be the first i lengths in the factorization y x. To each i, we associate v(i) the smallest integer

such that i ∑ k=1 l k y x ≤ v(i) ∑ k=1 l k y x,z .
We shall show that ∀i , v(i) ≤ i (which proves the claim by setting i = |L y x |). The proof is inductive:

For the base case (i = 1), v(1) = 1 because z could only bring a longer reference (so that l 1 y x ≤ l 1 y x,z ). Hence, v(1) = 1; For the inductive step, we shall show that:

v(i) ≤ i =⇒ v(i + 1) ≤ i + 1.
By definition of v(i), only two cases occur:

• (Case 1) i ∑ k=1 l k y x = v(i) ∑ k=1 l k y x,z .
In this case, the two factorizations are in sync and the exact next same longest substring is to be searched with or without knowledge of z. This next substring may only be longer in z so v(i + 1) = v(i) + 1, and (by the premise):

v(i + 1) = v(i) + 1 ≤ i + 1; • (Case 2) i ∑ k=1 l k y x < v(i)
∑ k=1 l k y x,z . In this case, the factorization y x, z is in advance over y x, and only two cases occur depending whether step i + 1 in y x will overtake step v(i) of y x, z or not:

(Case 2.1) i+1 ∑ k=1 l k y x ≤ v(i) ∑ k=1 l k y x,z . Then v(i + 1) = v(i)
and (also by the premise):

v(i + 1) = v(i) ≤ i < i + 1; (Case 2.2) i+1 ∑ k=1 l k y x > v(i)
∑ k=1 l k y x,z . In the worst case, the factorization step v(i) + 1 will see y x back in sync. Otherwise, y x, z keeps ahead. Hence,

l v(i)+1 y x,z ≥ i+1 ∑ n l n y x - v(i)
∑ n l n y x,z and v(i + 1) = v(i) + 1, which finally shows the claim similarly to Case 1 above. Definition 2.2 (Product of SALZA factorizations): Let y 1 and y 2 two strings being factorized, each with respective prior knowledge strings x 1,1 , . . . , x 1,m 1 and x 2,1 , . . . , x 2,m 2 . Let also:

y 1 x 1,1 , . . . , x 1,n 1 = (s 1,1 , l 1,1 , z 1,1 ) . . . (s 1,m 1 , l 1,m 1 , z 1,m 1
), and

y 2 x 2,1 , . . . , x 2,n 2 = (s 2,1 , l 2,1 , z 2,1 ) . . . (s 2,m 2 , l 2,m 2 , z 2,m 2 ).
We define their factorization product as the concatenation of their SALZA symbols:

y 1 x 1,1 , . . . , x 1,n 1 × y 2 x 2,1 , . . . , x 2,n 2 = (s 1,1 , l 1,1 , z 1,1 ) . . . (s 1,m 1 , l 1,m 1 , z 1,m 1 ) (s 2,1 , l 2,1 , z 2,1 ) . . . (s 2,m 2 , l 2,m 2 , z 2,m 2 ).
Definition 2.3 (SALZA joint factorization): The joint factorization of x 1 , . . . , x n ∈ A is defined as the following product of factorizations:

x 1 • . . . • x n = n ∏ i=1 x i x ≤i-1 .
Hence, x| / 0 may therefore be used to denote the usual LZ77 factorization of x and |x| = |L x| + / 0 |. Remark 3 (On concatenation and symmetry): Because the factorization is a sequential operation:

1) xy| / 0 = x| / 0 × y|x and therefore L xy| / 0 = L x| / 0 + L y|x ; 2) xy| + / 0 = x| + / 0 × y| + / 0 and therefore L xy| + / 0 = L x| + / 0 + L y| + / 0 ; 3) x • y = y • x in general.

B. Relative string similarity

If the multiset L y x 1 ,...,x n is seen as an integer composition of |y|, we can associate it to the integer partition of |y| obtained by sorting the lengths in L y x 1 ,...,x n in decreasing order. In the following, we assume L y x 1 ,...,x n is sorted in decreasing order and it is therefore a partition of |y|. When the context makes it clear, we either write S(L) or S(y x 1 , . . . , x n ) instead of S(L y x 1 ,...,x n ).

{1, 1, 1, 1, 1, 1, 1} {2, 1, 1, 1, 1, 1} {2, 2, 1, 1, 1} {2, 2, 2, 1} {3, 1, 1, 1, 1} {3, 2, 1, 1} {3, 2, 2} {3, 3, 1} {4, 1, 1, 1} {4, 2, 1} {4, 3} {5, 1, 1} {5, 2} {6, 1} {7}
Fig. 1: Binary tree of the integer partitions of 7. To each level corresponds the same number of summands. Whenever possible, growing downwards is done either by replacing two 1's with one 2 (left child), or by removing a 1 and increasing the smallest summand greater than 1 if it has multiplicity 1 (right child).

Examples of score functions S include counting the elements of the multiset and the length of the compressed string y according to the factors y x 1 , . . . , x n . We are interested in finding a score function that is as much precise as possible, even for short strings, so as to cover the full range of possible applications.

Even though no closed-form expression exists for the number of partitions of an integer 2 , they can be constructed along a binary tree [START_REF] Fenner | A Binary Tree Representation and Related Algorithms for Generating Integer Partitions[END_REF], see Fig. 1, that can be used to induce an order of interest (Proposition 2 and Table I). Such a binary tree is called the partition tree of the corresponding integer.

We slightly adapt notations from [START_REF] Fenner | A Binary Tree Representation and Related Algorithms for Generating Integer Partitions[END_REF] to our setting. A partition of a positive integer |y| is a m-tuple L = {l 1 , l 2 , . . . , l m } of positive integers (the summands of the partition) such that:

l 1 + l 2 + • • • + l m = |y|, l 1 ≥ l 2 ≥ • • • ≥ l m ≥ 1.
When it exists, the number of summands strictly greater than 1 is denoted by h, 1 ≤ h ≤ m, i.e.:

L = {l 1 , . . . , l h > 1, l h+1 = 1, . . . , l m }.
2 The partition function P(n) gives the number of partitions of an integer n. In 1918, Hardy and Ramanujan proposed the following asymptotic approximation:

P(n) ≈ 1 4n √ 3 exp π 2n 3 as n → ∞. {1, 1, 1, 1, 1, 1, 1} {2, 1, 1, 1, 1, 1} {2, 2, 1, 1, 1} {3, 1, 1, 1, 1} {2, 2, 2, 1} {3, 2, 1, 1} {4, 1, 1, 1} {3, 2, 2} {3, 3, 1} {4, 2, 1} {5, 1, 1} {4, 3} {5, 2} {6, 1} {7}
TABLE I: List of integer partitions of 7 obtained when traversing the tree of Fig. 1 in breadth-first order.

Given a partition L, one can define two operations acting on L, namely A and B, such as:

• A: Add one summand of value 2 and remove two of value 1;

• B: If the smallest summand strictly greater than 1 has multiplicity 1 (l h-1 > l h ), increase it by 1 and remove a summand of value 1.

When operation A or B can3 be applied to a partition L, the resulting partition is denoted by LA or LB. The binary tree enumerating all integer partitions of n is constructed starting from the root ⊥ = {l 1 = 1, . . . , l m } by applying A (resp. B) to get the left (resp. right) child of a partition whenever possible.

In the sequel, necessary conditions are assumed to hold true when the result of an operation on a partition is stated.

In [START_REF] Fenner | A Binary Tree Representation and Related Algorithms for Generating Integer Partitions[END_REF], Corollary 13 states that if the nodes of the partition tree are listed in top-down left-to-right level (breadth-first) order, the resulting list is in decreasing order of the number of summands and in lexicographic order for partitions having the same number of summands (see Table I).

1) ∀ 0 ≤ i ≤ |y|, S(⊥B i ) ≥ S(⊥A i+1 ) ; 2) ∀ 0 ≤ i ≤ |y|, S(LAB i ) ≥ S(LBA i ).
Then S complies with the breadth-first order.

Proof: Consider the binary alphabet B = {A, B}. We label any node x of the partition tree by the word on B composed of the sequence of operations performed from the root ⊥ to x and denote it by label(x). By convention, label(⊥) = / 0. Then, following the notations introduced before, the integer partition of |y| represented by the node x is denoted by r(x) = ⊥label(x).

Consider the ordered sequence {x k } k of nodes given by a breadth-first traversal of the integer partition tree of |y|. If A < B, then the nodes are ordered according to the radix order 4 (also called shortlex order) of their labels (see [START_REF] Knuth | The Art of Computer Programming[END_REF] 2) and given two consecutive nodes x k and

x k+1 , we have to show that S(r(x k )) ≥ S(r(x k+1 )). Two cases occur:

1) if x k and x k+1 are not on the same level of the tree, then there exist 0 ≤ i ≤ |y| such that label(x k ) = B i and label(x k+1 ) = A i+1 . We have S(r(x k )) = S(⊥B i ) and S(r(x k+1 )) = S(⊥A i+1 ), and the result follows from Condition (1); 2) otherwise x k and x k+1 are on the same tree level and there exists 0 ≤ i ≤ |y| and a word M on B such that label(x k ) = MAB i and label(x k+1 ) = MBA i . S(r(x k )) = S(⊥MAB i ) and S(r(x k+1 )) = S(⊥MBA i ), and the result follows from Condition (2) with L = ⊥M.

C. Order compliant score functions

One immediately verifies that S 1 (L) = |L| -1 complies with the breadth-first order. However, it is unable to make a difference between partitions located on the same tree level.

Definition 2.5 (SALZA score function S p ):

The SALZA score function of a partition L = L y x 1 ,...,x n = {l 1 , . . . , l m } is defined as:

S p (L) =      0 if L = / 0 |L| -1 |y|+1 -∑ p k=1 (|y| + 1) -k l k otherwise, (3)
where p ≥ 1 is a parameter allowing to tune the precision and

l k = 0 if k > m.
In other words, the fractional part of S p is built by writing the l k 's in radix (|y| + 1) (if not for the corrective term 1/(|y| + 1), see below). And S p may be seen as an arbitrary-precision version of S 1 . The value of p is set according to Appendix A-A.

Proposition 3 (SALZA score is compliant with breadth-first order):

1) ∀ 0 ≤ i < |y|, S p (⊥B i ) ≥ S p (⊥A i+1 ) ; 2) ∀ 0 ≤ i < |y|, S p (LAB i ) ≥ S p (LBA i ).
Proof: Let L = {l 1 , . . . , l h > 1, l h+1 = 1, . . . , l m } a partition of |y|. In the sequel, we shall repeatedly apply operations A and/or B to a partition, leading to different values of h (and l h )

for the resulting partitions. However, we want to express what happens to L after operations A and/or B have been applied, so we require that h and l h refer to that of L -and the actual h and l h of the resulting partitions can be easily deduced from the following expressions:

⊥A i = l 1 = 2, . . . , l i = 2, l i+1 = 1, . . . , l |y|-i , ⊥B i = l 1 = 1 + i, l 2 = 1, . . . , l |y|-i , LAB i = {l 1 , . . . , l h , l h+1 = 2 + i, l h+2 = 1, . . . , l m-i-1 } , LBA i = {l 1 , . . . , l h + 1, l h+1 = 2, . . . , l h+i = 2, l h+i+1 = 1, . . . , l m-i-1 } . Let ∆ 1 (i) = S p (⊥B i ) -S p (⊥A i+1
) and p = |y|i:

∆ 1 (i) = 1 - i -1 |y| + 1 + i+1 ∑ k=2 (|y| + 1) -k -(|y| + 1) -(|y|-i) .
The expression of ∆ 1 (i) is written such that decreasing p removes terms starting from the right. For instance, when p < |y| < i, the last term (|y| + 1) |y|-i disappears. Lowest values of ∆ 1 (i) are obtained for p = |y|i or p = 1 (in which case all positive terms of the sum disappear.)

For p = |y| -i, since 0 ≤ i < |y|, (|y| + 1) -(|y|-i) ≤ (|y| + 1) -1 , so that: i -1 |y| + 1 + 1 (|y| + 1) |y|-i ≤ i |y| + 1 < 1.
For p = 1, using the fact that i < |y|, we get i -1 < |y| + 1 and 1 -(i -1)/(|y| + 1) > 0. Hence:

∀ |y| ≥ 0, ∀ p ≥ 1, ∀ 0 ≤ i < |y|, ∆ 1 (i) ≥ 0. Let ∆ 2 (i) = S p (LAB i ) -S p (LBA i ) and p = m -i: ∆ 2 (i) = (|y| + 1) -h -i(|y| + 1) -(h+1) + h+i ∑ k=h+2 (|y| + 1) -k .
One has that: 1

(|y| + 1) h - i (|y| + 1) h+1 = |y| + 1 -i (|y| + 1) h+1 . Since i < |y|: ∀ |y| ≥ 0, ∀ p ≥ 1, ∀ 0 ≤ i < m ≤ |y|, ∆ 2 (i) ≥ 0.
Definition 2.5 provides an insight into the way SALZA implements string similarity: strings y and x 1 , . . . , x n will be all the more similar when (i) the associated factorization is shorter and

(ii) when it contains more longer factors.

When the set of prior knowledge strings is empty, two cases arise depending on the kind of factorization used:

• Case y| + / 0: no substring can ever be found, so we end up enumerating the characters of y and S p (y| + / 0) will merely reflect the length of y;

• Case y| / 0: substrings can only be found in the part of y that has already been parsed, and S p (y| / 0) can therefore be seen as a measure of self-similarity. Compared to [START_REF] Lempel | On the Complexity of Finite Sequences[END_REF] and [START_REF] Ziv | A Measure of Relative Entropy Between Individual Sequences with Application to Universal Classification[END_REF] where one counts the number of factors in a factorization (thus amounting to selecting a level on the partition tree with S 1 ), SALZA will provide a more finegrained measure because it is able to make a difference between partitions (and associated factorizations) located on the same tree level. Further, this difference is associated to a rather intuitive measure of similarity based on the length of the factors -which is backed by empirical evidence elsewhere [START_REF] Cancedda | Word-Sequence Kernels[END_REF], that longer matches help in building more discriminative features for string similarity (the Noisy Stemming Hypothesis in [START_REF] Cancedda | Word-Sequence Kernels[END_REF]).

Proposition 4 (S p (y x 1 , . . . , x n ) is positive and bounded):

0 ≤ S p (y x 1 , . . . , x n ) < |y|. (4)
Proof: When y = / 0, L y x 1 ,...,x n = / 0 and S p (y x1, . . . , x n ) = 0. When |y| > 0, Proposition 3 states that the minimum is reached when |L y x 1 ,...,x n | = 1 (y is a substring of any of the x 1 , . . . , x n ), and the maximum is reached when |L y x 1 ,...,x n | = |y| (no substring of any of the x 1 , . . . , x n matches any substring of y -this is the root of the partition tree).

Hence, for |y| > 1, the minimum is:

1 - 1 |y| + 1 - |y| |y| + 1 = 0 ≤ S p (y x 1 , . . . , x n ),
and the maximum is:

S p (y x 1 , . . . , x n ) ≤ |y| - 1 |y| + 1 - p ∑ k=1 (|y| + 1) -k < |y|.
Remark 4: For strings y of size 1, both S 1 and S p equal zero.

D. Comparing with compression

Most practical compressors may be seen as two-stage pieces of software engineering. The first stage will extract some regularity information from the string (LZ coding, Burrows-Wheeler transform (BWT) [START_REF] Burrows | A Block-Sorting Lossless Data Compression Algorithm[END_REF]), while the second stage will use entropy coding in order to converge faster towards entropy (Huffman [START_REF] Huffman | A Method for the Construction of Minimum-Redundancy Codes[END_REF] or arithmetic coding, possibly intertwined with Move-to-Front strategy). The intricacies between these two stages make modeling of the performances of real compressors a tedious task (most notably, LZ77-based compressors will further quantize reference lengths and distances prior to encoding them). In the literature, tractable expressions are obtained by implementing the second stage with universal coding of integers [START_REF] Elias | Universal Codeword Sets and Representations of the Integers[END_REF].

Here, we use an estimation of the compressed string length (à la LZ77) for a score function both for the length and offset of a reference, we shall spend log l k + 2 log(1 + log l k ) + log d + 2 log(1 + log d) bits for a reference 5 .

In this setting, the compressed string length based on the L = L y x 1 ,...,x n partition of |y| reads:

S c (L) = h ∑ k=1 log c(1 + log l k ) 2 l k + |L| ∑ k=h+1 log |A|, where c = log d(1 + log d) 2 . Let now ∆ c,1 (i) = S c (LB i ) -S c (LA i+1 ): ∆ c,1 (i) = log(c(l h + i)(1 + log(l h + i)) 2 ) + h+i+1 ∑ k=h+1 log |A| + m-i-1 ∑ k=h+i+2 log |A| + log |A| -log(cl h (1 + log l h ) 2 ) - h+i+1 ∑ k=h+1 log 8c - m-i-1 ∑ k=h+i+2 log |A| ∆ c,1 (i) = log |A| (l h + i)(1 + log(l h + i)) 2 l h (1 + log l h ) 2 + (i + 1) log |A| 8c . Let ∆ c,2 (i) = S c (LAB i ) -S c (LBA i ). Similarly: ∆ c,2 (i) = log l h (1 + log l h ) 2 (l h + 1)(1 + log(l h + 1)) 2 + log (2 + i)(1 + log(2 + i)) 2 8 + (i -1) log |A| 8c .
Neither ∆ c,1 nor ∆ c,2 is positive for all values of i ≥ 0 and |y| > 1, only when i is large enough -meaning that compressed string length will be a meaningful score only when the differences between the two factorizations will be large enough (the partitions corresponding to the two factorizations are on two sufficiently distant levels of the partition tree). Even using a sliding window of size w ≤ d would only marginally help (replacing d by w in c, although the actual coding scheme would still have to be defined).

The above should be mitigated given the following two observations:

1) A compressor is not designed to specifically comply with the two conditions of Proposition 2 -but even if this comparison may seem not so fair, it nevertheless suggests that compressors are likely to be less precise tools for assessing similarity in general;

2) Real compressors use several tricks that we are unable to model but that allow to converge faster towards entropy (they are, however, unable to handle conditioning with as much flexibility as SALZA).

III. SALZA AS AN IMPLEMENTATION OF INFORMATION THEORY

Our notion of string similarity is not defined in the asymptotic regime, which is how performances of compressors are evaluated. Now, even if information theory is primarily concerned with very long strings, we are nevertheless in position to list various information-theoretic properties for SALZA that hold strictly in practice for arbitrary string lengths.

The biggest departure from the probabilistic setting is that, because factorizations are sequential processes, we will be limited to asymmetric relationships. For a number of relevant applications, however, this is not an issue (see Sec. IV and V).

We shall use the subscript f in the notation S f to denote either S 1 or S p , since both were shown to comply with our notion of string similarity. When it is omitted, we assume S f = S p .

A. Basic tools and conditioning Definition 3.1 (SALZA joint and self measures): Given strings x 1 , . . . , x n ∈ A * and following Def. 2.3, the SALZA joint measure is defined as:

S f (x 1 • . . . • x n ) = n ∑ i=1 S f (x i x ≤i-1 ).
Remark 5 (On joint and self measures):

1) By Def. 2.4, the order of the x 1 , . . . , x n does matter;

2) When applied to a single string, the notation for joint measure may gracefully degrade into that of the LZ77-based computation of a measure on self (factorization |). Further, in the case of a single string, the context will help determine which kind of factorization is meant.

Definition 3.2 (SALZA conditional mutual measure, asymmetric version): Given strings x, y, z ∈ A * , the SALZA conditional mutual measure of x and y given z is defined as:

I f (x : y z) = S f (z • x) + S f (z • y) -S f (z • x • y) -S f (z). (5)
If I f (x : y z) = 0, x and y are said to be independent given z.

B. Properties

Proposition 5 (S f (y x 1 , . . . , x n ) is non-increasing by conditioning):

Let y ∈ A and {x 1 , . . . , x n+1 } a set of n + 1 strings in A .

∀ n ≥ 0, S f (y x ≤n+1 ) ≤ S f (y x ≤n ).

Proof: Proposition 1 only allows three cases:

• (Case 1) |L y x ≤n+1 | < |L y x ≤n |, in which case the integer partition associated to L y x ≤n+1 is located on a strictly lower level of the partition tree of |y| than that associated to L y x ≤n and Proposition 3 ensures that S f (y x ≤n+1 ) < S f (y x ≤n ) ;

• |L y x ≤n+1 | = |L y x ≤n | and the integer partition associated to L y x ≤n+1 is -(Case 2.1) the same as that associated to L y x ≤n (string x n+1 did not bring new knowledge on y) and S f (y x ≤n+1 ) = S f (y x ≤n ) ; -(Case 2.2) in strictly increasing lexicographic order than that associated to L y x ≤n (because at least one longer substring was found during the factorization y x ≤n+1 compared to y x ≤n ) and Proposition 3 ensures that S f (y x ≤n+1 ) ≤ S f (y x ≤n ).

Setting n = 0 shows that the maximum is reached for S f (y / 0) = S f (y).

Proposition 6 (S f is subadditive): Given strings x, y ∈ A ,

S f (x • y) ≤S f (x) + S f (y); S f (y • x) ≤S f (x) + S f (y).
Proof: By Def. 3.1:

S f (x • y) = S f (x) + S f (y x) ≤ S f (x) + S f (y) (by Proposition 5).
The case S f (y • x) is similar.

Proposition 7 (Kolmogorov's formula6 for SALZA): Given strings x, y, z ∈ A ,

I f (z • x : y) = I f (z : y) + I f (x : y z).
Proof: By Eq. 5 on all three terms:

I f (z • x : y) = S f (z • x) + S f (y) -S f (z • x • y) I f (z : y) = S f (z) + S f (y) -S f (z • y) I f (x : y z) = S f (z • x) + S f (z • y) -S f (z • x • y) -S f (z).
Proposition 8 (Fast computation of I f (x : y z)):

I f (x : y z) = S f (y z) -S f (y x, z). (7)
Proof: By Eq. 5 and Def. 3.1:

I f (x : y z) = S f (z • x) + S f (z • y) -S f (z • x • y) -S f (z) = S f (z) + S f (x z) + S f (z) + S f (y z) -S f (z) + S f (x z) + S f (y z, x) -S f (z) I f (x : y z) = S f (y z) -S f (y x, z).
We have used Remark 1 relative to the invariance with respect to permutation of the prior knowledge strings (S f (y z, x) = S f (y x, z)).

Proposition 9 (SALZA conditional mutual measure is positive): Given strings x, y, z ∈ A , 0 ≤ I f (x : y z).

Proof: By Eq. 7 and Proposition 5.

C. SALZA as an information measure Proposition 10 (Chain rule for SALZA conditional mutual measure, asymmetric version):

Given strings x, y, z,t ∈ A ,

I f (x : y • z t) = I f (x : y t) + I f (x : z t, y). (9)
Proof: By Eq. 5:

I f (x : y • z t) = S f (t • x) + S f (t • y • z) -S f (t • x • y • z) -S f (t) = S f (t • x) -S f (t) + S f (y t) -S f (x t)
-S f (y x,t) + S f (z t, y) -S f (z x, y,t).

Also by Eq. 5:

I f (x : y t) = S f (t • x) -S f (t) + S f (t • y • z) -S f (t • x • y • z) = S f (t • x) -S f (t) + S f (y t) -S f (x t)
-S f (y x,t).

By Eq. 7:

I f (x : z t, y) = S f (z t, y) -S f (z x, y,t).
Proposition 11 (Data processing inequality for SALZA, asymmetric version): Given strings

x, y, z ∈ A :

S f (y z) = 0 =⇒ I f (x : y z) = 0 =⇒ I f (x : y) ≤ I f (z : y).
Proof: The first implication is clear by the proof of Proposition 4 since S f (y z) = 0 either implies that y = / 0 or y is a substring of z. By Eq. 7, I f (x : y z) = 0. Also by Eq. 7, I f (x : y z) = 0 ⇐⇒ S f (y z) = S f (y x, z). The same equation further states that:

I f (x : y) = S f (y) -S f (y x), I f (z : y) = S f (y) -S f (y z).
Replacing S f (y z) with S f (y x, z) in the latter, one has that I f (z : y) = S f (y) -S f (y x, z). By Proposition 5, -S f (y x) ≤ -S f (y x, z) and the second implication holds.

We now show that S f satisfies the same requirements as in Sec. 7.1 of [START_REF] Steudel | Causal Markov Condition for Submodular Information Measures[END_REF] to qualify as a submodular information measure. Such an information measure is defined on the lattice of strings, using the lexicographical order as ≤. We shall rephrase here the approach of [START_REF] Steudel | Causal Markov Condition for Submodular Information Measures[END_REF] to approximate submodularity using joint factorizations.

Let s,t ∈ A and z the concatenation of the unique substrings shared by s and t. Let x (resp. y) the concatenation of the substrings in s (resp. t) that are not in z such that |s| = |zx| (resp. |t| = |zy|). Then, using joint factorizations, the lattice operators are computed like s ∨ t = z • x • y and s∧t = z. Now, the submodularity property (S f (s)+S f (t) ≥ S f (s∧t)+S f (s∨t)) translates into positivity of the conditional mutual information I f (x : y z). Note that this is only approximate submodularity, since we do not ensure commutativity of the lattice operators (S f (x • y) = S f (y • x) in general, see Fig. 2).

Proposition 12 (S f is an information measure in the sense of [START_REF] Steudel | Causal Markov Condition for Submodular Information Measures[END_REF]): Given strings x, y, z ∈ A : 1) Normalization: S f ( / 0) = 0; 2) Monotonicity (over prefixes [START_REF] Li | An Introduction to Kolmogorov Complexity and Its Applications[END_REF], p. 197): x ≤ y =⇒ S f (x) ≤ S f (y); L y / 0 ) and let ∆(x, y) = S p (L y / 0 ) -S p (L x / 0 ):

3) Approximate submodularity: S f (z • x) + S f (z • y) ≥ S f (z • x • y) + S f (z).
∆(x, y) = |L y / 0 | -|L x / 0 | + |z| (|x| + |z| + 1)(|x| + 1) + p ∑ k=1 (|x| + 1) -k l (x) k -(|x| + |z| + 1) -k l (y) k ≥ |L y / 0 | -|L x / 0 | + |z| (|x| + |z| + 1) 2 + p ∑ k=1 (|x| + 1) -k l (x) k -(|x| + |z| + 1) -k l (y) k
The first difference above is positive, as well as those in the sum for 1 The asymmetry of S p is tedious to derive analytically. We use four datasets, namely markov (see infra), and mammals, languages and music from [START_REF] Li | The Similarity Metric[END_REF], [START_REF] Cilibrasi | Clustering by Compression[END_REF] to compute departure from symmetry like:

≤ p ≤ |L x / 0 | (if
2|S p (x • y) -S p (y • x)| S p (x • y) + S p (y • x) .
The heatmap of Fig 2 represents departure from symmetry for all pairs of strings in the four datasets we consider. While the mean departure is 0.24%, only some strings (a handful of MIDI music scores) will cause more severe deviation.

From now on, we start specifying explicitly the kind of factorization suited to dedicated applications.

IV. NORMALIZED SALZA SEMI-SISTANCE

As a first application of SALZA, we now devise a universal semi-distance on strings, the Normalized SALZA semi-Distance (NSD). It is universal in the sense that it does not use modeling of the data. We shall compare it to the NCD [START_REF] Li | The Similarity Metric[END_REF], [START_REF] Cilibrasi | Clustering by Compression[END_REF] for classification purposes, since we seek to provide objective measures of performance (which can be hardly derived from the phylogenetic trees produced for clustering applications).

In doing so, we shall describe a straightforward and reproducible way of producing strings from continuously-valued attributes -hence encompassing both discrete and continuous datasets for clustering/classification applications. First, we review the NCD from the point of view of the conditioning.

A. Review of NCD implementation

An information distance between strings x and y is the length of the shortest program for a universal computer to transform x into y and y into x (see [START_REF] Li | An Introduction to Kolmogorov Complexity and Its Applications[END_REF], p. 642). This ultimately gave birth to the NID (Eq. 1). Using an actual compressor C as a computable proxy for the uncomputable K, one defines the NCD (Eq. 2) as a similarity metric on strings. Note that while the NID is actually a distance (up to negligible errors -see [START_REF] Li | An Introduction to Kolmogorov Complexity and Its Applications[END_REF], Theorem 8.4.1, p. 662), using an actual compressor may cause deviation from this ideal behaviour.

Let us now consider how the NCD works and focus on the numerator when a LZ77-based compressor is used (e.g., gzip or lzma). When computing C(xy), the compression is done sequentially. When the sliding window is moving over the boundary of x and y, it is essentially the same as if we were using the factorization y|x. But when the sliding window is entirely over y, no regularity information from x can ever be taken into account anymore. The same conclusions may be drawn when considering a block-BWT compressor (e.g., bzip2). SALZA, because it does not work by blocks or features a sliding window limited by the host memory, is largely immune to this issue. And it enforces a clean use of the x| + y conditioning without those side-effects.

Yet, the NCD (that assesses string similarity) is not the only approach to universal clustering/classification. Ziv and Merhav [START_REF] Ziv | A Measure of Relative Entropy Between Individual Sequences with Application to Universal Classification[END_REF] have proposed an approach that relies on estimating the KL-divergence (which is akin to assessing dissimilarity between strings) by using x| + y factorizations [START_REF] Wyner | Fixed Data Base Version of the Lempel-Ziv Compression Algorithm[END_REF]. We shall see that assessing either similarity or dissimilarity can lead to two semi-distances.

B. Definition of the NSD Definition 4.1 (Normalized SALZA Semi-Distance): Given two strings x, y defined on A such that |x|, |y| > 1, the normalized SALZA semi-distance between them is defined as:

NSD f (x, y) = max S f (x| + y), S f (y| + x) max S f (x), S f (y) . ( 10
)
Proposition 13 (NSD f is a normalized semi-distance):

1) NSD f (x, y) = NSD f (y, x); 2) 0 ≤ NSD f (x, y) ≤ 1; 3) NSD f (x, y) = 0 ⇐⇒ x = y. Proof:
1) This is clear by Def. 10;

2) By Proposition 4, NSD f (x, y) ≥ 0. By Proposition 5, NSD f (x, y) ≤ 1;

3) Restricting the identity of indiscernibles to cases |x|, |y| > 1 avoids the two pathological cases for which NSD f is indeterminate (after Remark 4). We shall now write the rest of this proof using the generic factorization symbol because either factorization kind would work. Setting NSD f (x, y) = 0 implies that S f (x y) = S f (y x) = 0, in which case

x is a substring of y and y is a substring of x, so x = y. Now setting x = y implies that S f (x y) = S f (y x) = 0, in which case NSD f (x, y) = 0.

Remark 6 (Changing the conditioning): We wrote the above proof using to emphasize that assessing similarity also leads to a semi-distance:

NSD sim f (x, y) max S f (x|y), S f (y|x) max S f (x), S f (y) . (11)
Remark 7 (True distance): It is easy to find counter-examples that violate the triangle inequality for NSD f (or Eq. 11). If it is an issue, then one can use an external procedure based on Dual-Rooted Prim Trees (DRPT) that will turn a semi-distance matrix into a true distance matrix [START_REF] Galluccio | Clustering with a New Distance Measure based on a Dual-Rooted Tree[END_REF], even if renormalization may be needed afterwards. We have implemented a naive, slow version of the DRPT procedure. For the record, we shall include results using the DRPT procedure in the benchmark of Sec. IV-E.

C. Datasets description

In contrast with [START_REF] Cilibrasi | Clustering by Compression[END_REF] where the authors focus on clustering, we would like to assess the performances of universal semi-/distances for classification. Another difference with [START_REF] Cilibrasi | Clustering by Compression[END_REF] is that we are seeking a way to naturally extend the application of those metrics to datasets containing continous values describing numeric attributes, not only discrete datasets.

Of course, there should be a conversion from the real domain to discrete values. A lot of works in the literature focus on computing thresholds to binarize each value, but this approach seems to be very data-dependent and therefore hard to reproduce -not to mention that the length of the strings produced is the number of attributes. Instead, we put forward the following simple idea that for each attribute we:

1) Select a different set of two characters ;

2) Normalize each attribute value;

3) Use N characters to actually write down an attribute value.

For example, using character set {a, b} and N = 10, a value of 0.3 will produce the string aaabbbbbbb. With a high number of attributes, character sets will ultimately repeat but this approach can already handle a few of them. Table II summarizes the characteristics of the datasets we consider. Obviously, the value of N should be set according to a tradeoff between the length of the strings produced and the precision of the numerical values for the attributes.

The markov dataset is the only purely discrete dataset. It was constructed using 6 probability transition matrices between 256 characters. For each matrix, 6 realizations where obtained to produce strings of 15001 characters. This length was chosen to fit the internal 32KiB buffer of gzip so that it will handle the concatenation of two strings correctly (for the NCD/gzip). The goal for classification is to group realizations by transition matrices. This dataset is the closest to those used in [START_REF] Cilibrasi | Clustering by Compression[END_REF] and compressors are expected to perform well.

The gaussian dataset was constructed in a more involved fashion, in order to assess the classification robustness against noise. We have chosen respectively {22, 20, 25, 26, 27, 24, 24} samples of 7 classes. Each sample contains values for 1000 attributes and 1000 noisy, noninformative values. Each informative attribute value is drawn according to N (µ k , σ 2 ), where µ k ∈ {20, . . . , 26} is the mean value of each class and σ 2 = 23.5 is their common variance.

Non-informative values are drawn according to N (19, σ 2 ). In total, each column of the sample matrix contains 2000 values and the rows are randomly permuted. The standard iris and wine datasets were downloaded from the UCI Machine Learning

Repository7 [START_REF] Dua | UCI machine learning repository[END_REF]. Each column contains respectively 4 and 13 attributes. The iris dataset is often advertised as being easily classified. Yet, it is anticipated that the small number of attributes it contains may represent a challenge for universal semi-/distances and the value of N may eventually be the most important parameter. For the sake of uniformity and after taking into account the various numerical precision of the attributes, we set N = 51 for all relevant datasets.

D. Methodology

Classification is performed on the semi-/distances matrices using spectral clustering. In order to compute accuracy figures (see below), the true and predicted labels are matched iteratively by selecting the most probable label each time.

Remark 8 (On the NCD): It turns out that the matrices produced by the NCD 8 are neither symmetric nor perfectly normalized. This does not come as a surprise (at least for the lack of normalization), as it was reported in [START_REF] Cilibrasi | Clustering by Compression[END_REF].

Yet, this prevents direct use of spectral clustering. Therefore, the NCD matrices are symmetrized by taking the maximum value of NCD(i, j) and NCD( j, i) and they are further properly normalized by using the largest value in the matrix.

For each dataset, we have computed the NCD using four standard compressors: gzip [START_REF] Deutsch | DEFLATE Compressed Data Format Specification, version 1.3[END_REF], [START_REF]GZIP File Format Specification[END_REF], lzma [START_REF] Pavlov | 7z Format[END_REF], bzip2 [START_REF] Seward | [END_REF] and ppmd [START_REF] Shkarin | PPM: One Step to Practicality[END_REF]. gzip and lzma are both based on LZ77: gzip uses a 32KiB sliding window with Huffman coding of the characters, the quantized lengths and distances, while lzma uses a default 1MiB sliding window with arithmetic coding of the characters, the quantized lengths and distances, intertwined with Move-To-Front of the distances. bzip2 is based on the BWT of 256KiB-blocks of the input string, encoded using Huffman coding of the characters, processed with Move-To-Front. ppmd uses prediction by partial matching to build a model of the next character given the previous ones, producing socalled context and transition structures that are encoded using range coding. Both bzip2 and ppmd are noted for performing best on natural language or source code inputs, while gzip and lzma offer more uniform performances -with lzma being the current state of the art in general-purpose compression.

E. Results

We shall evaluate the NSD's and the NCD's against the following criterions:

1) accuracy: the ratio of correctly classified samples ;

2) silhouette coefficient [START_REF] Rousseeuw | Silhouettes: A Graphical Aid to the Interpretation and Validation of Cluster Analysis[END_REF]: a measure of cohesion (the mean distance of a sample to others in the same cluster) vs. separation (the minimum distance between a sample and those in other clusters) -the highest the silhouette coefficient value, the cleanest the clustering ;

3) wall-clock time: the NCD relies on the potential parallelism implemented at the compression routine level (to the best of our knowledge, only lzma implements some form of parallelism), while the NSD natively breaks down all factorizations so they can be distributed on the available CPU cores. In this regard, since we did not try to rewrite the NCD to take advantage of parallelism, it should be noted that most running times of the NCD are actually worst-case figures.

Table III reports the overall results obtained for the four datasets and the ten semi-/distances we have benchmarked.

As for the accuracy, the NSD's gives more consistent results than the NCD: the NCD may outperform the NSD (only for the gaussian dataset, NCD/gzip outperforms NSD sim -yet not the NSD), but at the price of selecting the compressor best suited to the data at hand [START_REF] Cebrián | Common Pitfalls Using the Normalized Compression Distance: What to Watch Out for in a Compressor[END_REF] and with the additional issue that a given compressor may produce catastrophic results (gzip for markov or wine, lzma for gaussian or ppmd for iris). For the NSD, the following observations are in order: 1) the DRPT routine surprinsingly decreases the classification performances. Actually, the DRPT routine works by inferring the manifold on which the data lives -and we likely have not enough samples for this strategy to work ;

2) comparing NSD vs. NSD sim , the NSD consistently produces clusters with higher silhouette coefficients, yet the accuracy is sometimes better using NSD sim ;

3) comparing NSD 1 vs. NSD, the results are nearly identical -the strings are already long enough so that the added precision of S p does not help here. Also for the silhouette coefficient, the NSD's gives more consistent and better results than the NCD, meaning that the NSD clusters are more compact and distant from one another. Again, results vary greatly for the different compressors, with lzma producing better clusters more often than other compressors.

Up to the reservation explained above, the running time is better for the NCD on small datasets (with the notable exceptions of lzma and ppmd being consistently the slowest performers because of their entropy coding stage based on arithmetic coding). When the size of the dataset increases (markov and gaussian), the built-in parallelism of the NSD starts to outperform the various NCD's, often by orders of magnitude. Using NSD instead of NSD 1 is hardly noticeable from the running time point of view.

For the record, we include in Fig. 3 the phylogenetic trees obtained from computing both NSD sim and NSD on the language dataset. This figure reflects the general conclusion of Table III that measuring dissimilarity will often produce more compact and distant clusters than measuring similarity.

V. CAUSAL INFERENCE

Named after the initials of its inventors, the PC algorithm [START_REF] Spirtes | Causation, Prediction, and Search[END_REF] is the standard method to perform causal inference. It was first devised for the probabilistic setting using Markovian independence on nodes of the Directed Acyclic Graph (DAG) underlying the observations. Later, it was shown that algorithmic information is also a faithful measure to assess Markovian independence on a graph [START_REF] Janzing | Causal Inference Using the Algorithmic Markov Condition[END_REF], [START_REF] Lemeire | Replacing Causal Faithfulness with Algorithmic Independence of Conditionals[END_REF].

Starting from a complete, undirected graph, the greedy PC algorithm will collapse edges based on a local, Markovian measure of independence to produce a skeleton. By far, this is the most computational-intensive part since, for each node y, one has to compute independence with each connected nodes x, conditioned on an increasing number of other connected nodes. This process is repeated until either no more edges can be collapsed, or a hard limit is reached (see [START_REF] Colombo | Order-Independent Constraint-Based Causal Structure Learning[END_REF] for more details on computing the skeleton).

The historical, population, version of the PC algorithm is noted for its sensitivity to the order in which Markovian independence is tested for edge collapse. Recently, a stable version has been proposed [START_REF] Colombo | Order-Independent Constraint-Based Causal Structure Learning[END_REF], that requires (i) the edges to be repeatedly tested in the same order and (ii) edges to be collapsed only after all candidates have been tested. This makes the stable version more robust but slower compared to the population version. We have implemented both versions.

In the first step of computing the skeleton, further information is saved for the second, final step, that will produce a Completed Partially DAG from the skeletoni.e., a graph with as much oriented edges as possible [START_REF] Pearl | Causality, Models, Reasoning, and Inference[END_REF] that is the representative of the class of DAGs the true DAG underlying the data belongs to.

A. SALZA normalized independence criterion

Hence, all that is needed to use the PC algorithm with strings is an efficient routine to compute string independence in order to build the skeleton. The routine should be efficient because, although it is greedy, the PC algorithm has a worst-case complexity of O(n 3 ) (where n is the number of strings) -the true complexity depending on the DAG underlying the observed strings. Repeatedly testing for string independence using a compressor would inevitably lead to severe inaccuracies and would require a cumbersome management of the concatenation of strings (not to mention that the search structures would have to be built from scratch all too often). This is where the flexible conditioning of SALZA really shines.

In the probabilistic setting, Markovian independence is typically assessed within 95% confidence intervals. Such a notion of confidence is tedious to translate in the algorithmic setting and we therefore choose to assess independence up to some normalized value η (see Def. 5.1). The default value for η is set to 10 -2 .

Definition 5.1 (Normalized independence criterion): For strings x, y, z 1 , . . . , z n ∈ A , |y| > 1, x and y are said to be independent of z 1 , . . . , z n up to η, 0 ≤ η ≤ 1, iff:

I f (x : y|z 1 , . . . , z n ) S f (y) ≤ η. ( 12
)
Proposition 14 (The left-hand side of Eq. 12 is normalized):

0 ≤ I f (x : y|z 1 , . . . , z n ) S f (y) ≤ 1.
Proof: By Proposition 5, S f (y|x, z 1 , . . . , z n ) ≤ S f (y|z 1 , . . . , z n ) ≤ S f (y). Propositions 4 and 8 conclude.

Remark 9 (I +

f ): As a by-product of our implementation, our code also natively implements I + f similarly to Remark 6. Proposition 14 obviously also holds for I + f . However, the meaning of I + f is less clear.

B. Causal relationships on languages

From the point of view of causality, representing languages with phylogenetic trees like in Fig. 3 implicitly assumes a tree-like causal structure. Such an assumption may well be regarded as too restrictive, hiding the influences several languages may have had on another.

We include on Fig. 4 the resulting DAG after applying the PC algorithm on the languages dataset. Here, it is worth mentioning that the texts that are included in this dataset are translations of the Universal Declaration of Humain Rights. This may explain that some languages are unexpectedly connected: the lexical field of this text is mostly legal, so that some languages may have borrowed terms from another. For example, the English legal lexicon contains a few words from the French.

Lowering η will produce less connected components but will generate more information to orient the edges. We have observed that a few causal structures are preserved while varying η.

However, the very idea of using a global value for η may be seen as a limitation of the PC algorithm: further investigation should concentrate on selecting an adequate threshold locally.

  is a practical embodiment of the NID that computes a universal distance between sequences x and y. If xy is the concatenation of x and y and |x| is the length of x, then: NCD(x, y) = |C(xy)| -min{|C(x)|, |C(y)|} max{|C(x)|, |C(y)|} .

Remark 2 (

 2 Integer composition of |y|): L y x 1 ,...,x n is an integer composition of |y| because ∑ l∈L y x 1 ,...,x n l = |y|. Proposition 1 (|L y x 1 ,...,x n | is non-increasing by conditioning): For any three strings x, y, z ∈ A , |L y x,z | ≤ |L y x |. Proof: If y = / 0, L y x = L y x,z = / 0. The following assumes y ∈ A + . Let L ≤i y x = l 1 y x , . . . , l i y x

Definition 2 . 4 (

 24 Relative similarity measure): Let [|y|] the set of partitions of integer |y| and L y x 1 ,...,x n ∈ [|y|]. A relative similarity measure between y and prior knowledge strings x 1 , . . . , x n is defined as a score function S : [|y|] → R + ∪ {0}.

Proposition 2 (

 2 Score compliant with breadth-first order): Let [|y|] the set of integer partitions of |y|, ⊥ = {1, . . . , 1} (|⊥| = |y|) and S : [|y|] → R + ∪ {0} a score function that verifies the following conditions:

  on an integer partition of |y|. For actual compression, of course, we could not encode the factors sorted by decreasing length, but the estimated compressed string length would be the same. Since real LZ77-based compressors tend to transmit characters as such (and not as length-1 references), we shall spend log |A| bits for each of them. Let d = max {|y|, |x 1 |, . . . , |x n |}. Using δ-codes[START_REF] Elias | Universal Codeword Sets and Representations of the Integers[END_REF]

Proof: 1 )

 1 (Normalization) By Def. 2.5; 2) (Monotonicity over prefixes) Using the lexicographic order on strings, x ≤ y implies that one can write y = xz and therefore |y| = |x| + |z|. Hence, depending on the kind of factorization used: • Case |: L y| / 0 = L xz| / 0 = L x| / 0 + L z|x , • Case | + : L y| + / 0 = L xz| + / 0 = L x| + / 0 + L z| + / 0 . Monotonicity over prefixes trivially holds for S 1 in both cases (L y / 0 can only have more elements than L x / 0 ). For S p , denote by l (x) k (resp. l (y) k ) the (sorted) elements of L x / 0 (resp.

Fig. 2 :

 2 Fig. 2: Heatmap of the departure from symmetry. The mean (resp. maximum) departure is 0.24% (resp. 3.95%.)

Fig. 3 :

 3 Fig. 3: Phylogenetic trees computed on the languages datasets [9], [10] using the Ward distance (top: NSD sim , bottom: NSD).

Fig. 4 :

 4 Fig.4: Stable skeleton computed from the languages dataset. η = 1% (time: 9m44s.) We did not represent the CPDAG here because only three edges were actually oriented.

  Population skeleton, η = 0.02.

  Population skeleton, η = 0.01.

  Stable skeleton, η = 0.01.

  Stable skeleton, η = 0.009.

  Population skeleton, η = 0.008.

  Stable skeleton, η = 0.008.

Fig. 5 :

 5 Fig.5: Causal inference on the Toussaint drafts: population (resp. stable) skeleton and CPDAG mean computation time: 0.339s (resp. 0.405s).

  , Vol. 1, solution of Exercise 2.3-(15) p. 562 and Sec. 2.3.3 p. 351). Assuming Conditions (1) and (

TABLE II :

 II Datasets for classification performance assessment (N = 51).

TABLE III :

 III Classification performance of universal semi-/distances. The best results are in bold face and the worst are in italic face. The wall-clock time for DRPT variants is reported as the time the naive DRPT routine took in addition to that of the corresponding NSD semi-distance matrix computation. The NSD used 19 CPU cores.

The Kolmogorov prefix complexity has become the standard version of algorithmic complexity (see[START_REF] Li | An Introduction to Kolmogorov Complexity and Its Applications[END_REF], Chap. 3) because, relying on a Turing machine accepting programs that are not the prefix of others, it is able to circumvent issues that the original, plain Kolmogorov complexity did not consider at the time (subadditivity and non-monotonicity over prefixes), see infra.

The number of partitions of an integer is finite, therefore operation A or B cannot always be applied. This is made explicitly clear in the definition of operation B, but it is more implicit in that of operation A (which requires two summands of value 1 to be applied). Thus, the partition tree is not a full binary tree, see Fig.1.

The radix order is defined by x < y if (|x| < |y|) or (|x| = |y| and x = uAx and y = uBy ).

The z k 's in Def. 2.1 start at 1, so we do not extend the δ-code for offsets in s k 's to handle zeros.

For the origin of this name, see[START_REF] Wyner | A Definition of Conditional Mutual Information for Arbitrary Ensembles[END_REF].

http://archive.ics.uci.edu/ml/datasets

Software downloaded from https://complearn.org.

In Fig.5, most nodes are connected with two arrows, giving the wrong impression of short cycles. This is indeed a limitation of GraphViz wich allows either for directed or undirected graphs, having no provision for graphs containing both types of edges. When two nodes are connected to each other with two arrows, it only means the algorithm left the edge undirected.

This mandatory feature we need is even called an "undesirable aspect" of CPS1 (see[START_REF] Chen | Lempel-Ziv Factorization using Less Time and Space[END_REF], Sec. 2.2).

ACKNOWLEDGEMENTS

We are indebted to the anonymous reviewers for providing many helpful comments and pointing us to [13].

The Toussaint drafts are courtesy Prof. T. Lebarbé (Univ. Grenoble Alpes and TGIR Huma-Num).

Marion Revolle was supported by a French Ministry of Research PhD grant. GIPSA-Lab, CNRS UMR 5216.

C. An experiment in literature

Jean-Philippe Toussaint is a Belgian author who works by typesetting each paragraph, which then gets printed, manually annotated for improvements, typeset again and so on until the final version is reached. We had access to the transcripts of eight successive drafts of one single paragraph in Toussaint's book La Réticence. The transcripts each contain approximately 1600 characters, with great changes between the first and only tiny improvements in the last ones.

We depict in Fig 5 the CPDAGs produced by the PC algorithm for its two variants of the computation of the skeleton (population and stable), as well as for various values of η 9 . The PC algorithm occasionally gets the edges wrong (e.g., transcript 3 is sometimes seen as causing transcript 1, which is impossible), yet it generally finds the correct causal dependencies when the changes are smaller. Following these results, one may formulate the hypothesis that the second transcript is an abandoned trail in the early steps of the writing process followed by Toussaint.

SOFTWARE

The SALZA reference implementation (C code for GNU/Linux and Debian amd64 binaries repository) and the datasets used in this work are available from: https://forge.uvolante.org/code/salza/wikis APPENDIX A

IMPLEMENTATION NOTES

A. Adjusting the precision p

It is clear from Eq. 3 that the maximum useful value for p is a function of the machine epsilon ε and |y|. Hence, p is selected such that:

.

With standard IEEE-754 binary64 floating point double's, ε = 2.22 × 10 -16 . Therefore, SALZA will take into account at least the first three l k 's for strings of size |y| ≤ 1.3 × 10 5 , and it will use only l 1 starting from strings of size |y| ≥ 6.7 × 10 7 . Of course, |L| is the dominant part in S p (L). Yet, SALZA will automatically allocate as much precision as possible for shorter strings, where it is needed more.

B. Basic algorithms and search structures

Because SALZA looks for the next longest substring in a set of strings, it requires to be able to efficiently do so from an arbitrary position in the string y being factorized.

In the vast literature pertaining to computing the Lempel-Ziv factorization, much of the works take advantage of incrementally factorizing a string (meaning that the next Lempel-Ziv factors are known only for positions at which they start in the string being factorized). However, one such algorithm exists [START_REF] Chen | Lempel-Ziv Factorization using Less Time and Space[END_REF] that computes the next Lempel-Ziv factor from an arbitrary position: CPS1 10 (in linear time if the entire string were to be factorized). This algorithm relies on the suffix array of the string being factorized, which can be constructed in linear time and uses 8 bytes of memory per string character [START_REF] Kärkkäinen | Linear Work Suffix Array Construction[END_REF]. It also makes use of the longest common prefix array, the construction and storage of which share roughly the same time and space performances [START_REF] Kasai | Linear-Time Longest-Common-Prefix Computation in Suffix Arrays and Its Applications[END_REF] as that of the suffix array.

The suffix array is the only data structure required to efficiently perform cross-factorization of y given x, also from an arbitrary position in y [START_REF] Hoobin | Relative Lempel-Ziv Factorization for Efficient Storage and Retrieval of Web Collections[END_REF]. Hence, upon start-up, SALZA builds the suffix array of all strings, saving the construction of the longest common prefix array in cases it is not needed.

C. Parallel computations

Various methods have been developed to compute the Lempel-Ziv factorization in parallel (e.g. [START_REF] Shun | Practical Parallel Lempel-Ziv Factorization[END_REF]), but they often lack the flexibility we need and their granularity sometimes is too fine-grained for CPU threads.

Instead, we have devised the following strategy. Each string to be factorized is divided in blocks. Each block factorization parameters (block start, strings involved, factorization kind for each string, etc.) are prepared and stored into a work queue. Once the work queue is ready, threads keep picking up factorization parameters for the next block and the work queue eventually gets empty. In order not to freeze the machine, we leave one core free for the other processes.

At this point, most of the work has been done in parallel. But it may happen that the Lempel-Ziv factors overlap block boundaries (the last factor of a block spans over the beginning of the next block). Hence, we perform a few additional calls to the main factorization routine (starting from the end of the last factor in a block), until the factorizations of the two consecutive blocks are synchronized. In practice, such a "stitching" overhead typically means computing less than 3 additional factors per block. This way, it is guaranteed that the resulting factorization is the same as if it were performed sequentially. Note that it is crucial for this strategy to be able to compute Lempel-Ziv factors from an arbitrary position.

Among the applications listed, computing the NSD f matrix is done even more efficiently in parallel since all computations can be scheduled right from the start (which is less easier with the PC algorithm for example). In this case, the work queue is constructed only once. Whereas for all other applications, parallelism takes place mostly at the whole string factorization level (as described above). For example, computing one conditional mutual information measure leads to the creation of a work queue containing the parameters for all blocks of each of the two parts of Eq. 7 and we can still benefit from parallel computations in performing causal inference.