Kaczmarz 算法

本文的閱讀等級:中級

Kaczmarz 算法是線性方程 A\mathbf{x}=\mathbf{b} 的一種迭代解法,1937年由波蘭數學家科區馬茲 (Stefan Kaczmarz) 所提出[1]。1970年,戈登 ( Richard Gordon)、本德爾 (Robert Bender) 和赫爾曼 (Gabor Herman) 三人重新發現此法,稱之為代數重建技術 (algebraic reconstruction technique),主要應用於電腦斷層掃描的影像重建[2]。我們以包含兩個未知數和兩個方程式的線性方程組 A\mathbf{x}=\mathbf{b} 說明 Kaczmarz 算法的基本原理。考慮

\displaystyle \begin{bmatrix} a_{11}&a_{12}\\ a_{21}&a_{22} \end{bmatrix}\begin{bmatrix} x_1\\ x_2 \end{bmatrix}=\begin{bmatrix} b_1\\ b_2 \end{bmatrix}

其中係數 a_{ij} 和常數 b_i 是實數。當係數矩陣可逆時,線性方程的解為下列兩個超平面 (此例為 \mathbb{R}^2 平面的二直線) 的交點:

\displaystyle\begin{aligned} H_1&=\{(x_1,x_2)\vert a_{11}x_1+a_{12}x_2=b_1\}\\ H_2&=\{(x_1,x_2)\vert a_{21}x_1+a_{22}x_2=b_2\}.\end{aligned}

見下圖,給定任一初始點 \mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)})^T,連續交互正交投影至超平面 H_1H_2 可得一向量序列 \mathbf{x}^{(k)}k=0,1,2,\ldots。當 k\to\infty\mathbf{x}^{(k)}\to H_1\cap H_2,此即 A\mathbf{x}=\mathbf{b} 的解。

Kaczmarz method

Kaczmarz 算法

 
在向量空間 \mathbb{R}^n,任一點 \mathbf{v} 如何正交投影至超平面 H=\{\mathbf{x}\in\mathbb{R}^n\vert\mathbf{a}^T\mathbf{x}=b\}?非零向量 \mathbf{a} 是超平面 H 的法向量 (見“超平面”)。若 b\neq 0,超平面 H 未穿越原點,因此不是一個子空間,正式的名稱是仿射空間 (affine space)。超平面 H 上所有的點是 \mathbf{a}^T\mathbf{x}=b 的解。明顯地,b\mathbf{a}/\mathbf{a}^T\mathbf{a} 是一個特解,故仿射空間 H 可用通解表示為 (見“仿射組合與仿射空間”):

\displaystyle H=\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}+N(\mathbf{a}^T)=\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}+\text{span}\{\mathbf{a}\}^\perp

其中 \mathbf{a}^T 的零空間 (nullspace) N(\mathbf{a}^T) 等於 \mathbf{a} 所擴張的子空間 (\mathbb{R}^n 的一直線) 的正交補餘 \text{span}\{\mathbf{a}\}^\perp。類似一般子空間的正交投影,點 \mathbf{v}H 的仿射投影就是超平面 H 中距離 \mathbf{v} 最接近的點 \mathbf{p}。最直接的作法是將向量 \mathbf{v} 和仿射空間 H 都減去 b\mathbf{a}/\mathbf{a}^T\mathbf{a},如此一來,所有的物件又回到子空間 \text{span}\{\mathbf{a}\}^\perp (見“仿射投影”圖例)。既然 \mathbf{p}\mathbf{v}H 的仿射投影, \mathbf{p}-b\mathbf{a}/\mathbf{a}^T\mathbf{a} 也就是 \mathbf{v}-b\mathbf{a}/\mathbf{a}^T\mathbf{a}\text{span}\{\mathbf{a}\}^\perp 的正交投影。若 P 代表至 \text{span}\{\mathbf{a}\}^\perp 的正交投影矩陣,則

\displaystyle \mathbf{p}-\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}=P\left(\mathbf{v}-\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}\right)

因為 \mathbf{a}\mathbf{a}^T/\mathbf{a}^T\mathbf{a} 是子空間 \text{span}\{\mathbf{a}\} 的正交投影矩陣 (見“線代膠囊──正交投影矩陣”),可得 P=I-\mathbf{a}\mathbf{a}^T/\mathbf{a}^T\mathbf{a}。因此,

\displaystyle\begin{aligned} \mathbf{p}&=\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}+\left(I-\frac{\mathbf{a}\mathbf{a}^T}{\mathbf{a}^T\mathbf{a}}\right)\left(\mathbf{v}-\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}\right)\\ &=\frac{b\mathbf{a}}{\mathbf{a}^T\mathbf{a}}+\mathbf{v}-\frac{\mathbf{a}\mathbf{a}^T\mathbf{v}}{\mathbf{a}^T\mathbf{a}}\\ &=\mathbf{v}+\left(\frac{b-\mathbf{a}^T\mathbf{v}}{\Vert\mathbf{a}\Vert^2}\right)\mathbf{a}. \end{aligned}

 
考慮一般的線性方程 A\mathbf{x}=\mathbf{b},其中 A 是一 m\times n 階實矩陣,m\ge n\text{rank}A=n,且 \mathbf{b}=(b_1,\ldots,b_m)^T 屬於 A 的行空間 (column space)。在此條件下,A\mathbf{x}=\mathbf{b} 存在唯一解。令 \mathbf{a}_i\in\mathbb{R}^n 為係數矩陣 A 的第 i 個列向量 (row vector),i=1,\ldots,m。寫出

\displaystyle A\mathbf{x}=\begin{bmatrix} \mathbf{a}_1^T\\ \vdots\\ \mathbf{a}_m^T \end{bmatrix}\mathbf{x}=\begin{bmatrix} \mathbf{a}_1^T\mathbf{x}\\ \vdots\\ \mathbf{a}_m^T\mathbf{x} \end{bmatrix}=\begin{bmatrix} b_1\\ \vdots\\ b_m \end{bmatrix}=\mathbf{b}

此方程式的解即為 \mathbb{R}^nm 個超平面的交集 H_1\cap\cdots\cap H_m,其中

\displaystyle H_i=\{\mathbf{x}\in\mathbb{R}^n\vert\mathbf{a}_i^T\mathbf{x}=b_i\},~~i=1,\ldots,m

給定一初始向量 \mathbf{x}^{(0)}\in\mathbb{R}^n,連續執行仿射投影至每一超平面,可得向量序列:

\displaystyle\begin{aligned} \mathbf{x}^{(1)}&=\mathbf{x}^{(0)}+\frac{b_1-\mathbf{a}_1^T\mathbf{x}^{(0)}}{\Vert\mathbf{a}_1\Vert^2}\mathbf{a}_1\\ \mathbf{x}^{(2)}&=\mathbf{x}^{(1)}+\frac{b_2-\mathbf{a}_2^T\mathbf{x}^{(1)}}{\Vert\mathbf{a}_2\Vert^2}\mathbf{a}_2\\ &~~\vdots\\ \mathbf{x}^{(m)}&=\mathbf{x}^{(m-1)}+\frac{b_m-\mathbf{a}_m^T\mathbf{x}^{(m-1)}}{\Vert\mathbf{a}_m\Vert^2}\mathbf{a}_m. \end{aligned}

當所有的超平面都用罄後,重複此程序將 \mathbf{x}^{(m)} 投影至 H_1 可得 \mathbf{x}^{(m+1)},再將 \mathbf{x}^{(m+1)} 投影至 H_2 可得 \mathbf{x}^{(m+2)},餘此類推。Kaczmarz 算法可用雙迴圈表示如下:

標準 Kaczmarz 算法


設初始向量為 \mathbf{x}^{(0)}
對於 j=0,1,2,\ldots
    對於 i=1,2,\ldots,m

\displaystyle \mathbf{x}^{(mj+i)}\leftarrow\mathbf{x}^{(mj+i-1)}+\frac{b_i-\mathbf{a}_i^T\mathbf{x}^{(mj+i-1)}}{\Vert\mathbf{a}_i\Vert^2}\mathbf{a}_i


 
下面證明 \mathbf{x}^{(mj+i)} 收斂至 A\mathbf{x}=\mathbf{b} 的解。為簡化記號,設 \delta=(b_i-\mathbf{a}_i^T\mathbf{x}^{(mj+i-1)})/\Vert\mathbf{a}_i\Vert^2,則 \mathbf{x}^{(mj+i)}=\mathbf{x}^{(mj+i-1)}+\delta\mathbf{a}_i。計算

\displaystyle\begin{aligned} \left\|\mathbf{x}^{(mj+i)}-\mathbf{x}\right\|^2&=\left\|\mathbf{x}^{(mj+i-1)}+\delta\mathbf{a}_i-\mathbf{x}\right\|^2\\ &=\left\|\mathbf{x}^{(mj+i-1)}-\mathbf{x}\right\|^2+2\delta\mathbf{a}_i^T(\mathbf{x}^{(mj+i-1)}-\mathbf{x})+\delta^2\Vert\mathbf{a}_i\Vert^2, \end{aligned}

其中第二項可化簡為

\displaystyle\begin{aligned} \mathbf{a}_i^T(\mathbf{x}^{(mj+i-1)}-\mathbf{x})&=\mathbf{a}_i^T\mathbf{x}^{(mj+i-1)}-\mathbf{a}_i^T\mathbf{x}\\ &=\mathbf{a}_i^T\mathbf{x}^{(mj+i-1)}-b_i=-\delta\Vert\mathbf{a}_i\Vert^2. \end{aligned}

代回前一式可得

\displaystyle \left\|\mathbf{x}^{(mj+i)}-\mathbf{x}\right\|^2=\left\|\mathbf{x}^{(mj+i-1)}-\mathbf{x}\right\|^2-\delta^2\Vert\mathbf{a}_i\Vert^2

故知 \left\|\mathbf{x}^{(mj+i)}-\mathbf{x}\right\|\le\left\|\mathbf{x}^{(mj+i-1)}-\mathbf{x}\right\|,等號發生於 \delta=0,亦即 \mathbf{x}^{mj+i-1}\in H_{i-1}\cap H_i。所以,Kaczmarz 算法產生的單調遞減數列 \left\|\mathbf{x}^{(mj+i)}-\mathbf{x}\right\| 存在一極值,換句話說,每一列 i 所定義的 \delta 都趨於零,也就證明 \mathbf{x}^{(mj+i)} 收斂至 \mathbf{x}

 
對於大尺寸問題,標準 Kaczmarz 算法的收斂速度不甚理想。實驗印證隨機選擇投影超平面可以有效地提升收斂速度,作法如下[3]

隨機 Kaczmarz 算法


設初始向量為 \mathbf{x}^{(0)}
對於 k=0,1,2,\ldots
    隨機選擇 i\in\{1,2,\ldots,m\},機率大小正比於 \Vert\mathbf{a}_i\Vert^2

\displaystyle \mathbf{x}^{(k)}\leftarrow\mathbf{x}^{(k-1)}+\frac{b_i-\mathbf{a}_i^T\mathbf{x}^{(k-1)}}{\Vert\mathbf{a}_i\Vert^2}\mathbf{a}_i


 
參考來源:
[1] 維基百科:Kaczmarz method
[2] 維基百科:Algebraic Reconstruction Technique
[3] Thomas Strohmer and Roman Vershynin, A randomized Kaczmarz algorithm with
exponential convergence (PDF), Journal of Fourier Analysis and Applications, Vol. 15, pp 262-278, 2009.

相關閱讀:
This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , . Bookmark the permalink.

3 則回應給 Kaczmarz 算法

  1. 張盛東 說:

    老師,請問一下這種迭代算法對比其他的迭代算法有何優勢?

  2. edge 說:

    看了一下,覺得此法搞不好非常省暫存記憶體
    對於求解大型稀疏矩陣的話

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s