· 개발

HyperLogLog: 조화평균으로 고유값 개수를 추정하는 원리와 1.04/√m 오차의 조건

2007년 Flajolet 등의 원 논문을 기준으로 HyperLogLog가 레지스터 j와 ρ로 고유값 개수를 추정하는 방식, 조화평균을 쓰는 이유, 소범위·대범위 보정, 1.04/√m 오차가 전제하는 조건을 정리합니다.

algorithmsdata-structureshyperloglogcardinality-estimation

중복을 뺀 원소의 개수를 카디널리티(cardinality)라고 합니다. 논문이 드는 예로는 네트워크 패킷 헤더에 담긴 출발지–목적지 주소 쌍이 일정 시간 동안 서로 다른 것이 몇 개였는지가 있습니다. 정확히 세려면 이미 본 값을 기억해야 하므로, 필요한 저장 공간은 원소 수에 거의 비례합니다.

Philippe Flajolet, Éric Fusy, Olivier Gandouet, Frédéric Meunier의 2007년 논문 “HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm”(AofA 2007, DMTCS Proceedings AH권, 127–146쪽)은 m개의 작은 레지스터만으로 데이터를 한 번 훑어 카디널리티를 추정하는 알고리즘을 제시하고, 상대 오차(표준오차)가 대략 1.04/√m이라고 분석합니다. 이 글은 Flajolet의 INRIA 페이지에 공개된 이 논문의 PDF를 기준으로, 레지스터 번호 j와 ρ를 어떻게 구하는지, 왜 조화평균을 쓰는지, 1.04/√m이라는 숫자가 어떤 가정 위에 있는지, 그리고 추정값에서 끌어내면 안 되는 결론이 무엇인지 다룹니다. 이후에 나온 HyperLogLog++ 같은 변형이나 특정 데이터베이스의 구현은 다루지 않습니다.

관찰값 하나: ρ

해시값을 0과 1로 된 비트열로 보고, ρ(s)를 비트열 s에서 가장 왼쪽에 있는 1의 위치로 정의합니다. 위치는 1부터 셉니다. ρ(1…) = 1, ρ(0001…) = 4이고, 맨 앞에 이어지는 0의 개수에 1을 더한 값과 같습니다.

논문은 해시 비트가 서로 독립이고 각 비트가 1/2 확률로 1이라고 가정합니다. 그러면 ρ ≥ k일 확률은 2^(1−k)입니다(2.1절). 0^(ρ−1)1로 시작하는 값을 봤다면 고유값이 대략 2^ρ개 이상이라는 신호가 되는데, 논문도 이를 “more or less a likely indication” 정도로 조심스럽게 표현합니다. 같은 값은 해시값도 같으므로 ρ의 최댓값은 어떤 값이 몇 번 반복되었는지와 무관합니다. 카디널리티를 세는 데 이 관찰값을 쓸 수 있는 이유입니다.

관찰값 하나로는 변동이 너무 큽니다. 서로 다른 해시 함수 m개로 같은 실험을 m번 하면 평균의 표준편차가 σ/√m으로 줄겠지만, 논문은 이 방법에 두 가지 문제가 있다고 지적합니다. 원소마다 해시를 m번 계산해야 하고, 서로 독립인 해시 함수 집합을 만드는 방법은 알려져 있지 않습니다. 그래서 Flajolet–Martin의 확률적 계수(Probabilistic Counting)에서 쓴 확률적 평균(stochastic averaging)을 가져옵니다. 해시값의 앞 b비트로 m = 2^b개의 하위 스트림 중 하나를 고르고, 하위 스트림마다 관찰값을 따로 유지하는 방식입니다. 원소 하나당 해시 계산은 한 번이면 됩니다.

해시 함수 자체에 대해서는 해시 함수 글을 참고하세요.

레지스터 j와 ρ(w): 1은 어디서 더해지나

논문 그림 2의 알고리즘은 다음과 같습니다.

  1. m = 2^b(b는 양의 정수)로 두고 레지스터 M[1], …, M[m]을 −∞로 초기화합니다.
  2. 각 원소 v에 대해 x = h(v)를 계산합니다.
  3. j = 1 + ⟨x₁x₂…x_b⟩₂. x의 앞 b비트를 이진수로 읽은 값에 1을 더한 것이 레지스터 번호입니다.
  4. w = x_{b+1}x_{b+2}…. 앞 b비트를 떼어 낸 나머지에서 ρ(w)를 구하고, M[j] = max(M[j], ρ(w))로 갱신합니다.
  5. 모든 원소를 읽은 뒤 다음 절의 Z와 E를 계산합니다.

1이 더해지는 곳이 두 군데라서 구현할 때 헷갈리기 쉽습니다. j의 +1은 논문이 레지스터를 1번부터 세기 때문에 붙은 것이므로, 0부터 세는 배열이라면 앞 b비트 값을 그대로 인덱스로 쓰면 됩니다. ρ의 +1은 ρ의 정의에 들어 있으므로 인덱스 방식과 상관없이 유지해야 합니다. 또 ρ는 x 전체가 아니라 w에서 셉니다. 앞 b비트는 이미 레지스터를 고르는 데 썼기 때문입니다.

아래는 이 글이 만든 예입니다. 해시가 32비트, b = 4(m = 16), x = 0x502C0000이라고 합시다.

항목 값
x (2진수) 0101 0000 0010 1100 0000 0000 0000 0000
앞 4비트 0101 = 5
논문의 j (1부터 셈) 6
0부터 세는 배열의 인덱스 5
w (나머지 28비트) 0000 0010 1100 0000 0000 0000 0000
w 앞쪽의 0 개수 6
ρ(w) 7

해시 길이가 유한하면 w가 전부 0일 수도 있습니다. 32비트 해시를 전제로 한 논문 그림 3은 ρ(0^K) = K + 1로 정의하므로, 이 경우 ρ(w) = 32 − b + 1입니다. 그래서 레지스터 값은 0부터 L + 1 − log₂m 사이에 있고(L은 해시 비트 수), 논문은 32비트 해시라면 레지스터 하나에 5비트(“short bytes”)면 충분하다고 적습니다.

조화평균으로 합치기

모든 원소를 읽은 뒤에는 다음을 계산합니다(논문 식 (2), (3)).

Z   = ( Σ_{j=1..m} 2^(−M[j]) )^(−1)
E   = α_m · m² · Z
α_m = ( m · ∫₀^∞ ( log₂((2 + u)/(1 + u)) )^m du )^(−1)

m·Z는 m개의 값 2^M[j]의 조화평균입니다. 논문이 제시하는 직관은 이렇습니다. 고유값이 n개이면 하위 스트림마다 대략 n/m개씩 들어가므로 M[j]는 log₂(n/m)에 가깝고, 2^M[j]의 조화평균은 n/m 정도가 됩니다. 여기에 m을 한 번 더 곱한 m²Z가 n 정도이고, α_m은 m²Z에 남는 체계적인 곱셈 편향을 바로잡는 상수입니다.

조화평균을 고른 이유로 논문은 다음을 듭니다. 같은 관찰값을 쓰는 이전 알고리즘 LOGLOG는 사실상 기하평균에 해당하는 방식으로 레지스터를 합쳤습니다. 조화평균은 오른쪽 꼬리가 천천히 줄어드는 분포를 길들이는 효과가 있어서 여기서 분산을 줄이는 장치로 작동하고, 이 착상은 논문이 인용한 Chassaing과 Gérin의 노트에서 왔다고 밝힙니다. 이 글의 설명을 덧붙이면, 어떤 레지스터에 우연히 큰 ρ가 들어가면 2^M[j]는 매우 커지지만 조화평균에서는 그 역수 2^(−M[j])가 0에 가까워질 뿐이라 전체 값에 미치는 영향이 작습니다.

논문 그림 1은 LOGLOG의 정확도를 1.30/√m, HyperLogLog를 1.04/√m으로 비교하고, 같은 정확도를 원래 메모리의 64%로 낼 수 있다고 적습니다. 오차가 1/√m에 비례하므로 같은 오차에 필요한 m은 (1.04/1.30)² = 0.64배라는 계산(이 글의 확인)과 맞습니다.

정리 1이 말하는 것과 그 전제

논문의 분석(2절, 3절)은 해시값들의 모음이 “이상적인 다중집합(ideal multiset)“이라고 가정합니다. 정의 1에 따르면 이는 [0, 1] 구간 위의 독립이고 균등한 확률변수 n개에 임의의 복제와 순서 바꾸기를 적용해 얻은 수열입니다. 논문은 이것이 잘 설계된 해시 함수의 결과를 모형화하는 자연스러운 방법이라고 설명하고, 해시를 무한한 비트열로 보되 실제로는 32비트 해시로 10⁹이 넘는 카디널리티를 추정하기에 충분하다고 덧붙입니다.

이 모형에서 레지스터 m ≥ 3개로 얻은 추정값 E에 대해 정리 1은 다음을 말합니다.

  • (i) 편향: n → ∞일 때 E_n(E)/n = 1 + δ₁(n) + o(1)이고, m ≥ 16이면 |δ₁(n)| < 5·10⁻⁵입니다.
  • (ii) 표준오차: n → ∞일 때 √V_n(E)/n = β_m/√m + δ₂(n) + o(1)이고, m ≥ 16이면 |δ₂(n)| < 5·10⁻⁴입니다. β₁₆ ≈ 1.106, β₃₂ ≈ 1.070, β₆₄ ≈ 1.054, β₁₂₈ ≈ 1.046이고, 극한값은 β_∞ = √(3 log 2 − 1) ≈ 1.03896입니다.

읽을 때 주의할 점이 세 가지 있습니다. 첫째, 둘 다 n → ∞에서의 점근적 결과입니다. 둘째, 표준오차는 논문 표현대로 “평균제곱의 의미에서” 전형적으로 관찰되는 상대 오차이지, 개별 추정값이 넘지 않는 한계가 아닙니다. 셋째, δ₁, δ₂는 진폭이 아주 작은 진동 함수이고, 논문은 실용적으로 무시해도 된다고 봅니다. 1.04/√m은 β_m이 m이 클 때 1.039 근처로 수렴한다는 사실에서 나온 근사입니다. 같은 분석에서 α_m은 m이 커질수록 1/(2 log 2) ≈ 0.72134에 가까워집니다(식 (29), 여기서 log는 자연로그).

증명은 원소 수가 푸아송 분포를 따른다고 놓고 기댓값을 구한 뒤(poissonization), 멜린 변환으로 점근식을 얻고, 해석적 depoissonization으로 고정된 n의 경우로 되돌리는 순서로 진행합니다. 분산 분석은 같은 방식으로 간략히 제시됩니다.

원시 추정값과 두 가지 보정

그림 2는 분석을 위한 알고리즘이고, 실제 사용을 염두에 둔 것은 그림 3의 프로그램입니다. 그림 3은 32비트 해시, 최대 카디널리티 10⁹, m = 2⁴…2¹⁶(b ∈ [4..16])에 맞춰 치수를 정했고, 그림 2와 비교해 초기화 하나와 마지막 보정 두 가지가 다릅니다. α_m으로는 다음 상수를 씁니다.

α_16 = 0.673,  α_32 = 0.697,  α_64 = 0.709,  α_m = 0.7213 / (1 + 1.079/m)  (m ≥ 128)

레지스터를 0으로 초기화하고 그림 2와 같은 방식으로 갱신한 뒤
E = α_m · m² · ( Σ_{j=1..m} 2^(−M[j]) )^(−1)          원시(raw) 추정값

E ≤ (5/2)·m 이면
    V = 값이 0인 레지스터 수
    V ≠ 0 이면 E* = m · log(m/V), 아니면 E* = E      소범위 보정
E ≤ (1/30)·2³² 이면 E* = E                          중간 범위, 보정 없음
E > (1/30)·2³² 이면 E* = −2³² · log(1 − E/2³²)       대범위 보정

초기화. 그림 2처럼 −∞로 초기화하면 빈 레지스터가 하나라도 있을 때 2^(−M[j])가 무한대가 되어 E = 0이 됩니다. 논문은 쿠폰 수집 문제를 근거로 n ≪ m log m이면 E = 0을 예상해야 하고, 따라서 작은 카디널리티에서 크게 틀린다고 말합니다. 그림 3은 0으로 초기화합니다. 이렇게 해도 n ≫ m log m이면 모든 하위 스트림이 비어 있지 않을 확률이 압도적으로 높으므로 정리 1의 결론이 여전히 적용된다는 것이 논문의 설명입니다.

소범위 보정. 0으로 초기화해도 n이 작으면 문제가 남습니다. 논문은 광범위한 시뮬레이션으로 m ≥ 16일 때 n = (5/2)m 부근에서 점근적인 영역에 실질적으로 도달하고, 그보다 작으면 비선형 왜곡이 나타난다고 적습니다. 극단적으로 n = 0이면 원시 추정값은 α_m·m ≈ 0.7m입니다. 그래서 E가 m에 비해 작을 때는 빈 레지스터 수를 씁니다. 공 n개를 상자 m개에 무작위로 던지면 빈 상자는 약 m·e^(−n/m)개이므로, 빈 상자 V개를 관찰했다면 n ≈ m log(m/V)로 추정할 수 있습니다. 논문은 이것을 Whang 등의 Hit Counting과 같은 성질이라고 소개합니다. e^(−n/m)에서 나온 식이므로 여기의 log는 자연로그입니다(이 글의 해석). 문턱값 (5/2)m은 정리에서 유도된 값이 아니라 시뮬레이션에 근거한 값입니다.

대범위 보정. 해시가 L비트이면 서로 다른 해시값은 2^L개뿐이므로, n이 2^L에 가까워지면 해시 충돌이 늘어납니다. 논문은 이때 E가 추정하는 것은 서로 다른 해시값의 개수이고, 이는 높은 확률로 약 2^L(1 − e^(−n/2^L))이라고 설명합니다. 이 관계를 n에 대해 풀면 n = −2^L log(1 − E/2^L)이 되고, 그림 3은 L = 32로 이 식을 씁니다. 문턱값 1/30·2³²에 대해서는 본문에 별도의 유도가 없습니다. 또 이 글의 계산으로는 E ≥ 2³²이면 log의 인자가 0 이하가 되어 식이 정의되지 않는데, 그림 3은 애초에 10⁹까지를 대상으로 치수를 정한 프로그램입니다.

그림 3의 의사코드는 세 조건을 각각 별도의 if로 적어 두었습니다. 글자 그대로 순서대로 실행하면 E ≤ (5/2)m일 때 두 번째 if가 소범위 보정 결과를 E로 덮어씁니다. 주석이 “small range”, “intermediate range”, “large range”로 범위를 나누고 있으므로, 이 글은 세 조건을 서로 배타적인 범위로 읽었습니다. 아래 코드도 그렇게 옮겼습니다.

메모리와 오차를 숫자로

다음 표는 이 글이 계산한 것입니다. 오차 열은 1.04/√m이고, 메모리 열은 레지스터당 5비트로 계산한 레지스터 배열 크기입니다. m = 16처럼 m이 작으면 정리 1의 β₁₆ = 1.106을 쓴 1.106/√16 ≈ 27.7%가 더 정확합니다.

b m 1.04/√m 레지스터 배열 (5비트 × m)
4 16 26% 10바이트
8 256 6.5% 160바이트
10 1,024 3.25% 640바이트
11 2,048 약 2.3% 1,280바이트
16 65,536 약 0.41% 40,960바이트

논문 초록과 서론은 m = 2048, 32비트 해시, 5비트 레지스터로 10⁹이 넘는 카디널리티를 전형적인 정확도 2%, 메모리 1.5kB로 추정할 수 있다고 적습니다. 위 계산으로는 1.04/√2048 ≈ 2.3%, 레지스터 배열은 1,280바이트이므로, 논문의 2%와 1.5kB는 반올림한 대략의 숫자로 읽는 것이 맞습니다. 그림 4의 시뮬레이션 그래프도 예측 표준오차를 m = 256에서 ±7%, m = 1024에서 ±3%, m = 65536에서 ±0.5%로 표시합니다.

논문 4절은 제목의 “near-optimal”을 두 가지 사실로 설명합니다. 하나는 N까지의 ε-근사 개수를 유지하려면 Ω(log log N)비트가 필요하다는 하한입니다. 서론에 따르면 HyperLogLog의 레지스터 하나는 log₂log₂N + O(1)비트 이하이므로 이 하한과 같은 차수입니다. 다른 하나는 논문이 인용한 Chassaing–Gérin의 결과로, 순서통계량 기반 알고리즘의 정확도 하한이 1/√m에 가깝다는 것입니다. 논문은 이를 근거로 HyperLogLog가 그 하한에서 약 4% 떨어져 있다고 말합니다. 이 하한 결과 자체는 이 글에서 원문을 확인하지 않았으므로, 논문이 그렇게 인용한다는 데까지만 적습니다.

Python으로 옮긴 예

아래 코드는 그림 3을 Python으로 옮긴 것입니다. 논문과 다른 점은 이렇습니다. 레지스터 인덱스는 0부터 세므로 j에 1을 더하지 않습니다. 세 범위는 앞 절에서 설명한 대로 서로 배타적으로 처리합니다. 해시로는 SHA-1 다이제스트의 앞 32비트를 씁니다. 논문은 해시 가정을 뒷받침하는 실용적인 방법의 예로 sha1을 언급하는데, 이 코드는 그 언급을 따른 예시일 뿐 속도나 다른 성질에 대해서는 아무것도 주장하지 않습니다.

import hashlib
import math


def hll_estimate(items, b=11):
    """그림 3의 HyperLogLog 프로그램을 옮긴 것. items는 문자열의 반복 가능 객체."""
    if not 4 <= b <= 16:
        raise ValueError("그림 3은 b가 4..16인 경우를 다룬다")
    m = 1 << b
    alpha = {16: 0.673, 32: 0.697, 64: 0.709}.get(m, 0.7213 / (1 + 1.079 / m))
    registers = [0] * m  # 그림 3: −∞가 아니라 0으로 초기화
    for v in items:
        x = int.from_bytes(hashlib.sha1(v.encode("utf-8")).digest()[:4], "big")  # 32비트
        j = x >> (32 - b)                 # 앞 b비트. 논문의 j는 이 값 + 1
        w = x & ((1 << (32 - b)) - 1)     # 나머지 32 − b비트
        rho = (32 - b) - w.bit_length() + 1  # 가장 왼쪽 1의 위치. w = 0이면 32 − b + 1
        registers[j] = max(registers[j], rho)
    e = alpha * m * m / sum(2.0 ** -r for r in registers)  # 원시 추정값
    if e <= 2.5 * m:
        v = registers.count(0)
        return m * math.log(m / v) if v else e  # 소범위 보정
    if e <= 2**32 / 30:
        return e  # 중간 범위, 보정 없음
    return -(2**32) * math.log(1 - e / 2**32)  # 대범위 보정


print(round(hll_estimate(f"item-{i}" for i in range(100_000))))

출력되는 값은 해시와 입력에 따라 정해지는 하나의 추정값이며, 100,000과 같다는 보장은 없습니다. b = 11이면 정리 1이 말하는 전형적인 상대 오차는 약 2.3%이지만, 앞에서 본 것처럼 이것은 개별 실행의 한계가 아닙니다. 또 이 코드는 E ≥ 2³²인 경우를 따로 처리하지 않습니다.

추정값에서 끌어내면 안 되는 것

  • 결정적 보장으로 읽기. 1.04/√m은 표준오차입니다. 논문은 추정값이 근사적으로 정규분포를 따를 것으로 보고(그림 4의 시뮬레이션이 이를 뒷받침한다고 적습니다), σ ≈ 1.04/√m일 때 추정값이 실제 값의 σ, 2σ, 3σ 이내에 들어오는 비율을 각각 65%, 95%, 99%로 예상합니다. 거꾸로 말하면 3σ를 넘는 경우도 예상 범위 안에 있습니다.
  • 작은 n에도 정리 1이 적용된다고 보기. 정리 1은 n → ∞의 점근 결과입니다. 작은 카디널리티를 다루는 소범위 보정과 문턱값 (5/2)m은 시뮬레이션에 근거합니다.
  • 해시 가정 없이 숫자를 쓰기. 분석은 해시 비트가 독립이고 균등하다는 가정 위에 있습니다. 논문은 해시 함수가 이 가정을 만족하지 않는 경우나, 해시를 알고 일부러 고른 입력에 대해서는 분석하지 않습니다.
  • 해시 공간을 넘어서 쓰기. 32비트 해시로 n이 2³²에 가까워지면 E는 서로 다른 해시값의 개수를 추정하게 되고, 그림 3은 이를 대범위 보정으로 다룹니다. 그림 3의 치수와 상수는 10⁹까지, m = 2⁴…2¹⁶을 위한 것입니다.
  • 실험 결과를 일반화하기. 논문의 실험은 2007년의 것입니다. 그림 4는 무작위 균등 데이터로 만든 이상적인 다중집합에 대한 시뮬레이션이고(히스토그램은 500회, 250회), 그림 5는 m = 2048로 한 책의 PostScript 소스를 8진 덤프한 데이터를 약 40,000줄씩 나눈 458개 조각에서 얻은 결과입니다. 실행 시간이 해시 계산에 좌우되어 wc -l로 단순히 훑는 것보다 3~4배 느린 정도라는 서술도 그 시점의 관찰입니다.
  • 지원하지 않는 연산을 기대하기. 논문은 원본을 임의로 나눈 부분 파일들의 레지스터를 성분별 max로 합치면 된다고 적습니다(각주 4). 반면 그림 2와 그림 3에는 레지스터를 max로 키우는 연산만 있으므로, 원소를 삭제하는 연산은 없습니다(이 글의 해석). 교집합 크기를 구하는 방법도 이 논문은 다루지 않습니다.
  • 블룸 필터와 같은 문제로 보기. 블룸 필터는 어떤 값을 넣은 적이 있는지 묻는 소속 질의를 위한 구조이고, 고유값 개수를 추정하는 HyperLogLog와는 답하는 질문이 다릅니다.
  • 후속 변형의 성질을 이 논문에 기대기. HyperLogLog++를 비롯한 이후의 변형과 데이터베이스별 구현은 이 글의 범위 밖이며, 여기서 정리한 상수와 문턱값은 2007년 논문의 그림 3에 대한 것입니다.

출처

두 페이지 모두 2026-10-20에 확인했습니다. 후속 변형이나 데이터베이스 구현에 관한 자료는 사용하지 않았습니다.

블로그로