Flag Matrix Rank and Singular Value Decomposition

My bookIllustrated Linear Algebra: Intuitive Understanding the Strang Way.is a blog series.

Previous,Rank of a Mondrian matrixwe have examined the flags of different countries using the same rules. Rank is the maximum number of linearly independent rows. In the flag example, it is the number of row patterns that eliminate duplicates by just adding and exchanging colors.

By the way, to accurately measure the rank of linear algebra, I will write down the matrix transformation rules. The following row basis does not change the rank of the matrix. (I have also included an easy-to-understand interpretation when thinking in terms of flags in parentheses (). After replacing with this transformation, count the number of patterns that remain at the end and calculate the rank by yourself.

line base transformation

  • Subtract one line from another (if there are lines with the same pattern, they can be eliminated).
  • Subtract alpha times one row from the other rows (if there are different color rows with the same pattern, they can be eliminated).
  • Exchange one row for another (rows with close patterns are considered by moving them closer together).
  • Multiply a row by a non-zero alpha (change the color of a row. Cannot change to white).

Furthermore, you can make these transformations in columns instead of rows. The ranks are invariant even if you mix rows and columns in the transformations.

Fracne, Sweden, Greece

France has all the same lines repeated, with one pattern. Sweden has two line patterns. Greece's line patterns are a bit tricky, but there seem to be four. However, since the 2nd line + 3rd line = 1st line, the rank is 3. Alternatively, counting the column patterns is easier to understand, and there are clearly three patterns.

Japan

Now, the problem is Japan. In problems involving matrices like this, "diagonal lines" significantly affect the rank. They are difficult to eliminate with elementary row operations. Lines that perfectly align horizontally or vertically (horizontal/vertical) follow the same pattern as other rows or columns, making them less effective in increasing the rank. In contrast, recall triangular matrices and identity matrices. Due to the presence of diagonal lines, the rank becomes maximum (equal to the size of the diagonal).

Look at the figure below. Squares aligned vertically and horizontally have a minimum rank of 1, while squares tilted at 45 degrees, or diamonds, have a maximum rank of$n$. There is a large gap between them, with circles falling somewhere in between.

Prof. Gilbert Strang and Prof. Alex Townsend calculate the rank of this circle to be approximately$0.5858n $. The video "Why Do Singular Values Drop So Fast?" and the book "Strang: Linear Algebra for Everyone" include a problem about finding the rank of the Japanese flag. The solution is clearly explained in the video, and Solution Manual It's in Exercise 5 of Section 7.2 (p. 82).

Here's a summary of the solution.

  • まず、下のように円から正方形をくり抜きます。この正方形は上に取り出した細長い長方形の 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$ と同じと考えられます。なぜ?
  • In the case of the horizontal bow shape, we considered row elementary operations, but for this vertical bow shape, we consider column elementary operations. When considering columns, their inclination never becomes close to vertical, and they can all be considered independent columns. Since row and column elementary operations do not change the rank even when performed simultaneously, the rank of the entire matrix is twice the rank of one bow shape.
  • Thus, in conclusion, if the rank of the Japanese flag is $n$ is large, ignoring the small rank 1,

$$
2(n – 1/√2n) = (2 – √2)n ≈ 0.5858n
$$

This is how it will be. Now let's do a numerical experiment.

Japanese flag ranking calculation

キャンバスサイズを 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}$

Applying this formula to the radius $n=100$ set this time, we obtain

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

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

Numerical Experiment: Singular Value Trend and Reconstructed Image

I conducted a numerical experiment using Python to see how well this theoretical value holds up in actual digital images.

(1) Graph of Singular Value (Plummet to Zero)

A surprising result when plotting the calculated singular values on a logarithmic graph! The singular values maintained meaningful values up to around the 58th singular value (singular values are numbered from 0 in descending order), but the red dashed line shows Theoretical value 58.58 のラインを越えた直後、$10^{-13}$ 以下へと急落しました。

プログラムで $10^{-13}$ より大きい特異値をカウントすると、その数は60 pcs.The The results agreed with the predictions of the evaluation equation with extremely high accuracy.

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

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

② Transition of restored images (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: The impression of a circle is gone, and a faint cross-shaped shadow appears in the center.
  • Rank 2: A quadrilateral-like shape emerges. The majority of the information (energy) is already secured here.
  • Rank 3〜4: It changes to polygonal shapes such as rhombus and octagon, and gradually becomes "rounded".
  • Rank 10: At first glance, it looks like a "circle," but because the boundary is forcibly created with orthogonal components, ripple-like noise (ringing) is visible.
  • Rank 50: It looks like a near-perfect circle.
  • Rank 60: 特異値がゼロに落ち込んだ後の画像です。これ以降はいくら特異値を足しても、画像は1ピクセルたりとも変化しません(機械 $\epsilon $ 付近に入った)。これ以上の円はピクセルでは表現できません。

summary

From the perspective of image compression, most of the information in a "red circle" can be represented by just "the first few (ranks 3-4) singular values." However, in the world of pixels, fine singular values continue to be used in a long tail to sharpen the boundary of a "perfect circle."

And most beautiful is the place of the "last piece" to complete the perfect circle. The radius of the circle Approx. 58.58% If you sum singular values up to a certain number, you reach a "perfect circle" that can be represented digitally, and no further information is needed.

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

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

追記

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

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


Appendix: Python code used for verification

The graphs and images in this article are generated using the following code.

Decay of singular values by singular value decomposition
Python
import NumPy as No problem.
import matplotlib.pyplot as plt
# Settings
Height = 200
width = 300
n = 100 # radius
y, x = No problem..Ogrid[:Height, :width]
cy, cx = Height / 2, width / 2
mask = (x - cx)**2 + (y - cy)**2 Less than or equal to n**2
# Matrix Creation (Circle Itself=1, Background=0)
A = No problem..zeros((Height, width))
A[mask] = 1
# SVD calculation
U, S, VT = No problem..Linear Algebra.SVD(A, full_matrices=False)
# 予測される非零特異値の数
predicted_rank = (2 - No problem..square root(2)) * n
plt.figure(figsize=(10, 6))
# Plot of singular values
plt.Plot(S, marker=.;, linestyle=-;, color=b;, label=Singular Values (Circle);)
# line with (2-√2)n
plt.horizontal line(x=predicted_rank, color=R;, linestyle=--;, Line width=2,
label=f'(2-√2)n ≈ {predicted_rank:.2f} (n={n}It's;)
# zero threshold line
plt.axhline(y=1e-10, color=g;, linestyle=:;, label=Zero Threshold (1e-10);)
plt.Y-scale(log;)
plt.Title(Singular Values Dropping to Zero (Circle Radius n={n}, Flag Size 200x300;)
plt.X-axis label(Index;)
plt.ylabel(Singular Value (Log Scale);)
plt.xlim(0, 150) # Since n=100, I'll display it a bit wider.
plt.Ylim(1e-16, 1e4)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.save_fig(circle_rank_n100.png;)
Count the number of actual non-zero singular values #
actual_rank = No problem..sum(S > 1e-10)
print(Radius (n): {n}")
print(Predicted rank (2-√2)n: {predicted_rank:.2f}")
print(Actual non-zero singular values: {actual_rank}")
Flag image restoration
Python
import NumPy as No problem.
import matplotlib.pyplot as plt
# Canvas Size and Circle Radius Settings
Height = 200
width = 300
r = 50
# Create an RGB image array (background is white: 1.0, 1.0, 1.0)
image = No problem..ones((Height, width, 3))
#The red of the Japanese flag (# BC002D -> RGB: 188/255, 0, 45/255)
red color = No problem..array([188/255, 0, 45/255])
# Creation of a circular mask
y, x = No problem..Ogrid[:Height, :width]
cy, cx = Height / 2, width / 2
mask = (x - cx)**2 + (y - cy)**2 Less than or equal to r**2
# The mask part in red
image[mask] = red color
# Perform Singular Value Decomposition (SVD) for each RGB channel
You are, S_r, Vt_r = No problem..Linear Algebra.SVD(image[:,:,0], full_matrices=False)
U_g, Sg, Vt_g = No problem..Linear Algebra.SVD(image[:,:,1], full_matrices=False)
U_b, S_b, VT_b = No problem..Linear Algebra.SVD(image[:,:,2], full_matrices=False)
# Specified Rank
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):
# Reconstruct image with rank k for each channel
R_k = You are[:, :k] @ No problem..diag(S_r[:k]) @ Vt_r[:k, :]
G_k = U_g[:, :k] @ No problem..diag(Sg[:k]) @ Vt_g[:k, :]
B_k = U_b[:, :k] @ No problem..diag(S_b[:k]) @ VT_b[:k, :]
# Combine channels into RGB image (clip if value is less than 0 or greater than 1)
image_k = No problem..Stack((R_k, G_k, B_k), axis=2)
image_k = No problem..clip(image_k, 0, 1)
# Plot
axes[i].imshow(image_k)
axes[i].set_title(rank n={k}'.;)
axes[i].axis(off;)
# Thin border added for clarity
rectangle = plt.Rectangle((0, 0), width-1, Height-1, Fill=False, edgecolor=gray;, Line width=0.5)
axes[i].add_patch(rectangle)
plt.tight_layout()
plt.save_fig(svd_color_flag_reconstruction.png;, Dots Per Inch=150)

Leave a Reply