HyperLogLog：調和平均で異なる要素の数を推定する仕組みと1.04/√m誤差の前提
2007年のFlajoletらの原論文に基づき、HyperLogLogがレジスタ番号jとρで異なる要素の数を推定する方法、調和平均を使う理由、小範囲・大範囲の補正、1.04/√mの誤差が前提とする条件を整理します。
重複を除いた要素の数をカーディナリティ（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回走査し、カーディナリティを推定するアルゴリズムを示し、その相対誤差（標準誤差）がおおむね1.04/√mであることを解析しています。この記事は、FlajoletのINRIAのページで公開されているこの論文のPDFに基づいて、レジスタ番号jとρの求め方、調和平均を使う理由、1.04/√mという数字の前提、そして推定値から引き出してはいけない結論を扱います。その後に登場したHyperLogLog++などの変形や、特定のデータベースの実装は扱いません。
一つの観測量：ρ
ハッシュ値をビット列とみなし、ρ(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本の部分ストリームのうち一つを選び、部分ストリームごとに観測量を別々に保持する方法です。ハッシュの計算は要素1個につき1回で済みます。
ハッシュ関数そのものについては、ハッシュ関数の記事を参照してください。
レジスタ番号jとρ(w)：1はどこで足されるのか
論文の図2はアルゴリズムを次のように定めています。
- m = 2^b（bは正の整数）とし、レジスタM[1], …, M[m]を−∞で初期化します。
- 各要素vについてx = h(v)を計算します。
- j = 1 + ⟨x₁x₂…x_b⟩₂。xの先頭bビットを2進数として読んだ値に1を足したものがレジスタ番号です。
- w = x_{b+1}x_{b+2}…。先頭bビットを取り除いた残りからρ(w)を求め、M[j] = max(M[j], ρ(w))で更新します。
- すべての要素を読み終えたら、次の節の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ビットのハッシュならレジスタ1個あたり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で要素数を固定した場合に戻す、という順に進みます。分散の解析も同じ方針で、概略だけが示されています。
生の推定値と二つの補正
図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、メモリの列はレジスタ1個あたり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のレジスタ1個はlog₂log₂N + O(1)ビット以下なので、この下限と同じオーダーです。もう一つは、論文が引用するChassaingとGérinの結果で、順序統計量に基づく幅広いアルゴリズムでは達成できる精度の下限が1/√mに近いというものです。論文はこれを根拠に、HyperLogLogがその下限から約4%離れていると述べています。この記事ではChassaingとGérinの論文そのものは確認していないので、HyperLogLogの論文がそのように引用している、というところまでにとどめます。
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%ですが、前に見たとおり、これは1回ごとの実行の限界ではありません。また、このコードは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についてのものです。
出典
- Philippe Flajolet, Éric Fusy, Olivier Gandouet, Frédéric Meunier, “HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm,” AofA 2007, DMTCS Proceedings vol. AH, pp. 127–146. PDF: https://algo.inria.fr/flajolet/Publications/FlFuGaMe07.pdf
- Philippe Flajoletの論文リスト（INRIA）、上記論文の2007年の項目: https://algo.inria.fr/flajolet/Publications/publist.html
いずれのページも2026-10-20に確認しました。後続の変形やデータベース実装に関する資料は使っていません。