Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}];
 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);
      H = IdentityMatrix[n] - 2 x.Transpose[x];
      A = H.A.H;
     Hlist = Append[Hlist, H]; Print[Hlist];
     Alist = Append[Alist, A]; Print[Hlist];

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] :=
  - 2 Transpose[{v}] . {v} / (v.v)

End Aside

share|improve this question
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
up vote 1 down vote accepted

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:


Mathematica graphics


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

Mathematica graphics

share|improve this answer
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: – 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: – 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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.