Householder reflection

☆ Save On Wikipedia ↗

In linear algebra, a Householder transformation (also known as a Householder reflection or elementary reflector) is a linear transformation that describes a reflection about a plane or hyperplane containing the origin. The Householder transformation was used in a 1958 paper by Alston Scott Householder[1] as a way to speed up the stable solution of linear systems.

Definition

Operator and transformation

The Householder operator[2] may be defined over any finite-dimensional inner product space V {\displaystyle V} {\displaystyle V} with inner product ⟨ ⋅ , ⋅ ⟩ {\displaystyle \langle \cdot ,\cdot \rangle } {\displaystyle \langle \cdot ,\cdot \rangle } and unit vector u ∈ V {\displaystyle u\in V} {\displaystyle u\in V} as

H u ( x ) := x − 2 ⟨ x , u ⟩ u . {\displaystyle H_{u}(x):=x-2\,\langle x,u\rangle \,u\,.} {\displaystyle H_{u}(x):=x-2\,\langle x,u\rangle \,u\,.}[3]

As defined here, the inner product ⟨ ⋅ , ⋅ ⟩ {\displaystyle \langle \,\cdot {\text{,}}\,\cdot \rangle } {\displaystyle \langle \,\cdot {\text{,}}\,\cdot \rangle } is linear in its first argument, and antilinear in its second argument, so that if a and b are scalars, then ⟨ a x , b y ⟩ = a b ¯ ⟨ x , y ⟩ {\displaystyle \langle \,a\,x{\text{,}}\,b\,y\rangle \,=\,a\,{\overline {b}}\,\langle x{\text{,}}\,y\rangle } {\displaystyle \langle \,a\,x{\text{,}}\,b\,y\rangle \,=\,a\,{\overline {b}}\,\langle x{\text{,}}\,y\rangle }. Here b ¯ {\displaystyle {\overline {b}}} {\displaystyle {\overline {b}}} is the complex conjugate of b {\displaystyle b} {\displaystyle b}.

It is also common to choose a non-unit vector q ∈ V {\displaystyle q\in V} {\displaystyle q\in V}, and normalize it directly in the Householder operator's expression:[4]

H q ( x ) = x − 2 ⟨ x , q ⟩ ⟨ q , q ⟩ q . {\displaystyle H_{q}\left(x\right)=x-2\,{\frac {\langle x,q\rangle }{\langle q,q\rangle }}\,q\,.} {\displaystyle H_{q}\left(x\right)=x-2\,{\frac {\langle x,q\rangle }{\langle q,q\rangle }}\,q\,.}

Such an operator is linear and self-adjoint.

If V = C n {\displaystyle V=\mathbb {C} ^{n}} {\displaystyle V=\mathbb {C} ^{n}}, note that the reflection hyperplane can be defined by its normal vector, a unit vector v → ∈ V {\textstyle {\vec {v}}\in V} {\textstyle {\vec {v}}\in V} (a vector with length 1 {\textstyle 1} {\textstyle 1}) that is orthogonal to the hyperplane. The reflection of a point x {\textstyle x} {\textstyle x} about this hyperplane is the Householder transformation:

x → − 2 ⟨ x → , v → ⟩ v → = x → − 2 v → ( v → ∗ x → ) , {\displaystyle {\vec {x}}-2\langle {\vec {x}},{\vec {v}}\rangle {\vec {v}}={\vec {x}}-2{\vec {v}}\left({\vec {v}}^{*}{\vec {x}}\right),} {\displaystyle {\vec {x}}-2\langle {\vec {x}},{\vec {v}}\rangle {\vec {v}}={\vec {x}}-2{\vec {v}}\left({\vec {v}}^{*}{\vec {x}}\right),}

where x → {\displaystyle {\vec {x}}} {\displaystyle {\vec {x}}} is the vector from the origin to the point x {\displaystyle x} {\displaystyle x}, and v → ∗ {\textstyle {\vec {v}}^{*}} {\textstyle {\vec {v}}^{*}} is the conjugate transpose of v → {\textstyle {\vec {v}}} {\textstyle {\vec {v}}}.

The Householder transformation acting as a reflection of x {\displaystyle x} {\displaystyle x} about the hyperplane defined by v {\displaystyle v} {\displaystyle v}.

Householder matrix

The matrix constructed from this transformation can be expressed in terms of an outer product as:

P = I − 2 v → v → ∗ {\displaystyle P=I-2{\vec {v}}{\vec {v}}^{*}} {\displaystyle P=I-2{\vec {v}}{\vec {v}}^{*}}

is known as the Householder matrix, where I {\textstyle I} {\textstyle I} is the identity matrix.

Properties

The Householder matrix has the following properties:

  • it is Hermitian: P = P ∗ {\textstyle P=P^{*}} {\textstyle P=P^{*}},
  • it is unitary: P − 1 = P ∗ {\textstyle P^{-1}=P^{*}} {\textstyle P^{-1}=P^{*}} (via the Sherman-Morrison formula),
  • hence it is involutory: P = P − 1 {\textstyle P=P^{-1}} {\textstyle P=P^{-1}}.
  • A Householder matrix has eigenvalues ± 1 {\textstyle \pm 1} {\textstyle \pm 1}. To see this, notice that if x → {\textstyle {\vec {x}}} {\textstyle {\vec {x}}} is orthogonal to the vector v → {\textstyle {\vec {v}}} {\textstyle {\vec {v}}} which was used to create the reflector, then P v x → = ( I − 2 v → v → ∗ ) x → = x → − 2 ⟨ v → , x → ⟩ v → = x → {\textstyle P_{v}{\vec {x}}=(I-2{\vec {v}}{\vec {v}}^{*}){\vec {x}}={\vec {x}}-2\langle {\vec {v}},{\vec {x}}\rangle {\vec {v}}={\vec {x}}} {\textstyle P_{v}{\vec {x}}=(I-2{\vec {v}}{\vec {v}}^{*}){\vec {x}}={\vec {x}}-2\langle {\vec {v}},{\vec {x}}\rangle {\vec {v}}={\vec {x}}}, i.e., 1 {\textstyle 1} {\textstyle 1} is an eigenvalue of multiplicity n − 1 {\textstyle n-1} {\textstyle n-1}, since there are n − 1 {\textstyle n-1} {\textstyle n-1} independent vectors orthogonal to v → {\textstyle {\vec {v}}} {\textstyle {\vec {v}}}. Also, notice P v v → = ( I − 2 v → v → ∗ ) v → = v → − 2 ⟨ v → , v → ⟩ v → = − v → {\textstyle P_{v}{\vec {v}}=(I-2{\vec {v}}{\vec {v}}^{*}){\vec {v}}={\vec {v}}-2\langle {\vec {v}},{\vec {v}}\rangle {\vec {v}}=-{\vec {v}}} {\textstyle P_{v}{\vec {v}}=(I-2{\vec {v}}{\vec {v}}^{*}){\vec {v}}={\vec {v}}-2\langle {\vec {v}},{\vec {v}}\rangle {\vec {v}}=-{\vec {v}}} (since v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}} is by definition a unit vector), and so − 1 {\textstyle -1} {\textstyle -1} is an eigenvalue with multiplicity 1 {\textstyle 1} {\textstyle 1}.
  • The determinant of a Householder reflector is − 1 {\textstyle -1} {\textstyle -1}, since the determinant of a matrix is the product of its eigenvalues, in this case one of which is − 1 {\textstyle -1} {\textstyle -1} with the remainder being 1 {\textstyle 1} {\textstyle 1} (as in the previous point), or via the Matrix determinant lemma.

Example

Consider the normalization of a vector v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}} containing 1 {\displaystyle 1} {\displaystyle 1} in each entry,

v → = 1 2 [ 1 1 ] . {\displaystyle {\vec {v}}={\frac {1}{\sqrt {2}}}{\begin{bmatrix}1\\1\end{bmatrix}}.} {\displaystyle {\vec {v}}={\frac {1}{\sqrt {2}}}{\begin{bmatrix}1\\1\end{bmatrix}}.}

Then the Householder matrix corresponding to the vector v {\displaystyle v} {\displaystyle v} is

P v = [ 1 0 0 1 ] − 2 ( 1 2 [ 1 1 ] ) ( 1 2 [ 1 1 ] ) {\displaystyle P_{v}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-2\left({\frac {1}{\sqrt {2}}}{\begin{bmatrix}1\\1\end{bmatrix}}\right)\left({\frac {1}{\sqrt {2}}}{\begin{bmatrix}1&1\end{bmatrix}}\right)} {\displaystyle P_{v}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-2\left({\frac {1}{\sqrt {2}}}{\begin{bmatrix}1\\1\end{bmatrix}}\right)\left({\frac {1}{\sqrt {2}}}{\begin{bmatrix}1&1\end{bmatrix}}\right)}
= [ 1 0 0 1 ] − [ 1 1 ] [ 1 1 ] {\displaystyle \quad ={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-{\begin{bmatrix}1\\1\end{bmatrix}}{\begin{bmatrix}1&1\end{bmatrix}}} {\displaystyle \quad ={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-{\begin{bmatrix}1\\1\end{bmatrix}}{\begin{bmatrix}1&1\end{bmatrix}}}
= [ 1 0 0 1 ] − [ 1 1 1 1 ] {\displaystyle \quad ={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-{\begin{bmatrix}1&1\\1&1\end{bmatrix}}} {\displaystyle \quad ={\begin{bmatrix}1&0\\0&1\end{bmatrix}}-{\begin{bmatrix}1&1\\1&1\end{bmatrix}}}
= [ 0 − 1 − 1 0 ] . {\displaystyle \quad ={\begin{bmatrix}0&-1\\-1&0\end{bmatrix}}.} {\displaystyle \quad ={\begin{bmatrix}0&-1\\-1&0\end{bmatrix}}.}

Note that if we have another vector q → {\displaystyle {\vec {q}}} {\displaystyle {\vec {q}}} representing a coordinate in the 2D plane

q → = [ x y ] , {\displaystyle {\vec {q}}={\begin{bmatrix}x\\y\end{bmatrix}},} {\displaystyle {\vec {q}}={\begin{bmatrix}x\\y\end{bmatrix}},}

then in this case P v {\displaystyle P_{v}} {\displaystyle P_{v}} flips and negates the x {\displaystyle x} {\displaystyle x} and y {\displaystyle y} {\displaystyle y} coordinates, in other words we have

P v [ x y ] = [ − y − x ] , {\displaystyle P_{v}{\begin{bmatrix}x\\y\end{bmatrix}}={\begin{bmatrix}-y\\-x\end{bmatrix}},} {\displaystyle P_{v}{\begin{bmatrix}x\\y\end{bmatrix}}={\begin{bmatrix}-y\\-x\end{bmatrix}},}

which corresponds to reflecting the vector across the line y = − x {\displaystyle y=-x} {\displaystyle y=-x}, which our original vector v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}} is normal to.

Applications

Geometric optics

In geometric optics, specular reflection can be expressed in terms of the Householder matrix (see Specular reflection § Vector formulation).

Numerical linear algebra

Note that representing a Householder matrix requires only the entries of a single vector, not of an entire matrix (which in most algorithms is never explicitly formed), thereby minimizing the required storage and memory references needed to use them.

Further, multiplying a Householder matrix by a vector does not involve a full matrix-vector multiplication, but rather only one vector dot product, and then one axpy operation. This means its arithmetic complexity is of the same order of two low-level BLAS-1 operations. Therefore, Householder matrices are extremely arithmetically efficient.[5]

Finally, using ⋅ ^ {\displaystyle {\hat {\cdot }}} {\displaystyle {\hat {\cdot }}} to denote the computed value and ⋅ {\displaystyle \cdot } {\displaystyle \cdot } to denote the mathematically exact value, then for a given Householder matrix P {\displaystyle P} {\displaystyle P},

P b ^ = ( P + Δ P ) b {\displaystyle {\widehat {Pb}}=(P+\Delta P)b} {\displaystyle {\widehat {Pb}}=(P+\Delta P)b}

Where | | Δ P | | F ≤ γ n ~ := c n u 1 − c n u {\displaystyle \vert \vert \Delta P\vert \vert _{F}\leq {\tilde {\gamma _{n}}}:={\frac {cnu}{1-cnu}}} {\displaystyle \vert \vert \Delta P\vert \vert _{F}\leq {\tilde {\gamma _{n}}}:={\frac {cnu}{1-cnu}}} (where u {\displaystyle u} {\displaystyle u} is unit roundoff, n {\displaystyle n} {\displaystyle n} the size of the matrix P {\displaystyle P} {\displaystyle P}, and c {\displaystyle c} {\displaystyle c} some small constant). In other words, multiplications by Householder matrices are also extremely backwards stable.[6]

Since Householder transformations minimize storage, memory references, arithmetic complexity, and optimize numerical stability, they are widely used in numerical linear algebra, for example, to annihilate the entries below the main diagonal of a matrix,[7] to perform QR decompositions and in the first step of the QR algorithm. They are also widely used for transforming to a Hessenberg form. For symmetric or Hermitian matrices, the symmetry can be preserved, resulting in tridiagonalization.[8][9]

QR decomposition

Householder transformations can be used to calculate a QR decomposition. Consider a square matrix upper triangularized up to column i − 1 {\displaystyle i-1} {\displaystyle i-1}, then our goal is to construct such Householder matrices that act upon the principal submatrices of that matrix, which has the form

A ( i ) = [ a 11 a 12 ⋯ a 1 n 0 a 22 ⋯ a 1 n ⋮ ⋱ ⋮ 0 ⋯ 0 x 1 = a i i ⋯ a i n 0 ⋯ 0 ⋮ ⋮ 0 ⋯ 0 x n + 1 − i = a n i ⋯ a n n ] {\displaystyle A^{(i)}={\begin{bmatrix}a_{11}&a_{12}&\cdots &&&a_{1n}\\0&a_{22}&\cdots &&&a_{1n}\\\vdots &&\ddots &&&\vdots \\0&\cdots &0&x_{1}=a_{ii}&\cdots &a_{in}\\0&\cdots &0&\vdots &&\vdots \\0&\cdots &0&x_{n+1-i}=a_{ni}&\cdots &a_{nn}\end{bmatrix}}} {\displaystyle A^{(i)}={\begin{bmatrix}a_{11}&a_{12}&\cdots &&&a_{1n}\\0&a_{22}&\cdots &&&a_{1n}\\\vdots &&\ddots &&&\vdots \\0&\cdots &0&x_{1}=a_{ii}&\cdots &a_{in}\\0&\cdots &0&\vdots &&\vdots \\0&\cdots &0&x_{n+1-i}=a_{ni}&\cdots &a_{nn}\end{bmatrix}}}

The matrix A ( i ) {\displaystyle A^{(i)}} {\displaystyle A^{(i)}} has the block form

A ( i ) = [ T i − 1 U ∗ 0 P v ] {\displaystyle A^{(i)}={\begin{bmatrix}T_{i-1}^{U}&*\\0&P_{v}\end{bmatrix}}} {\displaystyle A^{(i)}={\begin{bmatrix}T_{i-1}^{U}&*\\0&P_{v}\end{bmatrix}}}.

Note that A ( n ) {\displaystyle A^{(n)}} {\displaystyle A^{(n)}} is the matrix R {\displaystyle R} {\displaystyle R} in the Q R {\displaystyle Q\,R} {\displaystyle Q\,R} decomposition. Here the submatrix T i − 1 U {\displaystyle T_{i-1}^{U}} {\displaystyle T_{i-1}^{U}} is an i-1 x i-1 square upper triangular matrix, 0 represents the n-i+1 x n-i+1 zero matrix, * represents an n-i+1 x n-i+1 matrix, and P v {\displaystyle P_{v}} {\displaystyle P_{v}} represents an n-i+1 x n-i+1 matrix which is to made upper triangular by working in its subspace without changing the block matrix form given above for A ( i ) {\displaystyle A^{(i)}} {\displaystyle A^{(i)}}.

If x 1 = a i i {\displaystyle x_{1}=a_{ii}} {\displaystyle x_{1}=a_{ii}} is zero and there is a non-zero element a j i {\displaystyle a_{ji}} {\displaystyle a_{ji}} with j>i, then the i-th and j-th rows may be interchanged by pre-multiplying by the row interchange matrix, which is unitary. If a j i = 0 {\displaystyle a_{ji}=0} {\displaystyle a_{ji}=0} for all j>=i, then proceed to the next column.

If x 1 = a i i {\displaystyle x_{1}=a_{ii}} {\displaystyle x_{1}=a_{ii}} is non-real complex given by x 1 = r 1 e i ϕ 1 {\displaystyle x_{1}=r_{1}e^{i\phi _{1}}} {\displaystyle x_{1}=r_{1}e^{i\phi _{1}}} for real r 1 , ϕ 1 {\displaystyle r_{1}{\text{,}}\,\phi _{1}} {\displaystyle r_{1}{\text{,}}\,\phi _{1}}, then this matrix may be pre-multiplied by the diagonal unitary matrix U with U i i = e − i ϕ 1 {\displaystyle U_{ii}=e^{-i\phi _{1}}} {\displaystyle U_{ii}=e^{-i\phi _{1}}} and with 1 for the other diagonal elements so as to make the new x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} non-zero real while preserving the block form. This ensures that ⟨ x , e → 1 ⟩ = ⟨ e → 1 , x ⟩ = x 1 {\displaystyle \langle \mathbf {x} ,{\vec {e}}_{1}\rangle \,=\,\langle {\vec {e}}_{1},\mathbf {x} \rangle \,=\,x_{1}} {\displaystyle \langle \mathbf {x} ,{\vec {e}}_{1}\rangle \,=\,\langle {\vec {e}}_{1},\mathbf {x} \rangle \,=\,x_{1}}. Here e → 1 {\displaystyle {\vec {e}}_{1}} {\displaystyle {\vec {e}}_{1}} is an n-dimensional vector with 1 in the i-th position and 0's elsewhere. (Note that we previously established that Householder transformations are unitary matrices, and since the multiplication of unitary matrices is itself a unitary matrix, this gives us the unitary matrix of the QR decomposition).

Let e 1 , e 2 , ... , e n {\displaystyle \mathbf {e} _{1}{\text{,}}\,\mathbf {e} _{2}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}} {\displaystyle \mathbf {e} _{1}{\text{,}}\,\mathbf {e} _{2}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}} be the basis vectors for the n-dimensional matrices A ( i ) {\displaystyle A^{(i)}} {\displaystyle A^{(i)}} and let e → 1 , e → 2 , ... , e → n − i + 1 {\displaystyle {\vec {e}}_{1}{\text{,}}\,{\vec {e}}_{2}{\text{,}}\,{\text{...}}{\text{,}}\,{\vec {e}}_{n-i+1}} {\displaystyle {\vec {e}}_{1}{\text{,}}\,{\vec {e}}_{2}{\text{,}}\,{\text{...}}{\text{,}}\,{\vec {e}}_{n-i+1}} be the basis vectors e i , e i + 1 , ... , e n {\displaystyle \mathbf {e} _{i}{\text{,}}\,\mathbf {e} _{i+1}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}} {\displaystyle \mathbf {e} _{i}{\text{,}}\,\mathbf {e} _{i+1}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}}, respectively.

If we can find a v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}} so that P v x → = α e → 1 {\displaystyle P_{v}{\vec {x}}=\alpha {\vec {e}}_{1}} {\displaystyle P_{v}{\vec {x}}=\alpha {\vec {e}}_{1}} we can extend the upper triangulization by one column. The vector v {\displaystyle \mathbf {v} } {\displaystyle \mathbf {v} } is to be in the subspace spanned by the vectors { e i , e i + 1 , ... , e n } {\displaystyle \{\mathbf {e} _{i}{\text{,}}\,\mathbf {e} _{i+1}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}\}} {\displaystyle \{\mathbf {e} _{i}{\text{,}}\,\mathbf {e} _{i+1}{\text{,}}\,{\text{...}}{\text{,}}\,\mathbf {e} _{n}\}} that the vector x {\displaystyle \mathbf {x} } {\displaystyle \mathbf {x} } is in. With x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} real, it will be seen that α {\displaystyle \alpha } {\displaystyle \alpha } is real too. Thinking geometrically, we are looking for a plane so that the reflection about this plane happens to land directly on the basis vector. In other words,

x → − 2 ⟨ x → , v → ⟩ v → = α e → 1 {\displaystyle {\vec {x}}-2\langle {\vec {x}},{\vec {v}}\rangle {\vec {v}}=\alpha {\vec {e}}_{1}} {\displaystyle {\vec {x}}-2\langle {\vec {x}},{\vec {v}}\rangle {\vec {v}}=\alpha {\vec {e}}_{1}} 1

for some constant α {\displaystyle \alpha } {\displaystyle \alpha }. However, for this to happen, we must have v → ∝ x → − α e → 1 . {\displaystyle {\vec {v}}\propto {\vec {x}}-\alpha {\vec {e}}_{1}{\text{.}}} {\displaystyle {\vec {v}}\propto {\vec {x}}-\alpha {\vec {e}}_{1}{\text{.}}} And since v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}} is a unit vector, this means that we must have

v → = ± x → − α e → 1 ‖ x → − α e → 1 ‖ 2 {\displaystyle {\vec {v}}=\pm {\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}} {\displaystyle {\vec {v}}=\pm {\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}} 2

We will find that there are two possible values for α {\displaystyle \alpha } {\displaystyle \alpha }. That value that makes ‖ x → − α e → 1 ‖ 2 {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}} {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}} the largest is the one that should be used for greatest accuracy.

Now if we apply equation (2) back into equation (1), we get x → − α e → 1 = 2 ⟨ x → , x → − α e → 1 ‖ x → − α e → 1 ‖ 2 ⟩ x → − α e → 1 ‖ x → − α e → 1 ‖ 2 {\displaystyle {\vec {x}}-\alpha {\vec {e}}_{1}=2\left\langle {\vec {x}},{\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}\right\rangle {\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}} {\displaystyle {\vec {x}}-\alpha {\vec {e}}_{1}=2\left\langle {\vec {x}},{\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}\right\rangle {\frac {{\vec {x}}-\alpha {\vec {e}}_{1}}{\|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}}}} Or, in other words, by comparing the scalars in front of the vector x → − α e → 1 {\displaystyle {\vec {x}}-\alpha {\vec {e}}_{1}} {\displaystyle {\vec {x}}-\alpha {\vec {e}}_{1}} we must have ‖ x → − α e → 1 ‖ 2 2 = 2 ⟨ x → , x → − α e 1 ⟩ . {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}^{2}=2\langle {\vec {x}},{\vec {x}}-\alpha e_{1}\rangle {\text{.}}} {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}^{2}=2\langle {\vec {x}},{\vec {x}}-\alpha e_{1}\rangle {\text{.}}} If x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} is real then alpha is also real and is obtained from the equation ‖ x → ‖ 2 2 − 2 α x 1 + α 2 = 2 ( ‖ x → ‖ 2 2 − α x 1 ) {\displaystyle \|{\vec {x}}\|_{2}^{2}-2\alpha x_{1}+\alpha ^{2}=2(\|{\vec {x}}\|_{2}^{2}-\alpha x_{1})} {\displaystyle \|{\vec {x}}\|_{2}^{2}-2\alpha x_{1}+\alpha ^{2}=2(\|{\vec {x}}\|_{2}^{2}-\alpha x_{1})} which means α = ± ‖ x → ‖ 2 {\displaystyle \alpha =\pm \|{\vec {x}}\|_{2}} {\displaystyle \alpha =\pm \|{\vec {x}}\|_{2}} If x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} is not real, then α {\displaystyle \alpha } {\displaystyle \alpha } is also not real and is obtained from the equation ‖ x → ‖ 2 2 − α x 1 ∗ − α ∗ x 1 + | α | 2 = 2 ( ‖ x → ‖ 2 2 − α ∗ x 1 ) {\displaystyle \|{\vec {x}}\|_{2}^{2}-\alpha \,x_{1}^{*}-\alpha ^{*}\,x_{1}+|\alpha |^{2}\,=\,2(\|{\vec {x}}\|_{2}^{2}-\alpha ^{*}\,x_{1})} {\displaystyle \|{\vec {x}}\|_{2}^{2}-\alpha \,x_{1}^{*}-\alpha ^{*}\,x_{1}+|\alpha |^{2}\,=\,2(\|{\vec {x}}\|_{2}^{2}-\alpha ^{*}\,x_{1})} Or, equivalently, | α | 2 = ‖ x → ‖ 2 2 + ( α x 1 ∗ − α ∗ x 1 ) {\displaystyle |\alpha |^{2}\,=\,\|{\vec {x}}\|_{2}^{2}+(\alpha \,x_{1}^{*}-\alpha ^{*}\,x_{1})} {\displaystyle |\alpha |^{2}\,=\,\|{\vec {x}}\|_{2}^{2}+(\alpha \,x_{1}^{*}-\alpha ^{*}\,x_{1})} The term in parentheses is pure imaginary, This equation requires that | α | = ‖ x → ‖ 2 {\displaystyle |\alpha |=\|{\vec {x}}\|_{2}} {\displaystyle |\alpha |=\|{\vec {x}}\|_{2}} and that arg ⁡ ( α ) = arg ⁡ ( x 1 ) {\displaystyle \arg(\alpha )\,=\,\arg(x_{1})} {\displaystyle \arg(\alpha )\,=\,\arg(x_{1})} within a multiple of π {\displaystyle \pi } {\displaystyle \pi }.

This completes the construction; however, in practice we want to avoid catastrophic cancellation in equation (2). To do so for real x 1 {\displaystyle x_{1}} {\displaystyle x_{1}}, we choose[5] the sign of α {\displaystyle \alpha } {\displaystyle \alpha } as α = − sgn ⁡ ( R e ( x 1 ) ) ‖ x → ‖ 2 {\displaystyle \alpha =-\operatorname {sgn}(\mathrm {Re} (x_{1}))\|{\vec {x}}\|_{2}} {\displaystyle \alpha =-\operatorname {sgn} (\mathrm {Re} (x_{1}))\|{\vec {x}}\|_{2}} and for complex x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} we choose the sign for α {\displaystyle \alpha } {\displaystyle \alpha } as α = − e i arg ⁡ ( x 1 ) ‖ x ‖ 2 {\displaystyle \alpha =-e^{i\,\arg(x_{1})}\,\|x\|_{2}} {\displaystyle \alpha =-e^{i\,\arg(x_{1})}\,\|x\|_{2}}, where i {\displaystyle i} {\displaystyle i} is the square root of − 1 {\displaystyle -1} {\displaystyle -1}. This agrees with the equation just given for α {\displaystyle \alpha } {\displaystyle \alpha } when x 1 {\displaystyle x_{1}} {\displaystyle x_{1}} is real. These choices of sign make ‖ x → − α e → 1 ‖ 2 {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}} {\displaystyle \|{\vec {x}}-\alpha {\vec {e}}_{1}\|_{2}} the largest. When | x 1 | ≪ ‖ x ‖ 2 {\displaystyle |x_{1}|\ll \|x\|_{2}} {\displaystyle |x_{1}|\ll \|x\|_{2}}, it makes little difference which sign is chosen.

Householder transformations can similarly be applied to a non-square complex matrix. The decompositions are somewhat different.

Tridiagonalization (Hessenberg)

Real symmetric matrix

This procedure is presented in Numerical Analysis by Burden and Faires for real symmetric matrices. In the non-symmetric case, it is still useful as a similar procedure can result in a Hessenberg matrix.

It uses a slightly altered sgn {\displaystyle \operatorname {sgn} } {\displaystyle \operatorname {sgn} } function with sgn ⁡ ( 0 ) = 1 {\displaystyle \operatorname {sgn} (0)=1} {\displaystyle \operatorname {sgn} (0)=1}.[10] In the first step, to form the Householder matrix in each step we need to determine α {\textstyle \alpha } {\textstyle \alpha } and r {\textstyle r} {\textstyle r}, which are:

α = − sgn ⁡ ( a 21 ) ∑ j = 2 n a j 1 2 ; r = 1 2 ( α 2 − a 21 α ) ; {\displaystyle {\begin{aligned}\alpha &=-\operatorname {sgn} \left(a_{21}\right){\sqrt {\sum _{j=2}^{n}a_{j1}^{2}}};\\r&={\sqrt {{\frac {1}{2}}\left(\alpha ^{2}-a_{21}\alpha \right)}};\end{aligned}}} {\displaystyle {\begin{aligned}\alpha &=-\operatorname {sgn} \left(a_{21}\right){\sqrt {\sum _{j=2}^{n}a_{j1}^{2}}};\\r&={\sqrt {{\frac {1}{2}}\left(\alpha ^{2}-a_{21}\alpha \right)}};\end{aligned}}}

From α {\textstyle \alpha } {\textstyle \alpha } and r {\textstyle r} {\textstyle r}, construct vector v {\textstyle v} {\textstyle v}:

v → ( 1 ) = [ v 1 v 2 ⋮ v n ] , {\displaystyle {\vec {v}}^{(1)}={\begin{bmatrix}v_{1}\\v_{2}\\\vdots \\v_{n}\end{bmatrix}},} {\displaystyle {\vec {v}}^{(1)}={\begin{bmatrix}v_{1}\\v_{2}\\\vdots \\v_{n}\end{bmatrix}},}

where v 1 = 0 {\textstyle v_{1}=0} {\textstyle v_{1}=0}, v 2 = a 21 − α 2 r {\textstyle v_{2}={\frac {a_{21}-\alpha }{2r}}} {\textstyle v_{2}={\frac {a_{21}-\alpha }{2r}}}, and

v k = a k 1 2 r {\displaystyle v_{k}={\frac {a_{k1}}{2r}}} {\displaystyle v_{k}={\frac {a_{k1}}{2r}}} for each k = 3 , 4 … n {\displaystyle k=3,4\ldots n} {\displaystyle k=3,4\ldots n}

Then compute:

P 1 = I − 2 v → ( 1 ) ( v → ( 1 ) ) T A ( 2 ) = P 1 A P 1 {\displaystyle {\begin{aligned}P^{1}&=I-2{\vec {v}}^{(1)}\left({\vec {v}}^{(1)}\right)^{\textsf {T}}\\A^{(2)}&=P^{1}AP^{1}\end{aligned}}} {\displaystyle {\begin{aligned}P^{1}&=I-2{\vec {v}}^{(1)}\left({\vec {v}}^{(1)}\right)^{\textsf {T}}\\A^{(2)}&=P^{1}AP^{1}\end{aligned}}}

Since P 1 {\displaystyle P^{1}} {\displaystyle P^{1}} is its own inverse (involutory), this is a similarity transformation. Having found P 1 {\textstyle P^{1}} {\textstyle P^{1}} and computed A ( 2 ) {\textstyle A^{(2)}} {\textstyle A^{(2)}} the process is repeated for k = 2 , 3 , … , n − 2 {\textstyle k=2,3,\ldots ,n-2} {\textstyle k=2,3,\ldots ,n-2} as follows:

α = − sgn ⁡ ( a k + 1 , k k ) ∑ j = k + 1 n ( a j k k ) 2 r = 1 2 ( α 2 − a k + 1 , k k α ) v 1 k = v 2 k = ⋯ = v k k = 0 v k + 1 k = a k + 1 , k k − α 2 r v j k = a j k k 2 r  for  j = k + 2 ,   k + 3 ,   … ,   n P k = I − 2 v → ( k ) ( v → ( k ) ) T A ( k + 1 ) = P k A ( k ) P k {\displaystyle {\begin{aligned}\alpha &=-\operatorname {sgn} \left(a_{k+1,k}^{k}\right){\sqrt {\sum _{j=k+1}^{n}\left(a_{jk}^{k}\right)^{2}}}\\[2pt]r&={\sqrt {{\frac {1}{2}}\left(\alpha ^{2}-a_{k+1,k}^{k}\alpha \right)}}\\[2pt]v_{1}^{k}&=v_{2}^{k}=\cdots =v_{k}^{k}=0\\[2pt]v_{k+1}^{k}&={\frac {a_{k+1,k}^{k}-\alpha }{2r}}\\v_{j}^{k}&={\frac {a_{jk}^{k}}{2r}}{\text{ for }}j=k+2,\ k+3,\ \ldots ,\ n\\P^{k}&=I-2{\vec {v}}^{(k)}\left({\vec {v}}^{(k)}\right)^{\textsf {T}}\\A^{(k+1)}&=P^{k}A^{(k)}P^{k}\end{aligned}}} {\displaystyle {\begin{aligned}\alpha &=-\operatorname {sgn} \left(a_{k+1,k}^{k}\right){\sqrt {\sum _{j=k+1}^{n}\left(a_{jk}^{k}\right)^{2}}}\\[2pt]r&={\sqrt {{\frac {1}{2}}\left(\alpha ^{2}-a_{k+1,k}^{k}\alpha \right)}}\\[2pt]v_{1}^{k}&=v_{2}^{k}=\cdots =v_{k}^{k}=0\\[2pt]v_{k+1}^{k}&={\frac {a_{k+1,k}^{k}-\alpha }{2r}}\\v_{j}^{k}&={\frac {a_{jk}^{k}}{2r}}{\text{ for }}j=k+2,\ k+3,\ \ldots ,\ n\\P^{k}&=I-2{\vec {v}}^{(k)}\left({\vec {v}}^{(k)}\right)^{\textsf {T}}\\A^{(k+1)}&=P^{k}A^{(k)}P^{k}\end{aligned}}}

Continuing in this manner, the tridiagonal and real symmetric matrix is formed. Since each Householder reflection is a similarity transformation, so is the overall transformation.

Hermitian matrix

Householder transformations can also transform a Hermitian matrix A {\displaystyle A} {\displaystyle A} into a Hermitian tri-diagonal matrix. We follow the approach in the University of Connecticut paper[11] except replace matrix transposes by matrix Hermitian conjugates for tri-diagonalization of complex Hermitian matrices rather than real symmetric matrices.

To start, let

A = ( a 11 a → 1 ∗ a → 1 A 1 ) {\displaystyle A=\left({\begin{array}{c|c}a_{11}&{\vec {a}}_{1}^{*}\\\hline {\vec {a}}_{1}&A_{1}\end{array}}\right)} {\displaystyle A=\left({\begin{array}{c|c}a_{11}&{\vec {a}}_{1}^{*}\\\hline {\vec {a}}_{1}&A_{1}\end{array}}\right)}

be an n × n {\displaystyle n\times n} {\displaystyle n\times n} Hermitian matrix. Here a 11 {\displaystyle a_{11}} {\displaystyle a_{11}} is a Hermitean 1 × 1 {\displaystyle 1\times 1} {\displaystyle 1\times 1} sub-matrix so is a real number, a → 1 {\displaystyle {\vec {a}}_{1}} {\displaystyle {\vec {a}}_{1}} is an ( n − 1 ) × 1 {\displaystyle (n-1)\times 1} {\displaystyle (n-1)\times 1} column vector, a → 1 ∗ {\displaystyle {\vec {a}}_{1}^{*}} {\displaystyle {\vec {a}}_{1}^{*}} is its Hermitian conjugate, which is a 1 × ( n − 1 ) {\displaystyle 1\times (n-1)} {\displaystyle 1\times (n-1)} row vector, and A 1 {\displaystyle A_{1}} {\displaystyle A_{1}} is an ( n − 1 ) × ( n − 1 ) {\displaystyle (n-1)\times (n-1)} {\displaystyle (n-1)\times (n-1)} Hermitian sub-matrix. Let the Householder transformation matrix P 1 {\displaystyle P_{1}} {\displaystyle P_{1}}have the form

P 1 = ( 1 0 → ∗ 0 → H 1 ) {\displaystyle P_{1}=\left({\begin{array}{c|c}1&{\vec {0}}^{*}\\\hline {\vec {0}}&H_{1}\end{array}}\right)} {\displaystyle P_{1}=\left({\begin{array}{c|c}1&{\vec {0}}^{*}\\\hline {\vec {0}}&H_{1}\end{array}}\right)}

where 0 → {\displaystyle {\vec {0}}} {\displaystyle {\vec {0}}} is an 1 × ( n − 1 ) {\displaystyle 1\times (n-1)} {\displaystyle 1\times (n-1)} column vector of zeros, 0 → ∗ {\displaystyle {\vec {0}}^{*}} {\displaystyle {\vec {0}}^{*}} is its Hermitian conjugate and so is an ( n − 1 ) × 1 {\displaystyle (n-1)\times 1} {\displaystyle (n-1)\times 1} row vector of zeros, and where H 1 = I − 2 v → v → ∗ {\displaystyle H_{1}=I\,-\,2\,{\vec {v}}\,{\vec {v}}^{*}} {\displaystyle H_{1}=I\,-\,2\,{\vec {v}}\,{\vec {v}}^{*}} and I {\displaystyle I} {\displaystyle I} is the ( n − 1 ) × ( n − 1 ) {\displaystyle (n-1)\times (n-1)} {\displaystyle (n-1)\times (n-1)} identity matrix. As before, v → ∗ {\displaystyle {\vec {v}}^{*}} {\displaystyle {\vec {v}}^{*}} is the matrix transpose of the complex conjugate or Hermitian conjugate of the column vector v → {\displaystyle {\vec {v}}} {\displaystyle {\vec {v}}}. Also, as before, the Householder transformation H 1 {\displaystyle H_{1}} {\displaystyle H_{1}} is a unitary hermitian involutory matrix. Then

P 1 A P 1 ∗ = ( a 11 ( H 1 a 1 ) ∗ H 1 a 1 H 1 A 1 H 1 ∗ ) {\displaystyle P_{1}\,A\,P_{1}^{*}=\,\left({\begin{array}{c|c}a_{11}&(H_{1}\,a_{1})^{*}\\\hline H_{1}\,a_{1}&H_{1}\,A_{1}\,H_{1}^{*}\end{array}}\right)} {\displaystyle P_{1}\,A\,P_{1}^{*}=\,\left({\begin{array}{c|c}a_{11}&(H_{1}\,a_{1})^{*}\\\hline H_{1}\,a_{1}&H_{1}\,A_{1}\,H_{1}^{*}\end{array}}\right)}

Since H 1 {\displaystyle H_{1}} {\displaystyle H_{1}} and therefore P 1 {\displaystyle P_{1}} {\displaystyle P_{1}} is hermitian and involutory H 1 = H 1 ∗ = H 1 − 1 P 1 = P 1 ∗ = P 1 − 1 {\displaystyle {\begin{aligned}H_{1}&=H_{1}^{*}=H_{1}^{-1}\\P_{1}&=P_{1}^{*}=P_{1}^{-1}\end{aligned}}} {\displaystyle {\begin{aligned}H_{1}&=H_{1}^{*}=H_{1}^{-1}\\P_{1}&=P_{1}^{*}=P_{1}^{-1}\end{aligned}}} and so this is a unitary similarity matrix transformation.

If we have H 1 a → 1 = α 1 e → 1 {\displaystyle H_{1}{\vec {a}}_{1}\,=\,\alpha _{1}\,{\vec {e}}_{1}} {\displaystyle H_{1}{\vec {a}}_{1}\,=\,\alpha _{1}\,{\vec {e}}_{1}}

then P 1 A P 1 ∗ = ( a 11 α 1 ∗ 0 → ∗ α 1 a 22 ( 1 ) a → 2 ∗ 0 → a → 2 A 2 ) {\displaystyle P_{1}\,A\,P_{1}^{*}\,=\,\left({\begin{array}{c|c|c}a_{11}&\alpha _{1}^{*}&{\vec {0}}^{*}\\\hline \alpha _{1}&a_{22}^{(1)}&{\vec {a}}_{2}^{*}\\\hline {\vec {0}}&{\vec {a}}_{2}&A_{2}\end{array}}\right)} {\displaystyle P_{1}\,A\,P_{1}^{*}\,=\,\left({\begin{array}{c|c|c}a_{11}&\alpha _{1}^{*}&{\vec {0}}^{*}\\\hline \alpha _{1}&a_{22}^{(1)}&{\vec {a}}_{2}^{*}\\\hline {\vec {0}}&{\vec {a}}_{2}&A_{2}\end{array}}\right)} where the ( n − 2 ) × ( n − 2 ) {\displaystyle (n-2)\times (n-2)} {\displaystyle (n-2)\times (n-2)} matrix A 2 {\displaystyle A_{2}} {\displaystyle A_{2}} is Hermitian since a unitary similarity transformation of a Hermitian matrix is again Hermitian. Now 0 → {\displaystyle {\vec {0}}} {\displaystyle {\vec {0}}} is an ( n − 2 ) × 1 {\displaystyle (n-2)\times 1} {\displaystyle (n-2)\times 1} column vector of zeros and 0 → ∗ {\displaystyle {\vec {0}}^{*}} {\displaystyle {\vec {0}}^{*}} is its Hermitian conjugate and is a 1 × ( n − 1 ) {\displaystyle 1\times (n-1)} {\displaystyle 1\times (n-1)} row vector of zeros. And a → 2 {\displaystyle {\vec {a}}_{2}} {\displaystyle {\vec {a}}_{2}} is a ( n − 2 ) × 1 {\displaystyle (n-2)\times 1} {\displaystyle (n-2)\times 1} column vector and a → 2 ∗ {\displaystyle {\vec {a}}_{2}^{*}} {\displaystyle {\vec {a}}_{2}^{*}} is its Hermitian conjugate and is a 1 × ( n − 2 ) {\displaystyle 1\times (n-2)} {\displaystyle 1\times (n-2)} row vector. We already know how to find the Householder transformation matrices H i {\displaystyle H_{i}} {\displaystyle H_{i}}.

Repeat this process for a total of n − 2 {\displaystyle n-2} {\displaystyle n-2} times for n > 2 {\displaystyle n>2} {\displaystyle n>2} to obtain the tri-diagonalization by a succession of unitary Hermitian involutory similarity transformations. For instance, the following example for a 4 × 4 {\displaystyle 4\times 4} {\displaystyle 4\times 4} real symmetric matrix requires two Householder transformation steps to be transformed into a tri-diagonal real symmetric matrix, a 3 × 3 {\displaystyle 3\times 3} {\displaystyle 3\times 3} real symmetric or Hermitian matrix requires only one step, and a 1 × 1 {\displaystyle 1\times 1} {\displaystyle 1\times 1} real matrix or a 2 × 2 {\displaystyle 2\times 2} {\displaystyle 2\times 2} real symmetric or Hermitian matrix is already so tri-diagonalized.

Since the product of unitary matrices is again a unitary matrix and since the similarity transformation of a similarity transformation is again a similarity transformation, the overall transformation is a unitary similarity transformation.

Examples

In this example, also from Burden and Faires,[10] the given matrix is transformed to the similar tridiagonal matrix A3 by using the Householder method.

A = [ 4 1 − 2 2 1 2 0 1 − 2 0 3 − 2 2 1 − 2 − 1 ] , {\displaystyle \mathbf {A} ={\begin{bmatrix}4&1&-2&2\\1&2&0&1\\-2&0&3&-2\\2&1&-2&-1\end{bmatrix}},} {\displaystyle \mathbf {A} ={\begin{bmatrix}4&1&-2&2\\1&2&0&1\\-2&0&3&-2\\2&1&-2&-1\end{bmatrix}},}

Following those steps in the Householder method, we have:

The first Householder matrix:

Q 1 = [ 1 0 0 0 0 − 1 3 2 3 − 2 3 0 2 3 2 3 1 3 0 − 2 3 1 3 2 3 ] , A 2 = Q 1 A Q 1 = [ 4 − 3 0 0 − 3 10 3 1 4 3 0 1 5 3 − 4 3 0 4 3 − 4 3 − 1 ] , {\displaystyle {\begin{aligned}Q_{1}&={\begin{bmatrix}1&0&0&0\\0&-{\frac {1}{3}}&{\frac {2}{3}}&-{\frac {2}{3}}\\0&{\frac {2}{3}}&{\frac {2}{3}}&{\frac {1}{3}}\\0&-{\frac {2}{3}}&{\frac {1}{3}}&{\frac {2}{3}}\end{bmatrix}},\\A_{2}=Q_{1}AQ_{1}&={\begin{bmatrix}4&-3&0&0\\-3&{\frac {10}{3}}&1&{\frac {4}{3}}\\0&1&{\frac {5}{3}}&-{\frac {4}{3}}\\0&{\frac {4}{3}}&-{\frac {4}{3}}&-1\end{bmatrix}},\end{aligned}}} {\displaystyle {\begin{aligned}Q_{1}&={\begin{bmatrix}1&0&0&0\\0&-{\frac {1}{3}}&{\frac {2}{3}}&-{\frac {2}{3}}\\0&{\frac {2}{3}}&{\frac {2}{3}}&{\frac {1}{3}}\\0&-{\frac {2}{3}}&{\frac {1}{3}}&{\frac {2}{3}}\end{bmatrix}},\\A_{2}=Q_{1}AQ_{1}&={\begin{bmatrix}4&-3&0&0\\-3&{\frac {10}{3}}&1&{\frac {4}{3}}\\0&1&{\frac {5}{3}}&-{\frac {4}{3}}\\0&{\frac {4}{3}}&-{\frac {4}{3}}&-1\end{bmatrix}},\end{aligned}}}

Used A 2 {\textstyle A_{2}} {\textstyle A_{2}} to form

Q 2 = [ 1 0 0 0 0 1 0 0 0 0 − 3 5 − 4 5 0 0 − 4 5 3 5 ] , A 3 = Q 2 A 2 Q 2 = [ 4 − 3 0 0 − 3 10 3 − 5 3 0 0 − 5 3 − 33 25 68 75 0 0 68 75 149 75 ] , {\displaystyle {\begin{aligned}Q_{2}&={\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&-{\frac {3}{5}}&-{\frac {4}{5}}\\0&0&-{\frac {4}{5}}&{\frac {3}{5}}\end{bmatrix}},\\A_{3}=Q_{2}A_{2}Q_{2}&={\begin{bmatrix}4&-3&0&0\\-3&{\frac {10}{3}}&-{\frac {5}{3}}&0\\0&-{\frac {5}{3}}&-{\frac {33}{25}}&{\frac {68}{75}}\\0&0&{\frac {68}{75}}&{\frac {149}{75}}\end{bmatrix}},\end{aligned}}} {\displaystyle {\begin{aligned}Q_{2}&={\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&-{\frac {3}{5}}&-{\frac {4}{5}}\\0&0&-{\frac {4}{5}}&{\frac {3}{5}}\end{bmatrix}},\\A_{3}=Q_{2}A_{2}Q_{2}&={\begin{bmatrix}4&-3&0&0\\-3&{\frac {10}{3}}&-{\frac {5}{3}}&0\\0&-{\frac {5}{3}}&-{\frac {33}{25}}&{\frac {68}{75}}\\0&0&{\frac {68}{75}}&{\frac {149}{75}}\end{bmatrix}},\end{aligned}}}

As we can see, the final result is a tridiagonal symmetric matrix which is similar to the original one. The process is finished after two steps.

Quantum computation

Picture showing the geometric interpretation of the first iteration of Grover's algorithm. The state vector | s ⟩ {\displaystyle |s\rangle } {\displaystyle |s\rangle } is rotated towards the target vector | ω ⟩ {\displaystyle |\omega \rangle } {\displaystyle |\omega \rangle } as shown.

As unitary matrices are useful in quantum computation, and Householder transformations are unitary, they are very useful in quantum computing. One of the central algorithms where they're useful is Grover's algorithm, where we are trying to solve for a representation of an oracle function represented by what turns out to be a Householder transformation:

{ U ω | x ⟩ = − | x ⟩ for  x = ω , that is,  f ( x ) = 1 , U ω | x ⟩ = | x ⟩ for  x ≠ ω , that is,  f ( x ) = 0. {\displaystyle {\begin{cases}U_{\omega }|x\rangle =-|x\rangle &{\text{for }}x=\omega {\text{, that is, }}f(x)=1,\\U_{\omega }|x\rangle =|x\rangle &{\text{for }}x\neq \omega {\text{, that is, }}f(x)=0.\end{cases}}} {\displaystyle {\begin{cases}U_{\omega }|x\rangle =-|x\rangle &{\text{for }}x=\omega {\text{, that is, }}f(x)=1,\\U_{\omega }|x\rangle =|x\rangle &{\text{for }}x\neq \omega {\text{, that is, }}f(x)=0.\end{cases}}}

(here the | x ⟩ {\displaystyle |x\rangle } {\displaystyle |x\rangle } is part of the bra-ket notation and is analogous to x → {\displaystyle {\vec {x}}} {\displaystyle {\vec {x}}} which we were using previously)

This is done via an algorithm that iterates via the oracle function U ω {\displaystyle U_{\omega }} {\displaystyle U_{\omega }} and another operator U s {\displaystyle U_{s}} {\displaystyle U_{s}} known as the Grover diffusion operator defined by

| s ⟩ = 1 N ∑ x = 0 N − 1 | x ⟩ . {\displaystyle |s\rangle ={\frac {1}{\sqrt {N}}}\sum _{x=0}^{N-1}|x\rangle .} {\displaystyle |s\rangle ={\frac {1}{\sqrt {N}}}\sum _{x=0}^{N-1}|x\rangle .} and U s = 2 | s ⟩ ⟨ s | − I {\displaystyle U_{s}=2\left|s\right\rangle \!\!\left\langle s\right|-I} {\displaystyle U_{s}=2\left|s\right\rangle \!\!\left\langle s\right|-I}.

Computational and theoretical relationship to other unitary transformations

The Householder transformation is a reflection about a hyperplane with unit normal vector v {\textstyle v} {\textstyle v}, as stated earlier. An N {\textstyle N} {\textstyle N}-by- N {\textstyle N} {\textstyle N} unitary transformation U {\textstyle U} {\textstyle U} satisfies U U ∗ = I {\textstyle UU^{*}=I} {\textstyle UU^{*}=I}. Taking the determinant ( N {\textstyle N} {\textstyle N}-th power of the geometric mean) and trace (proportional to arithmetic mean) of a unitary matrix reveals that its eigenvalues λ i {\textstyle \lambda _{i}} {\textstyle \lambda _{i}} have unit modulus. This can be seen directly and swiftly:

Trace ⁡ ( U U ∗ ) N = ∑ j = 1 N | λ j | 2 N = 1 , det ⁡ ( U U ∗ ) = ∏ j = 1 N | λ j | 2 = 1. {\displaystyle {\begin{aligned}{\frac {\operatorname {Trace} \left(UU^{*}\right)}{N}}&={\frac {\sum _{j=1}^{N}\left|\lambda _{j}\right|^{2}}{N}}=1,&\operatorname {det} \left(UU^{*}\right)&=\prod _{j=1}^{N}\left|\lambda _{j}\right|^{2}=1.\end{aligned}}} {\displaystyle {\begin{aligned}{\frac {\operatorname {Trace} \left(UU^{*}\right)}{N}}&={\frac {\sum _{j=1}^{N}\left|\lambda _{j}\right|^{2}}{N}}=1,&\operatorname {det} \left(UU^{*}\right)&=\prod _{j=1}^{N}\left|\lambda _{j}\right|^{2}=1.\end{aligned}}}

Since arithmetic and geometric means are equal if the variables are constant (see inequality of arithmetic and geometric means), we establish the claim of unit modulus.

For the case of real valued unitary matrices we obtain orthogonal matrices, U U T = I {\textstyle UU^{\textsf {T}}=I} {\textstyle UU^{\textsf {T}}=I}. It follows rather readily (see Orthogonal matrix) that any orthogonal matrix can be decomposed into a product of 2-by-2 rotations, called Givens rotations, and Householder reflections. This is appealing intuitively since multiplication of a vector by an orthogonal matrix preserves the length of that vector, and rotations and reflections exhaust the set of (real valued) geometric operations that render invariant a vector's length.

The Householder transformation was shown to have a one-to-one relationship with the canonical coset decomposition of unitary matrices defined in group theory, which can be used to parametrize unitary operators in a very efficient manner.[12]

Finally we note that a single Householder transform, unlike a solitary Givens transform, can act on all columns of a matrix, and as such exhibits the lowest computational cost for QR decomposition and tridiagonalization. The penalty for this "computational optimality" is, of course, that Householder operations cannot be as deeply or efficiently parallelized. As such Householder is preferred for dense matrices on sequential machines, whilst Givens is preferred on sparse matrices, and/or parallel machines.

See also

Notes

  1. Householder, A. S. (1958). "Unitary Triangularization of a Nonsymmetric Matrix" (PDF). Journal of the ACM. 5 (4): 339–342. doi:10.1145/320941.320947. MR 0111128. S2CID 9858625.
  2. Roman 2008, p. 243-244
  3. Methods of Applied Mathematics for Engineers and Scientist. Cambridge University Press. 28 June 2013. pp. Section E.4.11. ISBN 9781107244467.
  4. Roman 2008, p. 244
  5. Saad, Yousef (2003). Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics. pp. 11–14.
  6. Higham, Nicholas J. (2002). Accuracy and stability of numerical algorithms (2nd ed.). Philadelphia: Society for Industrial and Applied Mathematics. p. 358. ISBN 0-89871-521-0.
  7. Taboga, Marco. "Householder matrix, Lectures on matrix algebra".
  8. Schabauer, Hannes; Pacher, Christoph; Sunderland, Andrew G.; Gansterer, Wilfried N. (2010-05-01). "Toward a parallel solver for generalized complex symmetric eigenvalue problems". Procedia Computer Science. 1 (1): 437–445. doi:10.1016/j.procs.2010.04.047.
  9. Golub, Gene Howard; Van Loan, Charles F. (1996). Matrix computations (3rd ed.). Baltimore London: Johns Hopkins university press. p. 211. ISBN 0-8018-5414-8.
  10. Burden, Richard; Faires, Douglas; Burden, Annette (2016). Numerical analysis (10th ed.). Thomson Brooks/Cole. ISBN 9781305253667.
  11. Rozman. "Tridiagonalization" (PDF).
  12. Renan Cabrera; Traci Strohecker; Herschel Rabitz (2010). "The canonical coset decomposition of unitary matrices through Householder transformations". Journal of Mathematical Physics. 51 (8): 082101. arXiv:1008.2477. Bibcode:2010JMP....51h2101C. doi:10.1063/1.3466798. S2CID 119641896.

References