Householder变换

Householder变换

已知$\mathbf{x}$,做Householder变换，使得$\Vert \mathbf{x}\Vert =\Vert \mathbf{y}\Vert$,$\mathbf{Hx}=\mathbf{y}$,
要求$\mathbf{H}{\mathbf{H}}^H=\mathbf{I},{\mathbf{H}}^H=\mathbf{H}$
于是$\mathbf{H}=\mathbf{I}-\mathbf{u}{\mathbf{u}}^H$
其中$\mathbf{u}=\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}\Vert }$

推导：

$$\begin{array}{rl}\mathbf{y}& =\mathbf{x}-\left(\mathbf{x}-\mathbf{y}\right)\\ & =\mathbf{x}-2\frac{\mathbf{x}-\mathbf{y}}{2}\\ & =\mathbf{x}-2\mathbf{Px}\\ & =\mathbf{x}-2\mathbf{u}{\mathbf{u}}^H\mathbf{x}\\ & =\left(\mathbf{I}-2\mathbf{u}{\mathbf{u}}^H\right)\mathbf{x}\end{array}$$
其中$\mathbf{P}$是 span ⁡ { x − y } \operatorname{span}\lbrace \mathbf{x}-\mathbf{y} \rbrace span{x−y}的正交投影
span ⁡ { x − y } \operatorname{span}\lbrace \mathbf{x}-\mathbf{y} \rbrace span{x−y}的标准正交基$\mathbf{u}=\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}\Vert }$
$\mathbf{P}=\mathbf{u}{\mathbf{u}}^H$

变为e1

即$\mathbf{Hx}=\alpha {\mathbf{e}}_1$
容易求出来$\alpha =\pm \Vert \mathbf{x}\Vert$
可以取$\alpha =\mathrm{sign}\left(x_1\right)$，于是$\mathbf{u}=\frac{\mathbf{x}+\mathrm{sign}\left(x_1\right)\Vert \mathbf{x}\Vert {\mathbf{e}}_1}{\Vert \mathbf{x}+\mathrm{sign}\left(x_1\right)\Vert \mathbf{x}\Vert {\mathbf{e}}_1\Vert }$

减小误差

可能搜代码的时候会有这种
$$\mathbf{H}=\mathbf{I}-\tau \mathbf{u}{\mathbf{u}}^T,\tau =\frac{2}{{\mathbf{u}}^T\mathbf{u}}$$
其中$\mathbf{u}=\mathbf{x}-\alpha {\mathbf{e}}_1$,但是
当$x_1>0$时
$$u_1=x_1-\Vert \mathbf{x}\Vert =\frac{x_1^2-\Vert \mathbf{x}{\Vert }^2}{x_1+\Vert \mathbf{x}\Vert }=-\frac{x_2^2+\cdots +x_n^2}{x_1+\Vert \mathbf{x}\Vert }$$
当$x_1\le 0$时,依然用原来的方法算
$$u_1=x_1-\Vert \mathbf{x}\Vert$$
然后
$$\begin{gathered} \tau =\frac{2u_1^2}{\Vert \mathbf{u}{\Vert }^2} \\ \mathbf{u}=\frac{\mathbf{u}}{u_1} \end{gathered}$$
代码来自https://stackoverflow.com/questions/53489237/how-can-you-implement-householder-based-qr-decomposition-in-python

def householder(x: np.ndarray) -> Union[np.ndarray, int]:
    alpha = x[0]
    s = np.power(np.linalg.norm(x[1:]), 2)
    v = x.copy()

    if s == 0:
        tau = 0
    else:
        t = np.sqrt(alpha**2 + s)
        v[0] = alpha - t if alpha <= 0 else -s / (alpha + t)

        tau = 2 * v[0]**2 / (s + v[0]**2)
        v /= v[0]

    return v, tau

第二种

$$\mathbf{H}=\mathbf{I}-\tau \mathbf{u}{\mathbf{u}}^T,\tau =\frac{2}{{\mathbf{u}}^T\mathbf{u}}$$
其中$\mathbf{u}=\mathbf{x}+\mathrm{sign}\left(x_1\right)\Vert \mathbf{x}\Vert {\mathbf{e}}_1$
然后
$$\begin{gathered} \mathbf{u}=\frac{\mathbf{u}}{u_1} \\ \tau =\frac{2}{\Vert \mathbf{u}{\Vert }^2} \end{gathered}$$
https://stackoverflow.com/questions/53489237/how-can-you-implement-householder-based-qr-decomposition-in-python

def householder_vectorized(a):
    """Use this version of householder to reproduce the output of np.linalg.qr 
    exactly (specifically, to match the sign convention it uses)
    
    based on https://rosettacode.org/wiki/QR_decomposition#Python
    """
    v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
    v[0] = 1
    tau = 2 / (v.T @ v)
    
    return v,tau

QR分解

QR分解可以用施密特正交化做，也可以用householder变换
以4维为例

$$\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right)$$
目标是变为上三角矩阵
先利用Householder变换，变换第一列
设$\mathbf{x}=\left(\begin{array}{c}{a}_{11}\\ a_{21}\\ a_{31}\\ a_{41}\end{array}\right)$,$\mathbf{y}=\left(\begin{array}{c}{a}_{11}^{(1)}\\ 0\\ 0\\ 0\end{array}\right)$

根据Householder变换，$a_{11}^{(1)}=\Vert \mathbf{x}\Vert$
$$\mathbf{u}=\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}\Vert }$$
$${\mathbf{H}}^{(1)}=\mathbf{I}-\mathbf{u}{\mathbf{u}}^H$$
$${\mathbf{A}}^{(1)}={\mathbf{H}}^{(1)}\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& a_{14}^{(1)}\\ 0& a_{22}^{(1)}& a_{23}^{(1)}& a_{24}^{(1)}\\ 0& a_{32}^{(1)}& a_{33}^{(1)}& a_{34}^{(1)}\\ 0& a_{42}^{(1)}& a_{43}^{(1)}& a_{44}^{(1)}\end{array}\right)$$
接着变换第二列，此时要保证不改变第一行

设$\mathbf{x}=\left(\begin{array}{c}{a}_{22}^{(1)}\\ a_{32}^{(1)}\\ a_{42}^{(1)}\end{array}\right)$,$\mathbf{y}=\left(\begin{array}{c}{a}_{22}^{(2)}\\ 0\\ 0\end{array}\right)$
$$\mathbf{u}=\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}\Vert }$$
$${\tilde{\mathbf{H}}}^{(2)}=\mathbf{I}-\mathbf{u}{\mathbf{u}}^H$$
$${\mathbf{H}}^{(2)}=\left(\begin{array}{cc}1& 0^H\\ 0& {\tilde{\mathbf{H}}}^{(2)}\end{array}\right)$$
$${\mathbf{A}}^{(2)}={\mathbf{H}}^{(2)}{\mathbf{H}}^{(1)}\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& a_{14}^{(1)}\\ 0& a_{22}^{(2)}& a_{23}^{(2)}& a_{24}^{(2)}\\ 0& 0& a_{33}^{(2)}& a_{34}^{(2)}\\ 0& 0& a_{43}^{(2)}& a_{44}^{(2)}\end{array}\right)$$
接着变换第三列，同时保证第一行，第二行不改变
设$\mathbf{x}=\left(\begin{array}{c}{a}_{33}^{(2)}\\ a_{43}^{(2)}\end{array}\right)$,$\mathbf{y}=\left(\begin{array}{c}{a}_{33}^{(3)}\\ 0\end{array}\right)$
$$\mathbf{u}=\frac{\mathbf{x}-\mathbf{y}}{\Vert \mathbf{x}-\mathbf{y}\Vert }$$
$${\tilde{\mathbf{H}}}^{(3)}=\mathbf{I}-\mathbf{u}{\mathbf{u}}^H$$
$${\mathbf{H}}^{(3)}=\left(\begin{array}{cc}{\mathbf{I}}_{2\times 2}& 0^H\\ 0& {\tilde{\mathbf{H}}}^{(3)}\end{array}\right)$$
$$\mathbf{R}={\mathbf{A}}^{(3)}={\mathbf{H}}^{(3)}{\mathbf{H}}^{(2)}{\mathbf{H}}^{(1)}\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}^{(1)}& a_{12}^{(1)}& a_{13}^{(1)}& a_{14}^{(1)}\\ 0& a_{22}^{(2)}& a_{23}^{(2)}& a_{24}^{(2)}\\ 0& 0& a_{33}^{(3)}& a_{34}^{(3)}\\ 0& 0& 0& a_{44}^{(3)}\end{array}\right)$$
于是
$$\mathbf{Q}={\mathbf{H}}^{(1)}{\mathbf{H}}^{(2)}{\mathbf{H}}^{(3)}$$