平行演算法考量與數值穩定性 (Algorithm Considerations and Numerical Stability)

重點總覽 (Overview)

加總順序與平行 Reduction 的準確度 (Summation Order in Parallel Reduction)

矩陣乘法的 dot product、reduction 等都需要把大量數值加總。理想上順序無關(加法符合結合律),但在有限精度下,順序會改變最終準確度。

範例:5-bit 格式 (1-bit S, 2-bit E, 2-bit M)

加總 1.00B×2⁰ + 1.00B×2⁰ + 1.00B×2⁻² + 1.00B×2⁻²(即 1 + 1 + 0.25 + 0.25,精確值 = 2.5 = 1.01B×2¹)。

循序 (sequential, 左到右):

1.00×2⁰ + 1.00×2⁰              = 1.00×2¹            (= 2)
1.00×2¹ + 1.00×2⁻²            → 1.00×2¹            (0.25 對齊右移後消失)
1.00×2¹ + 1.00×2⁻²            → 1.00×2¹  = 2       (再次消失,誤差 0.5)

平行 (parallel reduction tree, 見第 10 章):

1.00×2⁰ + 1.00×2⁰  = 1.00×2¹           (= 2)
1.00×2⁻² + 1.00×2⁻² = 1.00×2⁻¹         (= 0.5)
1.00×2¹ + 1.00×2⁻¹ = 1.01×2¹  = 2.5    (精確!)

ASCII:兩種加總路徑

循序 (left-to-right)                     平行 (reduction tree)
  1   1  0.25 0.25                          1   1   0.25 0.25
   \ /                                       \ /     \ /
    2  +0.25 → 2  (0.25 消失)                  2       0.5
         \                                      \     /
          2 +0.25 → 2  (0.25 消失)               2 + 0.5
              = 2   (誤差 0.5)                     = 2.5  (精確)

提升準確度的技巧

線性解算器與數值穩定性 (Linear Solvers and Numerical Stability)

線性方程組解算器中,不同的輸入數值可能需要不同的運算順序才能求解 — 比 reduction 的順序問題更嚴重。

Gaussian Elimination 流程

以下系統(解為 X=1, Y=2, Z=3):

3X + 5Y + 2Z = 19
2X + 3Y +  Z = 11
 X + 2Y + 2Z = 11

矩陣 (augmented) 視角的消元 + 回代:

[ 3 5 2 | 19 ]   step1: 各 row 除以 lead 係數,使 X 係數一致
[ 2 3 1 | 11 ]   step2: 減去 top row → 消去 X
[ 1 2 2 | 11 ]   重複對 Y ... → upper-triangular（上三角）
                 → 對角線化為 identity matrix
                 → backward substitution: 由最後一式往回解 Z→Y→X

• 消元後成為 triangular matrix;最終對角線全為 1 即 identity matrix(任意矩陣乘 identity 等於自身)。
• 對矩陣做 Gaussian elimination 等價於乘上其反矩陣 (inverse)。
• 第二階段 backward substitution:由最後一式 (Z=3) 往前代回求 Y=2、X=1。
• 複雜度約 O(n³)。

平行化策略 (CUDA)

每個 thread 負責矩陣的一個 row;若系統可放進 shared memory,用一個 thread block 處理。所有 thread 同步迭代各 step:

// 一 thread 一 row;augmented matrix A 為 n×(n+1),放在 shared memory
__global__ void gaussianElim(float* A, int n) {
    int row = threadIdx.x;              // 每個 thread 擁有一個 row
    for (int k = 0; k < n; ++k) {       // 逐欄 k 消元
        // ── (pivoting 步驟會插在這裡:找 |A[i][k]| 最大者並 swap row) ──
        if (row == k) {                 // pivot row 正規化:除以 lead 係數
            float piv = A[k*(n+1) + k];
            for (int j = k; j <= n; ++j) A[k*(n+1)+j] /= piv;
        }
        __syncthreads();                // barrier:pivot row 已就緒
        if (row > k) {                  // 下方各 row 減去 pivot row → 消去變數
            float f = A[row*(n+1) + k];
            for (int j = k; j <= n; ++j)
                A[row*(n+1)+j] -= f * A[k*(n+1)+j];
        }
        __syncthreads();                // barrier:更新完成才進下一欄
    }
}

• division step 後 barrier,再做 subtraction step,再 barrier。
• 每完成一欄,負責該 row 的 thread 停止參與(其 row 在回代前已無工作),其餘繼續。
• 變數越多,平行化加速越顯著。

數值不穩定:lead 係數為 0

把第一式換成 5Y + 2Z = 16(X 係數 = 0):

0X + 5Y + 2Z = 16      ← step1 需除以 X 係數 0 → 失敗!
2X + 3Y +  Z = 11
 X + 2Y + 2Z = 11

此系統仍可解 (X=1, Y=2, Z=3),但 naive Gaussian elimination 除以零而失敗 → numerically unstable。

Pivoting:交換 row 求穩定

[ 0 5 2 | 16 ]   lead 係數 = 0 → 無法 divide
[ 2 3 1 | 11 ]   找 |lead| 最大的 row → swap row0 ↔ row1
[ 1 2 2 | 11 ]
        ↓
[ 2 3 1 | 11 ]   lead 係數 = 2 ≠ 0 → 正常進行 Gaussian elimination
[ 0 5 2 | 16 ]
[ 1 2 2 | 11 ]

Pivoting 的平行成本

兩種降低全域檢查成本的途徑:

附錄總結 (Appendix Summary)

考試/面試重點 (Exam / Test Patterns)

Related Notes

• 24-Numerical-Considerations/01-Floating-Point-Data-Representation
• 24-Numerical-Considerations/02-Representable-Numbers-Precision-and-Accuracy
• 10-Reduction/01-Reduction-Fundamentals-and-Simple-Kernel
• 13-Sorting/01-Sorting-Foundations-and-Parallel-Radix-Sort
• 03-Multidimensional-Grids-And-Data/04-Matrix-Multiplication-Kernel
• 19-Parallel-Programming-And-Computational-Thinking/02-Algorithm-Selection-and-Tradeoffs