# Decomposition of a unitary matrix via Householder matrices

If $U$ is unitary, how can I show that there exist $w_{1},w_{2},...,w_{k}\in \mathbb{C}^{n}$, $k\leq n$, and $\theta_{1},\theta_{2},...,\theta_{n}\in \mathbb{R}$ such that

$U=U_{w_{1}}U_{w_{2}}\cdots U_{w_{k}} \begin{pmatrix} e^{i\theta_{1}} & 0 & \cdots & 0 \\ 0 & e^{i\theta_{2}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & e^{i\theta_{n}} \end{pmatrix}$

where $U_{w_{i}}=I-\frac{2w_{i}w_{i}^{*}}{w_{i}^{*}w_{i}}$ are householder matrices.

-
Consider what happens if you attempt to perform QR decomposition on $\mathbf U$... – J. M. Sep 24 '11 at 6:34
@J.M. I know the matrix $diag(e^{\theta_{1}}\cdots e^{\theta_{n}})$ is unitary because I can show that, but is there some intuition as to how it acts geometrically? – Edison Sep 24 '11 at 14:34
In the real case, it corresponds to a diagonal matrix with entries of $\pm 1$, which of course corresponds to reflections. I don't know of a correspondingly intuitive interpretation for the "phase matrix". Sorry. :( – J. M. Sep 24 '11 at 14:48
You really haven't tried actually applying Householder reflectors to some unitary matrix and observing what happens? – J. M. Sep 28 '11 at 9:31
@J.M., I am not familiar with the term "Householder reflector" but I assume that refers to the QR decomposition you mentioned in an earlier comment? I actually do not know what that is, but this is something I will investigate. – Edison Sep 28 '11 at 12:53

a) Householder matrices are unitary; b) you can choose them to make the right-hand matrix upper-triangular; and c) a unitary upper-triangular matrix is diagonal. In more detail:

a) Householder matrices are unitary:

\begin{align} \left(I-2\frac{ww^\dagger}{w^\dagger w}\right)\left(I-2\frac{ww^\dagger}{w^\dagger w}\right)^\dagger & =\left(I-2\frac{ww^\dagger}{w^\dagger w}\right)\left(I-2\frac{ww^\dagger}{w^\dagger w}\right) \\ & =I-4\frac{ww^\dagger}{w^\dagger w}+4\frac{ww^\dagger ww^\dagger}{w^\dagger ww^\dagger w} \\ & =I-4\frac{ww^\dagger}{w^\dagger w}+4\frac{ww^\dagger}{w^\dagger w} \\ & =I\;. \end{align}

b) You can choose them to make the right-hand matrix upper-triangular:

It suffices to show that you can produce zeros below the diagonal in the first column; the remaining columns can be dealt with in the same way, with the identity in the rows already fixed to preserve the zeros there. So let $v$ be the first column of $U$, and choose $w=v+\lambda e_1$, where $e_1$ is the unit column vector with a $1$ in the first row. Then

\begin{align} \left(I-2\frac{ww^\dagger}{w^\dagger w}\right)v & =\left(I-2\frac{(v+\lambda e_1)(v+\lambda e_1)^\dagger}{(v+\lambda e_1)^\dagger(v+\lambda e_1)}\right)v \\ &=v-2\frac{v^\dagger v+\lambda v_1}{v^\dagger v+2\lambda v_1+\lambda^2}(v+\lambda e_1)\;. \end{align}

The fraction is $1$ for $\lambda=0$ and $0$ for $\lambda\to\infty$. Thus it takes the value $\frac12$ for some value of $\lambda$, and for that value, the resulting column vector is proportional to $e_1$, that is, it has zeros below the diagonal.

c) A unitary upper-triangular matrix is diagonal:

This follows directly from the fact that the scalar products of columns of a unitary matrix vanish.

The matrix obtained by left-multiplying $U$ by Householder matrices is upper-triangular; it is a product of unitary matrices, and thus itself unitary; thus it is a diagonal unitary matrix; and such a matrix must have the form of the right-most factor in the question.

-
...and this is exactly what a Householder QR decomposition algorithm does, premultiplying a given matrix by a sequence of unitary/orthogonal matrices to triangularize it. – J. M. Oct 3 '11 at 17:52
@J. M.: Yes -- I know the answer was basically there in the comments, but since the bounty was created after your first comment on QR decomposition, I figured the OP wanted a more detailed explanation. Also I followed your link in the last comment and couldn't view the crucial page there. – joriki Oct 3 '11 at 18:06

After some reading (I actually wanted to solve the exact same problem that Edison posted here), I came across a paper called "A basis-kernel representation of orthogonal matrices" (1995) by X Sun, C H Bischof.

The idea behind the paper is that orthogonal matrices can be expressed as Q = I - YSY'. If S is upper triangular, then the Y matrix contains the Householder vectors, with the diagonal of S containing scaling that needs to be applied on top of those.

After some experimentation, I was able to come up with the following MATLAB code that breaks an orthogonal matrix down to the same amount of Householder vectors as the dimension of the provided orthogonal matrix - assuming it has full rank...

function HHj=HH_decompose(AA)   % if AA not unitary, you'll get weird results...
ksize=length(AA);
A=eye(ksize)-AA;
B=eye(ksize);
for k=1:(ksize-1)
z=-A((k+1):ksize,k)'*inv(A((k+1):ksize,(k+1):ksize))';
B(k,(k+1):ksize)=z;
end
Binv=inv(B);
A=B*A*B';
C=diag(sqrt(2./diag(A)));
Binv=Binv*inv(C);
HHj=Binv;
end


PS: I have yet to make the code more robust (it works on most random unitary matrices I throw at it).

Here's a matrix I tested it with:

test =

0.2945   -0.1150    0.9052   -0.0571    0.2782
-0.0727   -0.6897    0.0288   -0.5835   -0.4215
-0.6306   -0.5117    0.0172    0.3441    0.4709
-0.4810    0.2050    0.4099    0.3243   -0.6734
0.5281   -0.4552   -0.1070    0.6578   -0.2642


And the Householder vectors that will generate it:

>> HH_decompose(test)

ans =

1.0000   -0.3749    0.7095    0.5321   -0.3321
0    0.9271    0.3454   -0.3124    0.2863
0         0    0.6142   -0.3259    0.0673
0         0         0    0.7163   -0.4137
0         0         0         0    0.7950

-
I voted you up here because I can appreciate when specific code is posted as it is always more precise than language. Along that spirit though, have you determined what kind of inputs can give problems? The original question is a simple one once the details of Householder reflections are taken as a given, and no unitary input is problematic. – adam W Apr 7 '13 at 16:06
I generated some random unitary matrices using ss in [ss vv dd]=svd(randn(5));. After reconstruction, the worst offset was after 7 decimal places. I considered pasting a bad-acting matrix here, but rounding makes things MUCH worse - the reconstructed matrix doesn't match at all. – Ru Hasha Apr 8 '13 at 17:55