1 | function [h] = householder_matrix(a, v) |
---|
2 | % h = householder_matrix(a, v) returns the Householder matrix that |
---|
3 | % will zero all elements of a except those corresponding to (any) |
---|
4 | % non-zero elements of v. If v has zero norm, an identity matrix is |
---|
5 | % returned. The matrix returned is either a left or right Householder |
---|
6 | % matrix dependent on whether the input parameters are column or row |
---|
7 | % vectors respectively. |
---|
8 | |
---|
9 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
10 | % See the file : Copyright.m for further details. |
---|
11 | |
---|
12 | % Reference: |
---|
13 | % |
---|
14 | % Sangwine, S. J. and Le Bihan, N., |
---|
15 | % Quaternion singular value decomposition based on bidiagonalization |
---|
16 | % to a real or complex matrix using quaternion Householder transformations, |
---|
17 | % Applied Mathematics and Computation, 182(1), 1 November 2006, 727-738, |
---|
18 | % DOI:10.1016/j.amc.2006.04.032. |
---|
19 | % |
---|
20 | % S. J. Sangwine and N. Le Bihan, |
---|
21 | % Quaternion Singular Value Decomposition based on Bidiagonalization |
---|
22 | % to a Real Matrix using Quaternion Householder Transformations, |
---|
23 | % arXiv:math.NA/0603251, 10 March 2006. Available at http://www.arxiv.org/. |
---|
24 | |
---|
25 | error(nargchk(2, 2, nargin)), error(nargoutchk(1, 1, nargout)) |
---|
26 | |
---|
27 | if size(a) ~= size(v) |
---|
28 | error('Input parameters must be vectors of the same size, both row, or both column') |
---|
29 | end |
---|
30 | |
---|
31 | [m, n] = size(a); |
---|
32 | |
---|
33 | col_vector = n == 1 & m == length(a); |
---|
34 | row_vector = m == 1 & n == length(a); |
---|
35 | |
---|
36 | if col_vector == 0 & row_vector == 0 |
---|
37 | error('Parameters must be either column or row vectors.'); |
---|
38 | end |
---|
39 | |
---|
40 | h = quaternion(eye(length(a))); |
---|
41 | |
---|
42 | n = norm(v); % If v has zero norm, we return an identity matrix. |
---|
43 | |
---|
44 | if n ~= 0 |
---|
45 | [u, zeta] = householder_vector(a, v ./ n); |
---|
46 | if size(u) ~= size(a) |
---|
47 | error('Result from function householder_vector is not of correct type (row, column)'); |
---|
48 | end |
---|
49 | if col_vector |
---|
50 | h = (1 ./ zeta) .* (h - u * u'); |
---|
51 | else |
---|
52 | h = (h - u.' * conj(u)) .* (1 ./ zeta); |
---|
53 | end |
---|
54 | end |
---|