# Mathematica Implementation of Householder’s Method

I typed the Householder code in this paper, which starts on page 7.

The code is:

 A = {{4., 1., -2., 2.}, {1., 2., 0, 1.}, {-2., 0, 3., -2.}, {2.,
1., -2., -1.}}; (* Made A numeric for speed-up purposes *)
n = Length[A[[1]]];
zeroVector = {};
For[i = 1, i <= n,
zeroVector = Append[zeroVector, {0}];
i++
];
Alist = {A};
Hlist = {};

For[j = 1, j <= n - 2,
If[A[[j + 1, j]] >= 0, c = 1, c = 2];
alpha = (-1)^c (Sum[A[[k, j]]^2, {k = j + 1, n}])^(1/2);
r = ((1/2) alpha^2 - (1/2) alpha A[[j + 1, j]])^(1/2);
x = zeroVector;
x[[j + 1, 1]] = (A[[j + 1, j]] - alpha)/(2 r);
For[k = j + 2, k <= n,
x[[k, 1]] = A[[k, j]]/(2 r);
k++];
H = IdentityMatrix[n] - 2 x.Transpose[x];
A = H.A.H;
Hlist = Append[Hlist, H]; Print[Hlist];
Alist = Append[Alist, A]; Print[Hlist];
j++];


Of course, the program is complaining with such things as "Sum::write: Tag Set in k=j+1 is Protected." I am just not seeing what is wrong. I added the prints as I want to see intermediate results. Can anyone spot the issue?

Is there a better solution or code for this?

Aside: I wish I could find a nice implementation of this somewhere, but unfortunately I have not been able to find in MMA or the web. I found the following snippet on Wolfram, but I cannot make use of it (maybe I am missing something).

 HouseholderMatrix[v_?VectorQ] :=
IdentityMatrix[Length[v]]
- 2 Transpose[{v}] . {v} / (v.v)


End Aside

• Sorry, but this is bad code :) zeroVector = {}; For[i = 1, i <= n, zeroVector = Append[zeroVector, {0}]; i++]; can be replaced by one Table[] command. I see you copied it from the paper/article. But papers are full of bad code as well. Looks like written by an ex-Fortran programmer. – Nasser Apr 14 '14 at 2:22
• I agree with your comment of bad code. That is why I have the aside in my post. However, I just want something quick to see HH transformations from step to step. Do you know of an alternate approach or code or something built-in to MMA? – Amzoti Apr 14 '14 at 2:24
• For the Sum::write: Tag Set in k=j+1 is error, simply change k= to k,. Now Alist is generated. I do not know if it is correct or not. – Nasser Apr 14 '14 at 4:02
• That did it! Knew it was something silly. Now I have to study the results! Thank you! Regards – Amzoti Apr 14 '14 at 4:06
• BTW: ReflectionMatrix[] is built-in; that can be repurposed to generate Householder matrices. – J. M.'s torpor May 16 '16 at 13:48

The implementation you linked to works fine. I do not understand why you could not make it work.

Verified it with Maple build-in function. Same result:

 LinearAlgebra[HouseholderMatrix](<1,2,3,4,5>);


Mathematica:

HouseholderMatrix[v_?VectorQ] := IdentityMatrix[Length[v]] -
2 Transpose[{v}].{v}/(v.v);
HouseholderMatrix[{1, 2, 3, 4, 5}];
MatrixForm[%]


• Maybe I am missing something, but how do you use it to get the same result as the $4x4$ in the paper? – Amzoti Apr 14 '14 at 3:25
• @Amzoti I have not looked at the paper result. I verified it using Maple. Same matrix resulted from same input. Will try to look at the paper. May be the paper is using different definition or has a bug :) – Nasser Apr 14 '14 at 3:30
• The example here uses a single vector, but we want to use it for a matrix. I suspect that you somehow call it for each row vector and somehow combine those. I am also concerned that this snippet does not use the sign test to determine the signs. Note, this paper is just an example, this is just Householder transforms. The paper was used only because it had code. Here is a Wiki example: en.wikipedia.org/wiki/Householder_transformation – Amzoti Apr 14 '14 at 3:32
• I see now. The Householder matrix is one step used in the Householder method. The code you linked to only generates the Householder matrix, it does not implement the full transformation. Here is a Mathematica implementaion I found: mathfaculty.fullerton.edu/mathews/n2003/HouseholderMod.html – Nasser Apr 14 '14 at 3:41
• I am very familiar with this site, unfortunately, all of the code has been removed and the stuff is not copyable. I would have to retype that code and hope no errors are made. Thanks anyway! Regards – Amzoti Apr 14 '14 at 3:46

The correct Mathematica code (v.10.2.0.0.) is:

A = N[{{4, 1, -2, 2}, {1, 2, 0, 1}, {-2, 0, 3, -2}, {2, 1, -2, -1}}];
(*A=N[{{-42,43,-2,28},{43,-98,72,-26},{-2,72,-96,53},{28,-26,53,54}}];*)
n = Length[A[[1]]];
zeroVector = {};
For[i = 1, i <= n, i++, zeroVector = Append[zeroVector, {0}]];
Alist = {A};
Hlist = {};

For[j = 1, j <= n - 2, j++,

If[A[[j + 1, j]] >= 0, c = 1, c = 2];
alpha = (-1)^c (Sum[A[[k, j]]^2, {k, j + 1, n}])^(1/2);
r = ((1/2) alpha^2 - (1/2) alpha A[[j + 1, j]])^(1/2);
x = zeroVector;
x[[j + 1, 1]] = (A[[j + 1, j]] - alpha)/(2 r);

For[k = j + 2, k <= n, k++,
x[[k, 1]] = A[[k, j]]/(2 r)];

H = IdentityMatrix[n] - 2 x.Transpose[x];
A = H.A.H;
Hlist = Append[Hlist, H];
Alist = Append[Alist, A];
]

MatrixForm[Chop[A]]


considering example MIHM from the paper, we find a the tridiagonal matrix A' of Matrix A:

     4  1 -2  2                 4      3          0          0
A =  1  2  0  1     ->    A' =  3      3.33333   -1.66667    0
-2  0  3 -2                 0     -1.66667   -1.32      -0.906667
2  1 -2 -1                 0      0         -0.906667   1.98667


This result is different from the result in the paper. However, we find the same eigenvalues/vectors for A and A', which is a strong indication that the code is correct.

By using the same code to calculate the tridiagonal matrix of example CNAHM on page 9 leads also to the correct answer.

It seems that the result for example MIHM in the paper is not the correct tridiagonal form of matrix A.