国旗行列のランクと特異値分解

前回、モンドリアン行列のランクを調べましたが、同じようなルールで各国の国旗を調べてみましょう。ランクとは、線形独立な最大な行数です。旗の例では、足し算や色の交換だけで重複を省いた行のパターン数です。

ちなみに、正確に線形代数のランクを測るための、行列の変換規則を書いておきます。以下の行基本を行なっても行列のランクは変化しません。（）の中に旗で考えるときのわかりやすい解釈も入れましたので、この変換で置き換えてから、最後に残るパターンの数を数え、自分でランクを計算してみてください。

行基本変形

ある行を他の行から引き算する（同じパターンの行があれば、それは消せる）。
ある行のα倍を他の行から引き算する（同じパターンで色違いの行があれば、それは消せる）。
ある行とある行を交換する（近いパターンの行は近くに移動させて考えられる）。
ある行を非零のα倍する（ある行の色を変える。白には変えられない）。

さらに、これらの変形を行でなく列にしても大丈夫。行と列を混ぜて変形してもランクは不変です。

Fracne, Sweden, Greece

France はすべてが同じ行の繰り返しでパターンが1つ。Sweden の行のパターンは2つですね。Greeceの行のパターンは少し難しいですが、4つありそうです。しかし、2行目+3行目=1行目なので、ランク3です。あるいは、列のパターン数を数えた方が分かりやすく、明らかに3パターンです。

Japan

さて、問題は日本です。このような行列の問題では「斜めのライン」がランクに大きく効きます。行基本変形で消えにくいのです。縦や横にぴったり沿ったライン（水平・垂直）は他の行や列と同じパターンになるので、ランクの増加に効きにくいのです。これに対して、三角行列や単位行列を思い出して下さい。斜めのラインがあるために、ランクは最大（斜め線のサイズ）になります。

下図を見て下さい。縦横に整列した正方形は最小ランクで1、45°傾いた正方形＝ダイヤモンドは最大ランクで$n$。と大きな開きがあり、円はその間の数になります。

Prof. Gilbert Strang と Prof. Alex Townsend は、この円のランクを $0.5858n$ 程度と計算しています。動画『なぜ特異値は急激に減少するのか』および書籍『ストラング：教養の線形代数』には、日本の国旗のランクを求める問題があります。解答は動画中で分かりやすく解説されているほか、 Solution Manual の 7.2 節の練習問題 5（p.82）にあります。

ここで、解答のあらすじを。

まず、下のように円から正方形をくり抜きます。この正方形は上に取り出した細長い長方形の rank1 行列を円から引いたと考えます。この細長い行は、円の45°位置と135°位置を結ぶ直線を取り出したもので、ランクは1で、 $n$ が大きくなると無視できます。
次に4つの弓形が円の周囲に残ります。上の1/4は、行数として $n – 1/\sqrt{2}n$ あり、傾斜の傾きが最大45°です。この部分は厚みの行数分だけのランク数と考えられます。真下の鏡像の1/4はこのコピーですので、ランクを増やしません。

次に左の縦の弓形です（右の弓形はこのコピーで無視できます）。この部分がすべて線形独立な（別パターンの）行と考えれば $1/\sqrt{2}n$ なのですが、実際には中央に近づくにつれ垂直に切り立っていき、同じ行とみなせるようになります。そして、この縦の弓形のランクは、先に計算した横の弓形のランク $n – 1/\sqrt{2}n$ と同じと考えられます。なぜ？
横の弓形では行基本変形で考えていましたが、この縦の弓形では列基本変形で考えるのです。列で考えるとその傾斜が垂直に近くなることがなく、すべて独立な列と考えられます。そして行基本変形と列基本変形は同時に行なってもランクを変えないので、1つの弓形の2倍のランクが全体のランクとなります。
よって、結論として、日本の国旗のランクは、$n$ が大きい場合、小さな rank 1を無視して、

$$
2(n – 1/\sqrt{2}n) = (2 – \sqrt{2})n \approx 0.5858n
$$

となります。では数値実験してみましょう。

日本国旗のランク計算

キャンバスサイズを 200×300 ピクセル、円の半径を $n=100$ と設定して計算します。ランクの計算には特異値分解を用います（ランクは非零特異値の数です）。

特異値分解(compact) は、行列 $A \in \bR^{m \times n}$ を $r$ をランクとして、2つの直交行列 $U \in \bR^{m \times r}, V \in \bR^{n \times r}$ と対角行列 $\Sigma \in \bR^{r \times r}$ の積に分解します。

$$
A = U \Sigma V\transp = \sigma_1 \bu_1 \bv_1\transp + \sigma_2 \bu_2 \bv_2\transp + \cdots \sigma_r \bu_r \bv_r\transp
$$

例えば $r=2$ の場合、こんな風に分解され、積と和の両方で表現できます。

$\Sigma$ の対角には、特異値（$\sigma_1 \cdots \sigma_n$ ）が単調減少で並びます。そして、非零の $\sigma$ の数がランクになります。特異値分解を使って特異値を取り出すことで、各ランクの「強さ」が正の実数として取り出せるのです。$\sigma_1$ は最強のシグナル、 $\sigma_2$ はその次、と徐々にシグナルが弱くなり、最後 $\sigma_{r+1}$ で完全に $0$ になります。

$$
\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r \ge \sigma_{r+1} = \cdots = 0
$$

特異値分解はその安定したアルゴリズムによって、特異値（実数）の列を作り、その非零の数（正の整数）がランクです。ランクは整数で非零 $\sigma$ 特異値の数ですが、零と非零の境界（ $\epsilon$ ）をどこにおくかで、現実的にその数は変わります。これには、(1)そもそも行列は整数の行数と列数を持つので円を精密には表現できないこと（実数の整数による近似）、(2)コンピュータ上で実数を浮動小数点で扱う場合避けられない微少量（機械 $\epsilon$ ： $1 + \epsilon = 1$ となる値）の複合問題です。

	32bit (float)	64bit (double)
機械 $\epsilon$	$\approx 10^{-7}$	$\approx 10^{−16}$

今回設定した半径 $n=100$ をこの公式に当てはめると、次のようになります。

$$(2-\sqrt{2}) \times 100 \approx 0.58578 \times 100 = 58.58$$

つまり、理論上は「58〜60個あたりの特異値を用意すれば完璧な円が構成でき、それ以降の特異値は、計算機上足しても変化しない量に落ち込むはずだ」という予測が立ちます。

数値実験：特異値の推移と復元画像

この理論値が実際のデジタル画像でどの程度成り立つのか、Pythonを用いて数値実験を行いました。

① 特異値のグラフ（ゼロへの急落）

計算された特異値を対数グラフにプロットしたところ、驚くべき結果です！特異値は58番目付近（特異値は大きい方から順に0から番号付けされている）までは意味のある数値を保っていましたが、赤い点線で示した 理論値 58.58 のラインを越えた直後、$10^{-13}$ 以下へと急落しました。

プログラムで $10^{-13}$ より大きい特異値をカウントすると、その数は60個。評価式の予測と極めて高い精度で一致しました。

この60という数を、機械ランク( Machine Rank )や実行ランク( Effective Rank )、数値ランク( Numerical Rank ) と呼びます。

Pythonの mumpy が用いる浮動小数点（np.float64）はC言語のdoubleに対応し、機械 $\epsilon \approx 10^{-16}$ 。この値は1.0 付近の足しても変化しない最小値と定義されています。実際の計算では足す相手の数によって足しても変化しない最小値は変わります。この実験では $10^{-14}$ あたりが最小値に見えます。

② 復元画像の推移（n=1, 2, 3…）

$$
A = U \Sigma V\transp = \sigma_1 \bu_1 \bv_1\transp + \sigma_2 \bu_2 \bv_2\transp + \cdots
$$

という減少する特異値列から得られる行列を、ある番号で打ち切ることで「近似」することができます（これが最良近似になっていることを示す Eckert-Young の定理があります）。含める特異値の数を $n=1, 2, 3…$ と徐々に増やしていくと、画像はどのように変化するのでしょうか。

Rank 1: 円の面影はなく、中央にぼんやりとした十字型の影が現れます。
Rank 2: 四角形のような形が浮かび上がります。情報量（エネルギー）の大部分はここで既に確保されています。
Rank 3〜4: ひし形や八角形のような多角形に変化し、徐々に「丸み」を帯びてきます。
Rank 10: パッと見は「円」ですが、境界部分を無理やり直交成分で作っているため、波紋状のノイズ（リンギング）が見えます。
Rank 50: ほぼ完璧な円に見えます。
Rank 60: 特異値がゼロに落ち込んだ後の画像です。これ以降はいくら特異値を足しても、画像は1ピクセルたりとも変化しません（機械 $\epsilon $ 付近に入った）。これ以上の円はピクセルでは表現できません。

まとめ

画像圧縮の観点では「最初の数個（ランク3〜4）の特異値」だけで、日の丸の大半の情報は表現できてしまいます。しかし、ピクセルの世界で「完璧な丸」の境界線を研ぎ澄ますためには、そのあとも細かな特異値がロングテールで使われ続けます。

そして最も美しいのは、完全な円を構成するための「最後のピース」の場所です。円の半径の 約58.58% の数まで特異値を足し合わせると、デジタル上で表現可能な「完璧な円」に到達し、それ以上の情報は一切不要になります。

$(2-\sqrt{2})n \approx 0.5858$ という数は、線形代数（特異値）、数値計算（浮動小数の限界）と画像処理（整数個のピクセル）が交差する点です。

シンプルな国旗の図形の中に、このような深い数学の結びつきが隠されていることに感動しました。

追記

この記事の元になった、モンドリアンの絵画のようなブロック状行列を、特異値分解でなくガウスの消去法や $CR$ 分解によってランクを求める記事については、「モンドリアン行列のランク」を参照してください。

（※Edit 2026/4/4 : この記事を Prof. Gilbert Strang に送ったところ、Prof. Alex Townsend ご本人にも記事を紹介して頂き、さっそくミスを１つ指摘して頂けました（ダイヤモンドのランクは $n$ ）。また、この契機に、WordPress に MathJax Plugin を入れていつも使っている TeX コマンドをカスタマイズし、また、TranslatePress Pluigin を入れて英語への自動変換を右下のフローティングボタンから可能にして、英語圏の方に参照しやすくしました。）

付録：検証に使用したPythonコード

本記事のグラフと画像は、以下のコードで生成しています。

特異値分解による特異値の減衰

import numpy as np
import matplotlib.pyplot as plt
# 設定
height = 200
width = 300
n = 100 # 半径
y, x = np.ogrid[:height, :width]
cy, cx = height / 2, width / 2
mask = (x - cx)**2 + (y - cy)**2 <= n**2
# 行列の作成（円自体=1、背景=0）
A = np.zeros((height, width))
A[mask] = 1
# SVDの計算
U, S, Vt = np.linalg.svd(A, full_matrices=False)
# 予測される非零特異値の数
predicted_rank = (2 - np.sqrt(2)) * n
plt.figure(figsize=(10, 6))
# 特異値のプロット
plt.plot(S, marker='.', linestyle='-', color='b', label='Singular Values (Circle)')
# (2-√2)n のライン
plt.axvline(x=predicted_rank, color='r', linestyle='--', linewidth=2, 
            label=f'(2-√2)n ≈ {predicted_rank:.2f} (n={n})')
# ゼロの閾値ライン
plt.axhline(y=1e-10, color='g', linestyle=':', label='Zero Threshold (1e-10)')
plt.yscale('log')
plt.title(f'Singular Values dropping to Zero (Circle Radius n={n}, Flag Size 200x300)')
plt.xlabel('Index')
plt.ylabel('Singular Value (Log Scale)')
plt.xlim(0, 150) # n=100なので少し広めに表示
plt.ylim(1e-16, 1e4)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig('circle_rank_n100.png')
# 実際の非零特異値の数をカウント
actual_rank = np.sum(S > 1e-10)
print(f"Radius (n): {n}")
print(f"Predicted rank (2-√2)n: {predicted_rank:.2f}")
print(f"Actual non-zero singular values: {actual_rank}")

国旗画像の復元

import numpy as np
import matplotlib.pyplot as plt
# キャンバスサイズと円の半径の設定
height = 200
width = 300
r = 50
# RGBの画像配列を作成（背景は白: 1.0, 1.0, 1.0）
img = np.ones((height, width, 3))
# 日の丸の赤色（#BC002D -> RGB: 188/255, 0, 45/255）
red_color = np.array([188/255, 0, 45/255])
# 円形マスクの作成
y, x = np.ogrid[:height, :width]
cy, cx = height / 2, width / 2
mask = (x - cx)**2 + (y - cy)**2 <= r**2
# マスク部分を赤色に
img[mask] = red_color
# RGB各チャンネルごとに特異値分解 (SVD) を実行
U_r, S_r, Vt_r = np.linalg.svd(img[:,:,0], full_matrices=False)
U_g, S_g, Vt_g = np.linalg.svd(img[:,:,1], full_matrices=False)
U_b, S_b, Vt_b = np.linalg.svd(img[:,:,2], full_matrices=False)
# 指定されたランク
ranks = [1, 2, 3, 4, 5, 10, 20, 30, 40, 50]
fig, axes = plt.subplots(2, 5, figsize=(15, 6))
axes = axes.flatten()
for i, k in enumerate(ranks):
    # 各チャンネルごとにランクkで画像を再構成
    R_k = U_r[:, :k] @ np.diag(S_r[:k]) @ Vt_r[:k, :]
    G_k = U_g[:, :k] @ np.diag(S_g[:k]) @ Vt_g[:k, :]
    B_k = U_b[:, :k] @ np.diag(S_b[:k]) @ Vt_b[:k, :]
    
    # チャンネルを結合してRGB画像に（値が0未満や1を超過する場合はクリップ）
    img_k = np.stack((R_k, G_k, B_k), axis=2)
    img_k = np.clip(img_k, 0, 1)
    
    # プロット
    axes[i].imshow(img_k)
    axes[i].set_title(f'Rank n={k}')
    axes[i].axis('off')
    
    # わかりやすいように薄い枠線を追加
    rect = plt.Rectangle((0, 0), width-1, height-1, fill=False, edgecolor='gray', linewidth=0.5)
    axes[i].add_patch(rect)
plt.tight_layout()
plt.savefig('svd_color_flag_reconstruction.png', dpi=150)

an Agile Way

Software is made of conversations

国旗行列のランクと特異値分解

行基本変形

Fracne, Sweden, Greece

Japan

日本国旗のランク計算