[37] | 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 |
---|