特殊矩陣 (4)：Householder 矩陣

Posted on 09/14/2009 by ccjou

本文的閱讀等級：中級

在幾何向量空間 $\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 矩陣即為鏡射矩陣

鏡射變換也可以用正交投影推導出來。先求出給定向量 $\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}=I$ ， $H$ 稱為對合 (involutory) 矩陣。一函數 $f$ 若同時為其反函數，也就是說，對於任意 $x$ ， $f(f(x))=x$ ， $f$ 稱為對合函數。

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}_i$ ， $i=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$ 。

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

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_2$ 是 $2\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)=1$ 若 $x_1\ge 0$ ， $\hbox{sign}(x_1)=-1$ 若 $x_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.

繼續閱讀：

This entry was posted in 特殊矩陣, 線性代數專欄 and tagged 特殊矩陣, 鏡射矩陣, 超平面, Householder 矩陣, QR 分解, QR 演算法, 分塊矩陣, 對合矩陣, 正交矩陣, 三對角化. Bookmark the permalink.

12 Responses to 特殊矩陣 (4)：Householder 矩陣

Vincent says:

06/15/2011 at 7:46 am

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.

Reply
ccjou says:

06/15/2011 at 8:41 am

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.

Reply
- Vahi Chen says:
  
  09/18/2015 at 1:49 am
  
  准确地说， $\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。
  
  Reply
  - ccjou says:
    
    09/18/2015 at 7:31 am
    
    謝謝你的補充解說，非常的精確。我將它貼到上文的附註裡。
    
    Reply
Kara says:

07/13/2011 at 12:48 am

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

Reply
ccjou says:

07/13/2011 at 6:11 am

確實我們無法想像 $\mathbb{R}^n (n>3)$ 的情況，所以還是由 $\mathbb{R}^3$ 取得幾何解釋。在 $\mathbb{R}^3$ 中穿越原點的一直線 $L$ 構成一子空間， $\mathrm{dim}L=1$ ，與此直線垂直的所有向量構成一穿越原點的平面 $P$ 。幾何學稱 $L$ 指出 $P$ 的法向量 (normal vector)，線性代數語言則為 $P$ 中任意向量正交於 $L$ 上所有向量，我們稱 $P$ 是 $L$ 的正交補集 (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/

Reply
npes_87184 says:

07/01/2013 at 8:37 am

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

Reply
- ccjou says:
  
  07/01/2013 at 9:33 am
  
  如果是矩陣表達的反射，必定通過原點，但仿射變換則否。請查閱相關文章。
  
  Reply
wonderlandtommy says:

09/23/2015 at 12:47 pm

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

Reply
wonderlandtommy says:

09/23/2015 at 12:49 pm

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

Reply
- ccjou says:
  
  09/23/2015 at 2:02 pm
  
  我在上文增添了推導細節，請查閱。
  
  Reply
  - wonderlandtommy says:
    
    09/24/2015 at 12:18 pm
    
    謝謝ccjou老師
    
    Reply

	陳倍恩 on 線性代數的第一堂課──矩陣乘法的定義
	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨