特殊矩陣 (4):Householder 矩陣

本文的閱讀等級:中級

在幾何向量空間 \mathbb{R}^n,鏡射 (reflection) 超平面 (hyperplane) 由單位法向量 \mathbf{v} 決定 (超平面是維數等於 n-1 的子空間)。對於任一點 (或向量) \mathbf{x},從空間幾何可推論點 \mathbf{x} 的鏡射為 \mathbf{x}-2(\mathbf{v}^T\mathbf{x})\mathbf{v},其中 \mathbf{v}^T\mathbf{x} 是點 \mathbf{x} 在法向量 \mathbf{v} 的投影量 (見下圖)。將純量 \mathbf{v}^T\mathbf{x} 和向量 \mathbf{v} 交換位置,並移除括弧,鏡射點可表示為

\mathbf{x}-2\mathbf{v}\mathbf{v}^T\mathbf{x}=(I-2\mathbf{v}\mathbf{v}^T)\mathbf{x}

單位法向量 \mathbf{v} 所定義的超平面的鏡射變換矩陣稱為 Householder 矩陣,如下:

H=I-2\mathbf{ v}\mathbf{v}^T

其中 \Vert\mathbf{v}\Vert=1。(對於 \mathbb{R}^n 的一般子空間的鏡射矩陣請見“旋轉與鏡射”。)

householder-matrix

Householder 矩陣即為鏡射矩陣

 
鏡射變換也可以用正交投影推導出來。先求出給定向量 \mathbf{x} 至法向量 \mathbf{v} 的正交投影 P\mathbf{x},其中 P 是正交投影矩陣 (見“正交投影──威力強大的線代工具”),如下:

\displaystyle P=\frac{\mathbf{v}\mathbf{v}^T}{\mathbf{v}^T\mathbf{v}}=\mathbf{v}\mathbf{v}^T

因為鏡射超平面是法向量的正交補餘 (orthogonal complement),\mathbf{x} 至鏡射面的正交投影即為 \mathbf{x}-P\mathbf{x}=(I-P)\mathbf{x},所以 \mathbf{x} 的鏡射點是

(I-P)\mathbf{x}-P\mathbf{x}=(I-2P)\mathbf{x}=(I-2\mathbf{v}\mathbf{v}^T)\mathbf{x}

 
Householder 矩陣 H 具有以下性質:H 是對稱矩陣也是正交矩陣。直接計算並利用 \mathbf{v}^T\mathbf{v}=1 可證明:

H^TH=(I-2\mathbf{v}\mathbf{v}^T)(I-2\mathbf{v}\mathbf{v}^T)=I-4\mathbf{v}\mathbf{v}^T+4\mathbf{v}\mathbf{v}^T\mathbf{v}\mathbf{v}^T=I

因此 H=H^T=H^{-1},故 H^2=HH^{-1}=IH 稱為對合 (involutory) 矩陣。一函數 f 若同時為其反函數,也就是說,對於任意 xf(f(x))=xf 稱為對合函數。

Alston Scott Householder (1904-1993) From http://apprendre-math.info/history/photos/Householder_2.jpeg

 
利用 Householder 矩陣的幾何意涵很容易可以求得其特徵值。首先考慮向量 \mathbf{v} 的鏡射:

H\mathbf{v}=(I-2\mathbf{v}\mathbf{v}^T)\mathbf{v}=\mathbf{v}-2\mathbf{v}=-\mathbf{v}

故知 H 有一特徵向量 \mathbf{v} 對應特徵值 -1。因為鏡射超平面的維數等於 n-1,令此平面上的任意 n-1 個線性獨立向量為 \mathbf{u}_ii=1,2,\ldots,n-1,它們必定滿足 \mathbf{v}^T\mathbf{u}_i=0,即得

H\mathbf{u}_i=(I-2\mathbf{v}\mathbf{v}^T)\mathbf{u}_i=\mathbf{u}_i

這說明了 Hn-1 個重複特徵值 1,故 \mathrm{det}H=-1H 是可逆矩陣。

 
Householder 矩陣的潛藏功能在於一旦選擇了適當的超平面法向量 \mathbf{v},鏡射向量 H\mathbf{x} 將指向標準基底向量 \mathbf{e}_1=(1,0,\ldots,0)^T,也就是說,除了第一個元外,H\mathbf{x} 的每個元皆為 0。詳細的算法見註解[1],下面給出一個簡約版本。令 \sigma=\Vert\mathbf{x}\Vert。設

\displaystyle  \mathbf{v}=\frac{\mathbf{x}+\sigma\mathbf{e}_1}{\Vert\mathbf{x}+\sigma\mathbf{e}_1\Vert}

則鏡射向量為

\displaystyle  H\mathbf{x}=\mathbf{x}-2\mathbf{v}\mathbf{v}^T\mathbf{x}=\mathbf{x}-(\mathbf{x}+\sigma\mathbf{e}_1)\frac{2(\mathbf{x}+\sigma\mathbf{e}_1)^T\mathbf{x}}{(\mathbf{x}+\sigma\mathbf{e}_1)^T(\mathbf{x}+\sigma\mathbf{e}_1)}

利用 \mathbf{x}^T\mathbf{x}=\sigma^2,可推得 2(\mathbf{x}+\sigma\mathbf{e}_1)^T\mathbf{x}=(\mathbf{x}+\sigma\mathbf{e}_1)^T(\mathbf{x}+\sigma\mathbf{e}_1)。所以,

H\mathbf{x}=\mathbf{x}-(\mathbf{x}+\sigma\mathbf{e}_1)=-\sigma\mathbf{e}_1

到目前為止,讀者或許還看不出來這其實是 Householder 矩陣最具實用價值的性質。

 
Householder 變換在數值線性代數有兩個主要應用:執行 QR 分解與計算特徵值的 QR 迭代法的第一個步驟。詳細的過程將另文說明,在此我們先研究 Householder 變換的基本應用技巧:將對稱矩陣 A 予以三對角化 (tridiagonalization),目的是盡可能產生零元:

\begin{bmatrix}    \ast&\ast&0&0&0\\    \ast&\ast&\ast&0&0\\    0&\ast&\ast&\ast&0\\    0&0&\ast&\ast&\ast\\    0&0&0&\ast&\ast    \end{bmatrix}

 
看下面的例子[2]

A=\left[\!\!\begin{array}{rcrr}    4&1&-2&2\\    1&2&0&1\\    -2&0&3&-2\\    2&1&-2&-1    \end{array}\!\!\right]

我們的想法是對 A 執行相似變換 P^{-1}AP,變換矩陣 P^{-1}AP 保有與 A 相同的特徵值。同時限制 P 是正交矩陣,這麼做的好處是逆矩陣可以立刻求得 P^{-1}=P^T,而且 A 的一些重要性質也被保留下來,譬如,當 A 是對稱矩陣時,P^{-1}AP 也是對稱矩陣:

(P^{-1}AP)^T=P^TA^TP=P^{-1}AP

使用前述鏡射性質,我們設計下面的主對角分塊正交矩陣:

P_1=\begin{bmatrix}    1&\vline& 0&0&0\\ \hline    0&\vline&~&~&~\\    0&\vline&~&H_1&~\\    0&\vline&~&~&~    \end{bmatrix}

其中 H_1 是一 3\times 3 階 Householder 矩陣。考慮分塊矩陣乘法

P_1AP_1=\begin{bmatrix}  1&0\\  0&H_1  \end{bmatrix}\begin{bmatrix}  a_{11}&\mathbf{a}^T\\  \mathbf{a}&\tilde{A}  \end{bmatrix}\begin{bmatrix}  1&0\\  0&H_1  \end{bmatrix}=\begin{bmatrix}  a_{11}&\mathbf{a}^TH_1\\  H_1\mathbf{a}&H_1\tilde{A}H_1  \end{bmatrix}

其中 \mathbf{a}A=[a_{ij}] 的第一行的後面三個元構成的向量,即 \mathbf{a}=\begin{bmatrix}  a_{21}\\  a_{31}\\  a_{41}  \end{bmatrix}\tilde{A}3\times 3 階分塊。我們的目標是使 H_1\mathbf{a} 平行於 \mathbf{e}_1=\begin{bmatrix}  1\\  0\\  0  \end{bmatrix},故設 \mathbf{x}_1=\mathbf{a}=\left[\!\!\begin{array}{r}    1\\    -2\\    2    \end{array}\!\!\right]\Vert\mathbf{x}_1\Vert=3。將 \mathbf{x}_1 代入前述公式算出 \mathbf{v}_1=\displaystyle\frac{1}{\sqrt{6}}\left[\!\!\begin{array}{r}    2\\    -1\\    1    \end{array}\!\!\right],再代入 Householder 矩陣公式可得

H_1=\left[\!\!\begin{array}{rcr}    -\frac{1}{3}&\frac{2}{3}&-\frac{2}{3}\\ [0.3em]    \frac{2}{3}&\frac{2}{3}&\frac{1}{3}\\ [0.3em]    -\frac{2}{3}&\frac{1}{3}&\frac{2}{3}    \end{array}\!\!\right]

因為 H=H^T=H^{-1},直接計算可知 P_1=P_1^T=P_1^{-1}P_1 也是對合矩陣,並得到

A_1=P_1AP_1=\left[\!\!\begin{array}{rrrr}    4&-3&0&0\\ [0.3em]    -3&\frac{10}{3}&1&\frac{4}{3}\\ [0.3em]    0&1&\frac{5}{3}&-\frac{4}{3}\\ [0.3em]    0&\frac{4}{3}&-\frac{4}{3}&-1    \end{array}\!\!\right]

 
重複上述步驟,但僅考慮 P_1AP_1 的右下 3\times 3 分塊,H_22\times 2 矩陣,P_2 則為

P_2=\left[\!\!\begin{array}{ccrr}    1&0&0&0\\    0&1&0&0\\ [0.3em]    0&0&-\frac{3}{5}&-\frac{4}{5}\\ [0.3em]    0&0&-\frac{4}{5}&\frac{3}{5}    \end{array}\!\!\right]

再計算

A_2=P_2A_1P_2=\left[\!\!\begin{array}{rrrr}    4&-3&0&0\\ [0.3em]    -3&\frac{10}{3}&-\frac{5}{3}&0\\ [0.3em]    0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\\ [0.3em]    0&0&\frac{68}{75}&\frac{149}{75}    \end{array}\!\!\right]

所以,A_2=P^{-1}AP 為一三對角矩陣,其中 P=P_1P_2 是正交矩陣,因為 P^{-1}=P_2^{-1}P_1^{-1}=P_2^TP_1^T=(P_1P_2)^T=P^T

 
註解
[1] 設非零向量 \mathbf{u}\in\mathbb{R}^n,將 Householder 矩陣改寫為

\displaystyle  H(\mathbf{u})=I-2\frac{\mathbf{u}\mathbf{u}^T}{\mathbf{u}^T\mathbf{u}}

使之符合實際應用場合。考慮 \mathbf{x},\mathbf{y}\in\mathbb{R}^n,且 \Vert\mathbf{x}\Vert=\Vert\mathbf{y}\Vert。令 \mathbf{u}=\mathbf{y}+\mathbf{x}\mathbf{w}=\mathbf{y}-\mathbf{x}。因此,\mathbf{u}^T\mathbf{w}=(\mathbf{y}+\mathbf{x})^T(\mathbf{y}-\mathbf{x})=\Vert\mathbf{y}\Vert^2-\Vert\mathbf{x}\Vert^2=0,並有

\displaystyle  H(\mathbf{u})\mathbf{u}=\left(I-2\frac{\mathbf{u}\mathbf{u}^T}{\mathbf{u}^T\mathbf{u}}\right)\mathbf{u}=-\mathbf{u}

\displaystyle  H(\mathbf{u})\mathbf{w}=\left(I-2\frac{\mathbf{u}\mathbf{u}^T}{\mathbf{u}^T\mathbf{u}}\right)\mathbf{w}=\mathbf{w}

所以,

\displaystyle\begin{aligned}  H(\mathbf{u})\mathbf{y}-H(\mathbf{u})\mathbf{x}&=H(\mathbf{u})(\mathbf{y}-\mathbf{x})=H(\mathbf{u})\mathbf{w}=\mathbf{w}=\mathbf{y}-\mathbf{x}\\  H(\mathbf{u})\mathbf{y}+H(\mathbf{u})\mathbf{x}&=H(\mathbf{u})(\mathbf{y}+\mathbf{x})=H(\mathbf{u})\mathbf{u}=-\mathbf{u}=-\mathbf{y}-\mathbf{x}.\end{aligned}

令上面兩式相減可得 H(\mathbf{u})\mathbf{x}=-\mathbf{y}。使用相同方法,亦可推得 H(\mathbf{w})\mathbf{x}=\mathbf{y}。兩個 Householder 矩陣的差異在於目標向量的方向相反。

 
給定非零向量 \mathbf{x}=(x_1,\ldots,x_n)^T\in\mathbb{R}^n,我們希望找出 \mathbf{v}\in\mathbb{R}^n 使得 H(\mathbf{v})\mathbf{x} 與標準單位向量 \mathbf{e}_1 平行。定義 \hbox{sign}(x_1)=1x_1\ge 0\hbox{sign}(x_1)=-1x_1<0。令 \mathbf{y}=\hbox{sign}(x_1)\Vert\mathbf{x}\Vert\mathbf{e}_1。沿用前述記號,\mathbf{u}=\mathbf{y}+\mathbf{x}=\hbox{sign}(x_1)\Vert\mathbf{x}\Vert\mathbf{e}_1+\mathbf{x}\mathbf{w}=\mathbf{y}-\mathbf{x}=\hbox{sign}(x_1)\Vert\mathbf{x}\Vert\mathbf{e}_1-\mathbf{x}。當 x_1\neq 0\Vert\mathbf{u}\Vert>\Vert\mathbf{w}\Vert。為了讓數值計算較為穩定,我們挑選向量長度大者,故設 \mathbf{v}=\mathbf{u}。因此,H(\mathbf{v})\mathbf{x}=-\mathbf{y}=-\hbox{sign}(x_1)\Vert\mathbf{x}\Vert\mathbf{e}_1

[2] 引用來源 Richard L. Burden and J. Douglas Faires, Numerical Analysis, 9th edition, 2010.

繼續閱讀:
Advertisements
本篇發表於 特殊矩陣, 線性代數專欄 並標籤為 , , , , , , , , , 。將永久鏈結加入書籤。

12 則回應給 特殊矩陣 (4):Householder 矩陣

  1. Vincent 說道:

    hello. I think the reflection vector, where you wrote as V=(x+sigma e1)/||x+ sigma e1||, it should be (x- sigma e1)/||x-sigma e1||?

    because V should be perpendicular to the half angle line of x and sigma*e1, am I right?

    Thank you.

  2. ccjou 說道:

    The question you raised is quite interesting. Technically, both x-\sigma{e}_1 and x+\sigma{e}_1 will serve our purpose, i.e., to eleminate all the entries except the first entry of x. See the following article for details
    https://ccjou.wordpress.com/2011/05/24/householder-%E8%AE%8A%E6%8F%9B%E6%96%BC-qr-%E5%88%86%E8%A7%A3%E7%9A%84%E6%87%89%E7%94%A8/

    I will give a brief explanation here. Hope it can clear the doubt. Let x=(3,4) and e_1=(1,0). If you choose v_1=x-\sigma e_1=(3,4)-5(1,0)=(-2,4) and normalize it, then H_1 will reflect x and make it align with e_1. But, if you choose v_2=x+\sigma e_1=(8,4), then H_2 will reflect x and make it align with -e_1! In other words, H_2 simply reflects x to the opposite direction of e_1. It should be helpful if you draw a figure to see how the above reflections work.

    • Vahi Chen 說道:

      准确地说,\sigma前面使用加号还是减号,会决定\sigma本身的正负性,即:

      1. 若取\mathbf{v}=\dfrac{\mathbf{x}-\sigma\mathbf{e}_1}{\lVert\mathbf{x}-\sigma\mathbf{e}_1\rVert},则定义\sigma=-\text{sign}(x_1)\lVert\mathbf{x}\rVert_2

      2. 若取\mathbf{v}=\dfrac{\mathbf{x}+\sigma\mathbf{e}_1}{\lVert\mathbf{x}+\sigma\mathbf{e}_1\rVert},则定义\sigma=\text{sign}(x_1)\lVert\mathbf{x}\rVert_2

      理由是在构造\mathbf{v}的过程中,避免同符号的数相减引起有效位数损失,导致\mathbf{v}不能定义或计算结果不精确。对于情形1,\mathbf{v}的分母\mathbf{x}-\sigma\mathbf{e}_1=\mathbf{x}+\text{sign}(x_1)\lVert\mathbf{x}\rVert_2\mathbf{e}_1=\mathbf{w},容易验证,\lVert\mathbf{w}\rVert^2_2=2\lVert\mathbf{x}\rVert_2\left(\lVert\mathbf{x}\rVert_2+\lvert{x_1}\rvert\right)>0,同理可证情形2。

      考虑:\mathbf{x}=(3,0),\mathbf{e}_1=(1,0),若定义\mathbf{v}=\dfrac{\mathbf{x}-\sigma\mathbf{e}_1}{\lVert\mathbf{x}-\sigma\mathbf{e}_1\rVert_2},则必须定义\sigma=-\text{sign}(x_1)\lVert\mathbf{x}\rVert_2=-3,得到\mathbf{v}=\dfrac{(3,0)+3(1,0)}{\lVert(3,0)+3(1,0)\rVert_2},从而保证分母不为零。

      最后需要注意的是,若x_1=0,可以取\text{sign}(x_1)为1或-1。

  3. Kara 說道:

    您好, 感谢您的文章,非常清楚!
    能否请教一个问题: 为何在n维空间里取一个平面,则此平面是n-1维的?(直观上看确实如此,但是如何用数学语言叙述?)
    谢谢!

  4. ccjou 說道:

    確實我們無法想像 \mathbb{R}^n (n>3) 的情況,所以還是由 \mathbb{R}^3 取得幾何解釋。在 \mathbb{R}^3 中穿越原點的一直線 L 構成一子空間,\mathrm{dim}L=1,與此直線垂直的所有向量構成一穿越原點的平面 P。幾何學稱 L 指出 P 的法向量 (normal vector),線性代數語言則為 P 中任意向量正交於 L 上所有向量,我們稱 PL 的正交補集 (orthogonal complement),記為 P=L^{\perp},則有 \mathrm{dim}L+\mathrm{dim}P=3,因為 \mathrm{dim}\mathbb{R}^3=3。接下來,將 3 換成 n,將平面換成超平面即可。

    正交補集的討論請參考
    https://ccjou.wordpress.com/2011/05/19/%E6%AD%A3%E4%BA%A4%E8%A3%9C%E9%9B%86%E8%88%87%E6%8A%95%E5%BD%B1%E5%AE%9A%E7%90%86/

  5. npes_87184 說道:

    他所對應的對稱,那個超平面必定要通過原點?

  6. wonderlandtommy 說道:

    請問在三對角化的過程中,3階householder矩陣H1必須由A的第一行的後三元生成?

  7. wonderlandtommy 說道:

    請問在三對角化的過程中,為什麼H1必須由A的第一行的後三元生成?

發表迴響

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

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s