Kolmogorov Complexity
The Information in One Object

Toss a fair coin 32 times. Each string below has the same probability, 2⁻³²:

A: 01010101010101010101010101010101
B: 00000000000000000000000000000000
C: 10011101001011100010110111010010

Yet everyone says C "looks random" and A and B do not. Probability cannot explain this. Every single string is equally unlikely.

Descriptions as programs

prog = 'print("ab" * 50000)'
print(len(prog), "chars of program for", len(samples["'ab' repeated"]), "bytes of output")

19 chars of program for 100000 bytes of output

Python is a fine universal machine for intuition. The theory uses a fixed universal Turing machine instead. The invariance theorem says the choice barely matters.

A handy trick. To describe x, name a set S that holds it. Then give x's position (index) inside S.

cost ≈ K(S) + log₂|S| bits

The set is the pattern. The index is the leftover noise. Here n = 1000, and c is a small fixed cost for the decoder.

Objects other than strings

Encode them as strings first. A number n is its binary form. A pair is a self-delimiting encoding of both parts. Then K(n), K(x, y) and K(graph) all make sense.

Examples (up to an additive constant)

Logs are base 2. "± c" or O(1) means a constant that depends on U but not on x.

Proof

1. Take a shortest M-description p of x. So M(p) = x and |p| = C_M(x).
2. By universality, U(⟨M⟩ p) = M(p) = x.
3. So C_U(x) ≤ |⟨M⟩| + |p| = C_M(x) + c_M, with c_M = |⟨M⟩|.
4. For two universal machines, apply this in both directions. ∎

One detail: U must know where ⟨M⟩ ends and p starts. So the machine codes ⟨M⟩ are chosen prefix-free.

C is not monotone in prefixes. A long string can be simpler than one of its prefixes. For example, 0ⁿ with n = 2^2k is simpler than most shorter runs of zeros.

Proof (pigeonhole)

1. A description shorter than n − c bits is one of 1 + 2 + 4 + … + 2^n−c−1 = 2^n−c − 1 strings.
2. Each description yields at most one output.
3. So at most 2^n−c − 1 strings get a short description. ∎

With c = 0: at least one string of every length has C(x) ≥ n. That string is truly incompressible.

def fraction_compressible(n, c):
    """Fraction of n-bit strings that have a description of < n - c bits."""
    programs = 2 ** (n - c) - 1          # all strings of length 0 .. n-c-1
    return programs / 2 ** n
for c in [1, 2, 8, 10, 20]:
    print(f"save more than {c:>2} bits: at most {fraction_compressible(64, c):.7f} of 64-bit strings")

save more than  1 bits: at most 0.5000000 of 64-bit strings
save more than  2 bits: at most 0.2500000 of 64-bit strings
save more than  8 bits: at most 0.0039062 of 64-bit strings
save more than 10 bits: at most 0.0009766 of 64-bit strings
save more than 20 bits: at most 0.0000010 of 64-bit strings

Recipe. Pick an incompressible object. Assume the claim fails. Use the failure to describe the object in too few bits. Contradiction.

Claim. No DFA accepts L = { 0ⁿ1ⁿ : n ≥ 0 }.

The idea in one line: a DFA's memory is a constant number of bits. The language needs it to remember a number that is not a constant.

Compare the pumping lemma. It proves the same thing, but you must juggle "for all splits". The incompressibility proof has no case analysis.

Proof

1. Suppose program Kc computes K.
2. Write Berry(L): "try strings in order, return the first x with Kc(x) > L". One exists, by counting.
3. Its size is |Kc| + log L + c. So the output x has K(x) ≤ log L + c′.
4. But K(x) > L by choice. For large L, L > log L + c′. Contradiction. ∎

With a computable proxy, the search is fine

from itertools import product
def first_hard(k):
    """First bit-string (by length, then lexicographic) whose zlib size is >= k bytes."""
    for n in range(1, 64):
        for bits in product(b"01", repeat=n):
            s = bytes(bits)
            if len(zlib.compress(s, 9)) >= k:
                return s
print(first_hard(20).decode())

00001011100110

This is no paradox. The program above is short, and it describes its output. That just proves zlib size is not K. zlib calls a 14-character bit string "20 bytes complex".

Halting ⇒ K

For L = 0, 1, 2, …: list every program of length L. Ask the oracle which ones halt. Run just those. The first L where some program outputs x is K(x).

K ⇒ Halting

Fix n and let m = n + c. Using K, list every string with K(x) < m (a finite set). Dovetail all programs until each of those strings has been printed. Call that time T. Suppose some program p with |p| < n halts at time t > T. Then "run p to learn t; run every program for t steps; print the first string of length m not yet printed" is a program of length n + O(1). Its output has K ≥ m = n + c. For large c that is a contradiction. So every halting program shorter than n halts by time T. Run them all for T steps and you know which ones halt.

Why bother?

• K(x, y) ≤ K(x) + K(y) + O(1), with no log term.
• ∑_x 2^−K(x) ≤ 1. So 2^−K(x) acts like a probability. This gives Ω and the coding theorem.

def elias_gamma(n: int) -> str:
    b = bin(n)[2:]
    return "0" * (len(b) - 1) + b

def elias_delta(n: int) -> str:
    b = bin(n)[2:]
    return elias_gamma(len(b)) + b[1:]

for n in [1, 2, 5, 17, 1000]:
    print(f"{n:>5}  gamma={elias_gamma(n):<20} delta={elias_delta(n)}")

codes = [elias_gamma(n) for n in range(1, 2000)]
print("prefix-free:", is_prefix_free(codes),
      " Kraft sum:", round(sum(2 ** -len(w) for w in codes), 6))

    1  gamma=1                    delta=1
    2  gamma=010                  delta=0100
    5  gamma=00101                delta=01101
   17  gamma=000010001            delta=001010001
 1000  gamma=0000000001111101000  delta=0001010111101000
prefix-free: True  Kraft sum: 0.999489

Gamma costs 2 log n + 1 bits. Delta costs log n + 2 log log n + O(1). is_prefix_free sorts the words and checks neighbours. The Kraft sum tends to 1 as more words are added.

Gamma: say the length in unary, then the number

Write n in binary. If it has L bits, put L − 1 zeros in front.

Kraft sum for 1..7: 1/2 + 2/8 + 4/32 = 0.875 ≤ 1. Good.

Delta: say the length in gamma

13 = 1101 has 4 bits. So write gamma(4) = 00100, then the bits after the leading 1: 101. Result 00100101 (8 bits). Gamma used 7 here, but delta wins for large n.

Levin calls m the "a priori probability". Solomonoff used the same idea for prediction (see later).

Proof sketch

• Left: the programs {p_x} form a prefix-free code. Shannon's source coding theorem says no prefix code beats H on average.
• Right: build a Shannon–Fano code for P, with lengths ⌈−log P(x)⌉. A program with ⟨P⟩ can decode it. So K(x) ≤ −log P(x) + K(P) + O(1). Take the average.

For n i.i.d. draws, divide by n. The K(P)/n term vanishes. Per symbol, expected K equals the entropy rate.

def H(p): return -p * math.log2(p) - (1 - p) * math.log2(1 - p)

print(f"{'p':>5}{'H(p)':>8}{'lzma bits/bit':>15}")
for p in [0.5, 0.25, 0.1, 0.01]:
    data = bernoulli_bytes(p, 400_000, seed=7)
    rate = 8 * len(lzma.compress(data, preset=9)) / 400_000
    print(f"{p:>5}{H(p):>8.3f}{rate:>15.3f}")

    p    H(p)  lzma bits/bit
  0.5   1.000          1.001
 0.25   0.811          0.838
  0.1   0.469          0.502
 0.01   0.081          0.107

bernoulli_bytes(p, n, seed) packs n biased coin flips into bytes. lzma is a general tool, not tuned for this source. It still lands within about 0.03 bits of the entropy.

Compare Shannon: H(X,Y) = H(X) + H(Y|X). The algorithmic version is the same law, for single objects, with the small K(x) correction.

Estimate it with a compressor

def C_given(x: bytes, y: bytes) -> int:
    """C(x | y) ~ C(y + x) - C(y)."""
    return C(y + x) - C(y)
x = docs["dnaA'"]
print("C(x) =", C(x), "  C(x | dnaA) =", C_given(x, docs["dnaA"]),
      "  C(x | dnaB) =", C_given(x, docs["dnaB"]))

C(x) = 1178   C(x | dnaA) = 631   C(x | dnaB) = 1109

• dnaA' is dnaA with 10% of its letters changed.
• Knowing dnaA nearly halves the cost. What is left is mostly the mutations.
• Knowing an unrelated dnaB barely helps.

C(x) is the smallest of zlib, lzma and bz2 output, in bytes. The docs are built on the NCD slide. The trick "compress y then x" is the idea behind NCD.

Real uses: mtDNA trees of mammals, language family trees from the UN Human Rights text, music genres, plagiarism and malware families, and anomaly detection.

def Z(x: bytes) -> int:
    return len(lzma.compress(x, preset=9 | lzma.PRESET_EXTREME))

def ncd(x: bytes, y: bytes) -> float:
    zx, zy, zxy = Z(x), Z(y), Z(x + y)
    return (zxy - min(zx, zy)) / max(zx, zy)

def mutate(s: bytes, rate: float, seed: int) -> bytes:
    """Replace a fraction `rate` of bytes with other bytes from s itself."""
    r = random.Random(seed)
    return bytes(r.choice(s) if r.random() < rate else c for c in s)

r = random.Random(1)
dna_a = bytes(r.choice(b"ACGT") for _ in range(4000))
dna_b = bytes(r.choice(b"ACGT") for _ in range(4000))
eng  = b" ".join(r.choice([b"it was the best of times", ...]) for _ in range(150))
code = b"\n".join(r.choice([b"def f(n):", b"    if n < 2:", ...]) for _ in range(150))
docs = {"dnaA": dna_a, "dnaA'": mutate(dna_a, 0.10, 2),
        "dnaB": dna_b, "dnaB'": mutate(dna_b, 0.10, 3),
        "eng":  eng,   "eng'":  mutate(eng, 0.10, 4),
        "code": code,  "code'": mutate(code, 0.10, 5)}
names = list(docs)
D = {(a, b): ncd(docs[a], docs[b]) for a in names for b in names if a != b}

Each eng and code list has 8 phrases or lines. The full lists are in the test file.

         dnaA  dnaA'   dnaB  dnaB'    eng   eng'   code  code'
dnaA        -   0.29   0.84   0.83   0.95   0.95   0.95   0.95
dnaA'    0.29      -   0.84   0.83   0.95   0.95   0.95   0.95
dnaB     0.84   0.84      -   0.28   0.95   0.95   0.95   0.95
dnaB'    0.84   0.84   0.25      -   0.95   0.95   0.95   0.95
eng      0.93   0.93   0.93   0.93      -   0.48   0.75   0.91
eng'     0.93   0.93   0.92   0.93   0.83      -   0.91   0.87
code     0.93   0.93   0.92   0.93   0.74   0.91      -   0.47
code'    0.93   0.93   0.92   0.93   0.90   0.86   0.81      -

NCD(x,y) = (Z(xy) − min{Z(x),Z(y)}) / max{Z(x),Z(y)}

Read it as: "How many new bytes does y add once x is known, compared to the bigger file?" Sizes are lzma bytes from the demo.

• dnaA, dnaA': 90% of the copy is shared. lzma only pays for the 10% of changed letters.
• dnaA, eng: nothing to reuse. Z(xy) ≈ Z(x) + Z(y), so NCD is near 1.
• eng, code: both are made of short English-like words. Some phrases help each other.

def single_link(names, D):
    clusters = [[n] for n in names]
    while len(clusters) > 1:
        i, j = min(((i, j) for i in range(len(clusters))
                           for j in range(i + 1, len(clusters))),
                   key=lambda ij: min(D[a, b] for a in clusters[ij[0]]
                                               for b in clusters[ij[1]]))
        d = min(D[a, b] for a in clusters[i] for b in clusters[j])
        print(f"merge {clusters[i]} + {clusters[j]}  at {d:.2f}")
        clusters[i] = clusters[i] + clusters[j]; del clusters[j]
single_link(names, D)

merge ['dnaB'] + ["dnaB'"]  at 0.28
merge ['dnaA'] + ["dnaA'"]  at 0.29
merge ['code'] + ["code'"]  at 0.47
merge ['eng'] + ["eng'"]  at 0.48
merge ['eng', "eng'"] + ['code', "code'"]  at 0.75
merge ['dnaA', "dnaA'"] + ['dnaB', "dnaB'"]  at 0.83
merge ['dnaA', "dnaA'", 'dnaB', "dnaB'"] + ['eng', "eng'", 'code', "code'"]  at 0.95

What makes an infinite 0/1 sequence random? Frequency rules like "half zeros" are not enough. 0101010… passes that test. We want a sequence to pass every effective test.

• Each test catches a set of measure 0. There are only countably many tests. So almost every sequence is ML-random.
• There is a single universal test that catches everything any test catches.

So δ is a universal test statistic. A deficiency of 30 bits is like a p-value of 2⁻³⁰. And it holds for every computable test at once, not just the one you picked.

def deficiency_bits(x: bytes) -> int:
    """n - C(x) in bits: a lower bound on the true deficiency n - K(x), up to O(1)."""
    return 8 * len(x) - 8 * C(x)

tests = {
    "fair coin":      bernoulli_bytes(0.5, 80_000, seed=11),
    "biased p=0.1":   bernoulli_bytes(0.1, 80_000, seed=12),
    "period 7":       (b"\x5a\x13\x77\x00\xff\x42\x99" * 1500)[:10_000],
}
for name, x in tests.items():
    print(f"{name:<14} n={8*len(x)} bits   n - C(x) = {deficiency_bits(x):>6} bits")

fair coin      n=80000 bits   n - C(x) =    -88 bits
biased p=0.1   n=80000 bits   n - C(x) =  36128 bits
period 7       n=80000 bits   n - C(x) =  79608 bits

• Fair coin: the compressor loses 88 bits of header. We find no evidence against "uniform", as expected.
• Biased coin: 36128 bits of evidence that this is not a fair coin. That fits n(1 − H(0.1)) ≈ 42500, minus lzma's overhead.
• Period 7: almost all of it is structure.

Ω depends on U. Calude, Hertling, Khoussainov, Wang and Kučera–Slaman showed that the left-c.e. ML-random reals are exactly the Ωs of the different universal machines.

Proof (Berry again)

1. For a number L, write program P_L. It lists all proofs of F in order. It prints the first x with a proof of "K(x) > L".
2. |P_L| ≤ log L + c_F, since F's axioms are a fixed program.
3. If it halts, K(x) ≤ log L + c_F. By soundness, also K(x) > L.
4. For L > log L + c_F, that is a contradiction. So P_L never halts. ∎

In practice, L_F is a few thousand bits for ZFC. You cannot prove that any specific megabyte file is incompressible, even though nearly all of them are.

• Ideal MDL uses K(M) + K(D | M). Practical MDL restricts to a model class. It uses (k/2) log n bits for k real parameters.
• The (k/2) log n penalty equals BIC. Better versions include normalized maximum likelihood and Bayesian mixture codes.
• Kolmogorov structure function: the best "sufficient statistic" of x is the smallest model that makes x look like a typical member.

Data: 20,000 bits from a source that repeats its last bit with probability 0.9. That is order-1 Markov. Which order k does MDL pick?

s = markov_source(20_000, seed=3)
print(f"{'k':>2}{'L(model)':>10}{'L(data|model)':>15}{'total':>10}")
best = None
for k in range(6):
    m, d = two_part_bits(s, k)
    print(f"{k:>2}{m:>10.0f}{d:>15.0f}{m + d:>10.0f}")
    if best is None or m + d < best[1]: best = (k, m + d)
print("MDL picks k =", best[0])

 k  L(model)  L(data|model)     total
 0         7          20000     20007
 1        14           9307      9321
 2        29           9304      9333
 3        57           9299      9356
 4       114           9293      9408
 5       229           9279      9508
MDL picks k = 1

two_part_bits counts next bits per context. It charges the empirical entropy for the data and 2^k · ½ log n for the model. Higher orders fit a little better, but not enough to pay for their parameters.

It is uncomputable, like K. AIXI (Hutter) adds actions and rewards on top. Practical versions limit the model class, as in our toy.

A toy: all "repeat this pattern" programs

def solomonoff_next(observed: str, max_period=12):
    """Mix all 'repeat pattern p' programs, prior 2^-(len(code))."""
    w = {"0": 0.0, "1": 0.0}
    for per in range(1, max_period + 1):
        for bits in product("01", repeat=per):
            pat = "".join(bits)
            gen = (pat * (len(observed) // per + 2))
            if gen.startswith(observed):
                prior = 2.0 ** -(2 * per + 1)   # per ones, a zero, then per bits
                w[gen[len(observed)]] += prior
    z = w["0"] + w["1"]
    return {b: round(v / z, 4) for b, v in w.items()}

for obs in ["0", "01", "0101", "011011", "0110110110"]:
    print(f"{obs:<11} -> P(next) = {solomonoff_next(obs)}")

0           -> P(next) = {'0': 0.7501, '1': 0.2499}
01          -> P(next) = {'0': 0.7501, '1': 0.2499}
0101        -> P(next) = {'0': 0.9723, '1': 0.0277}
011011      -> P(next) = {'0': 0.9925, '1': 0.0075}
0110110110  -> P(next) = {'0': 0.0001, '1': 0.9999}

After "0", the period-1 pattern "repeat 0" (a 3-bit code) has the biggest weight. After "0110110110", only period 3 and its multiples still fit, plus a few long patterns. They all but agree on "1".

K ranks a random string as the most complex. But a random string is boring. A DNA genome or a proof feels "complex" in a different way. It is structured and hard to make.

Related: Gell-Mann & Lloyd's effective complexity, and the Kolmogorov structure function.

Both measures depend on a significance parameter s or c. Both are even less computable than K.

def C(x: bytes) -> int:
    """Best of three real compressors, in bytes. An upper bound on K(x) + O(1)."""
    return min(len(zlib.compress(x, 9)), len(lzma.compress(x, preset=9)), len(bz2.compress(x, 9)))

samples = {
    "zeros":        bytes(N),
    "'ab' repeated": b"ab" * (N // 2),
    "English-ish":  (b"the quick brown fox jumps over the lazy dog " * 3000)[:N],
    "pi digits":    pi_digits(N).encode(),
    "random bytes": rng.randbytes(N),
}
for name, x in samples.items():
    print(f"{name:<15}{len(x):>8}{len(zlib.compress(x, 9)):>8}"
          f"{len(lzma.compress(x, preset=9)):>8}{len(bz2.compress(x, 9)):>8}")
x = samples["pi digits"]
print(f"pi: {8 * C(x) / len(x):.3f} bits/digit  vs  log2(10) = {math.log2(10):.3f}")

input               raw    zlib    lzma     bz2
zeros            100000     120     148      47
'ab' repeated    100000     122     152      43
English-ish      100000     358     192     171
pi digits        100000   48174   43732   43229
random bytes     100000  100041  100064  100808
pi: 3.458 bits/digit  vs  log2(10) = 3.322

N = 100_000, rng = random.Random(42). pi_digits uses Machin's formula with Python integers. It is 10 lines, in the test file.