On Mon, Mar 29, 2010 at 12:54 PM, Jitse Niesen
<jitse@xxxxxxxxxxxxxxxxx> wrote:
Hello,
I have some questions about the function householderCoefficients() in HessenbergDecomposition, which returns the vector of coefficients for the Householder reflections which reduce the given matrix to Hessenberg form.
Why does it return a vector of length N-1 (where N is the size of the matrix whose Hessenberg decomposition is being computed)? I think that you need only N-2 Householder reflections to reduce a matrix to Hessenberg form. The last two columns of the matrix are already in the correct form, so only N-2 columns need to be transformed. See also for instance Algorithm 7.4.2 in Golub & Van Loan. This should perhaps also be taken into account when computing the Hessenberg decomposition.
For real it's true, but not for complexes where we ensure the sub diagonal entries are real. So I think it is better to always have a vector of length N-1. Well at least it seems simpler. Of course if you think avoiding the last useless reflector for real has a non negligible impact on the performance, then why not...
While I was thinking about this, I started to wonder whether we should store the coefficients at all. For real matrices, the coefficients are given by 2 / |v|^2 where v is the Householder vector so it is not necessary to store them. On the other hand, we do save some effort if we do not recompute the coefficients - though probably not very much - and the memory requirements for storing the coefficients is only O(N). In the complex case, the coefficients are not uniquely determined by the vectors IIRC, and it makes sense to use the same algorithm for both real and complex scalars if possible. So there seems to be a good argument for sticking to the current strategy of storing the coefficients.
Indeed. Moreover computing |v|^2 is rather expensive, so it's better to store them.
The EigenSolver class uses its own routine for reducing the matrix to Hessenberg form. I assume that this is simply because nobody has had time to rewrite it to use HessenbergDecomposition.
Exactly ;) The EigenSolver has to be entirely rewritten from scratch, though the different heuristics which are involved might be good to keep.
gael
Cheers,
Jitse