1 | function [u, zeta] = householder_vector(a, v) |
---|
2 | % [u, zeta] = householder_vector(a, v) calculates a Householder vector, |
---|
3 | % u, with a norm of sqrt(2), and a value zeta, from the vectors a and v. The |
---|
4 | % results may be used to construct a left or right Householder matrix and they |
---|
5 | % depend on whether the input parameters are column or row vectors respectively. |
---|
6 | % The result returned in u has the same type (row, column) as the parameter a. |
---|
7 | % zeta will be real, complex or quaternion, depending on the type of a. |
---|
8 | % The parameter v must be real (mathematical limitation). |
---|
9 | |
---|
10 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
11 | % See the file : Copyright.m for further details. |
---|
12 | |
---|
13 | error(nargchk(2, 2, nargin)), error(nargoutchk(2, 2, nargout)) |
---|
14 | |
---|
15 | if size(a) ~= size(v) |
---|
16 | error('Input parameters must be vectors of the same size, both row, or both column') |
---|
17 | end |
---|
18 | |
---|
19 | if isa(v, 'quaternion') |
---|
20 | error('Parameter v may not be a quaternion vector (mathematical limitation)') |
---|
21 | end |
---|
22 | |
---|
23 | [m, n] = size(a); |
---|
24 | |
---|
25 | col_vector = n == 1 & m == length(a); |
---|
26 | row_vector = m == 1 & n == length(a); |
---|
27 | |
---|
28 | if col_vector == 0 & row_vector == 0 |
---|
29 | error('Parameters must be either column or row vectors.'); |
---|
30 | end |
---|
31 | |
---|
32 | % References: |
---|
33 | % |
---|
34 | % Alston S. Householder, Unitary Triangularization of a Nonsymmetric Matrix, |
---|
35 | % J. ACM, 5 (4), 1958, pp339--342. |
---|
36 | % DOI:10.1145/320941.320947 |
---|
37 | % |
---|
38 | % David D. Morrison, |
---|
39 | % Remarks on the Unitary Triangularization of a Nonsymmetric Matrix, |
---|
40 | % J. ACM, 7 (2), 1960, pp185--186. |
---|
41 | % DOI:10.1145/321021.321030 |
---|
42 | % |
---|
43 | % Sangwine, S. J. and Le Bihan, N., |
---|
44 | % Quaternion singular value decomposition based on bidiagonalization |
---|
45 | % to a real or complex matrix using quaternion Householder transformations, |
---|
46 | % Applied Mathematics and Computation, 182(1), 1 November 2006, 727-738, |
---|
47 | % DOI:10.1016/j.amc.2006.04.032. |
---|
48 | % |
---|
49 | % S. J. Sangwine and N. Le Bihan, |
---|
50 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
51 | % to a Real Matrix using Quaternion Householder Transformations, |
---|
52 | % arXiv:math.NA/0603251, 10 March 2006. Available at http://www.arxiv.org/. |
---|
53 | |
---|
54 | alpha = norm(a); |
---|
55 | |
---|
56 | if alpha == 0 |
---|
57 | u = a - a; % This ensures that u is zero, even if a is a |
---|
58 | % complex quaternion vector with zero norm. |
---|
59 | zeta = 1; |
---|
60 | return; |
---|
61 | end |
---|
62 | |
---|
63 | if col_vector |
---|
64 | romega = a.' * v; |
---|
65 | else |
---|
66 | romega = v * a.'; |
---|
67 | end |
---|
68 | |
---|
69 | r = abs(romega); |
---|
70 | |
---|
71 | if r ~= 0 |
---|
72 | zeta = -romega ./ r; |
---|
73 | else |
---|
74 | zeta = 1; |
---|
75 | end |
---|
76 | |
---|
77 | mu = sqrt(alpha .* (alpha + r)); |
---|
78 | |
---|
79 | if col_vector |
---|
80 | u = (1 ./ mu) .* (a - zeta .* v .* alpha); |
---|
81 | else |
---|
82 | u = (1 ./ mu) .* conj(alpha .* v .* zeta - a); |
---|
83 | end |
---|