Mercurial > matrix-functions

comparison toolbox/house.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function [v, beta] = house(x)
2 %HOUSE Householder matrix.
3 % If [v, beta] = HOUSE(x) then H = EYE - beta*v*v' is a Householder
4 % matrix such that Hx = -sign(x(1))*norm(x)*e_1.
5 % NB: If x = 0 then v = 0, beta = 1 is returned.
6 % x can be real or complex.
7 % sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
8
9 % Theory: (textbook references Golub & Van Loan 1989, 38-43;
10 % Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
11 % Hx = y: (I - beta*v*v')x = -s*e_1.
12 % Must have |s| = norm(x), v = x+s*e_1, and
13 % x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
14 % So take s = sign(x(1))*norm(x) (which avoids cancellation).
15 % v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
16 % = 2*norm(x)*(norm(x) + |x(1)|).
17 %
18 % References:
19 % G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
20 % Johns Hopkins University Press, Baltimore, Maryland, 1989.
21 % G.W. Stewart, Introduction to Matrix Computations, Academic Press,
22 % New York, 1973,
23 % J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
24 % Press, 1965.
25
26 [m, n] = size(x);
27 if n > 1, error('Argument must be a column vector.'), end
28
29 s = norm(x) * (sign(x(1)) + (x(1)==0)); % Modification for sign(0)=1.
30 v = x;
31 if s == 0, beta = 1; return, end % Quit if x is the zero vector.
32 v(1) = v(1) + s;
33 beta = 1/(s'*v(1)); % NB the conjugated s.
34
35 % beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
36 % But beta as above can be non-real (due to rounding) only when x is complex.