9.1 浮點數的真相#

本節目標#

• 了解 IEEE 754 浮點數的運作方式，以及為什麼它「天生不精確」
• 親眼驗證 \(0.1 + 0.2 \neq 0.3\) ，並理解原因
• 看到浮點誤差對線性代數的實際衝擊： \(AA^{-1} \neq I\)
• 學會用容差（tolerance）正確地比較浮點數
• 做一個捨入誤差累積實驗，誠實面對「電腦算的也可能錯」

電腦裡的數字是有限的#

數學上的實數有無窮多個，但電腦的記憶體是有限的。 float64（NumPy 的預設型別）用 64 個位元 表示一個數： 1 個位元的正負號、11 個位元的指數、52 個位元的小數部分（尾數）。

這代表：連續兩個可表示的浮點數之間有縫隙。 在 1.0 附近，這個縫隙大約是 \(2^{-52} \approx 2.2 \times 10^{-16}\) ， 這個量級有個名字叫「機器 epsilon」：

import numpy as np

print(np.finfo(float).eps)   # 2.220446049250313e-16
print(np.spacing(1.0))       # 1.0 旁邊下一個浮點數的距離，同一件事

數學上實數是連續的；電腦上的數字是一格一格的刻度。 每一次運算，結果都會被「吸附」到最近的刻度上，這就是捨入誤差的來源。

現場驗證： \(0.1 + 0.2 \neq 0.3\)#

這是浮點數最著名的「案發現場」。0.1 在二進制下是無限循環小數 （就像 1/3 在十進制下寫不完），所以存進電腦的那一刻就已經被截斷了：

print(0.1 + 0.2)          # 0.30000000000000004
print(0.1 + 0.2 == 0.3)   # False —— 這不是 bug，是浮點數的本性

# 看得更清楚：用 17 位有效數字印出內部儲存值
print(f"{0.1:.17f}")      # 0.10000000000000001
print(f"{0.2:.17f}")      # 0.20000000000000001
print(f"{0.1 + 0.2:.17f}")  # 0.30000000000000004

0.1 和 0.2 各自帶了一點點誤差，相加之後誤差沒有抵銷、反而現形， 所以結果比 0.3 多了一點點。這不是 Python 的錯，所有遵循 IEEE 754 的語言都一樣（包括 C、Java、Excel）。

浮點誤差在線性代數的衝擊#

數學課本告訴我們 \(AA^{-1} = I\) 。在電腦上呢？我們來實測：

import numpy as np

rng = np.random.default_rng(42)
A = rng.standard_normal((5, 5))
Ainv = np.linalg.inv(A)

E = A @ Ainv
print(np.abs(E - np.eye(5)).max())
# 典型輸出：1.1e-15 左右 —— 不是 0！

# 更慘烈的例子：接近奇異的矩陣
A = rng.standard_normal((8, 8))
A[7] = A[0] + A[1] + 1e-10 * rng.standard_normal(8)  # 幾乎線性相依
E = np.abs(A @ np.linalg.inv(A) - np.eye(8))
print(E.max())   # 誤差可能暴增好幾個數量級

第一次實驗的誤差約 \(10^{-15}\) ，聽起來很小，但注意兩件事：

1. 它不是零。如果你的程式碼寫 if (A @ Ainv == I)，永遠不會成立。
2. 當矩陣接近奇異時（第二次實驗），誤差會被放大到完全不可接受—— 這是下一節「條件數」要處理的問題。

比較浮點數一定要用容差#

既然相等比較不可靠，正確的做法是「夠接近就算相等」：

a, b = 0.1 + 0.2, 0.3

print(np.isclose(a, b))       # True —— 絕對容差約 1e-08
print(np.allclose([0.1+0.2], [0.3]))  # 陣列版：全部接近才回傳 True

# 自己指定容差（嚴格的科學計算常需要更緊的標準）
print(np.allclose(A @ Ainv, np.eye(8), atol=1e-8))

# 相對容差的觀念：比較 1e9 附近的數和比較 1e-9 附近的數，
# 合理的容差應該差九個數量級。np.isclose 預設混用兩者：
# |a - b| <= atol + rtol * |b|
print(np.isclose(1e9, 1e9 + 1))        # True（相對誤差 1e-9，很小）
print(np.isclose(0.0, 1e-9))           # False（b 是 0 時只剩絕對容差）

準則：任何時候想對浮點結果寫 ==，停下來，改用 np.allclose。

捨入誤差累積實驗#

單次誤差 \(10^{-16}\) 很小，但累積起來呢？我們把 0.1 連加一萬次：

import numpy as np

# 方法一：迴圈逐次累加 —— 誤差會滾雪球
total = 0.0
for _ in range(10000):
    total += 0.1
print(repr(total))    # 999.9999999999062，不是 1000！

# 方法二：np.sum —— 成對求和（pairwise summation），誤差小得多
arr = np.full(10000, 0.1)
print(repr(arr.sum()))    # 1000.0000000000001

# 方法三：math.fsum —— 精確求和，拿來當「標準答案」
import math
print(repr(math.fsum([0.1] * 10000)))  # 1000.0000000000001（真正的和）

print(abs(total - 1000))    # 約 9.4e-11
print(abs(arr.sum() - 1000))  # 約 1e-12

# 更極端：加一千萬次，迴圈版的誤差會達到 1e-6 量級
total = 0.0
for _ in range(10_000_000):
    total += 0.1
print(abs(total - 1_000_000))   # 你自己看結果，會嚇一跳

有趣的事實：np.sum 用「成對求和」策略——先把相鄰的數兩兩相加再往上合併， 讓誤差成長速度從 \(O(n)\) 降到 \(O(\log n)\) 。連求和的順序都會影響答案， 這就是浮點世界的日常。

延伸陷阱：矩陣乘法中不同的演算法（分塊、平行化）加法順序不同， 所以同一個矩陣乘積在不同 BLAS 函式庫上結果可能在小數第 12 位以後不同。 這不是誰算錯，是浮點數的必然。

誠實面對：「電腦算的也可能錯」#

很多人以為「電腦算的就是對的」。這一章的第一課要打破這個迷信：

好消息是：這些誤差通常在 \(10^{-16}\) 量級，對絕大多數應用無害。 壞消息是：當你不知道它存在，它就會在最關鍵的時候咬你—— 例如比對結果是否收斂、判斷某個值是否為零、長時間累積的模擬。

三個務實的習慣：

1. 比較一律用容差：np.allclose 而非 ==。
2. 對誤差有感：知道你的計算大約會累積多少誤差（本章教 you 估算）。
3. 需要更高精度時有工具：np.float128（延伸精度）、math.fsum（精確求和）、 或重新設計演算法（下一節的條件數就是「重新設計」的起點）。

# 最後一個小實驗：加法結合律在浮點數下失效
a, b, c = 1e16, 1.0, -1e16
print((a + b) + c)   # 0.0   —— b 被吞掉了
print(a + (b + c))   # 1.0   —— 正確答案

\(10^{16}\) 加上 1，在 float64 裡根本沒地方放那個 1（刻度間距大於 2）， 於是它被「捨」掉了。大數吃小數是浮點誤差最常見的形態之一。

小結與練習#

• float64 只有 64 個位元，數字之間有縫隙，每次運算都產生約 \(10^{-16}\) 的捨入誤差
• \(0.1 + 0.2 \neq 0.3\) ：二進制表示 + 誤差現形，所有 IEEE 754 系統皆然
• \(AA^{-1} \neq I\) ：線性代數的結果永遠帶誤差，接近奇異時誤差被放大
• 浮點比較必用容差：np.isclose / np.allclose
• 誤差會累積：求和順序影響結果，np.sum 的成對策略比樸素迴圈好

練習：

1. 求 np.finfo(float).max，並驗證 max + 1 == max 成立（大數吃小數）。
2. 用迴圈把 1e-8 連加 \(10^8\) 次，再與 np.sum 的結果比較，誤差差幾個數量級？
3. 寫一個函式 my_allclose(a, b, rtol=1e-5, atol=1e-8)， 依照 \(\lvert a-b \rvert \leq \mathrm{atol} + \mathrm{rtol} \cdot \lvert b \rvert\) 實作，並與 np.allclose 對照測試。
4. 隨機生成 100 個矩陣，統計 \(\lVert AA^{-1} - I \rVert_\infty\) 的最大值，感受「典型誤差量級」。