1 | function R = inv(a) |
---|
2 | % INV Inverse of a quaternion (matrix). |
---|
3 | % (Quaternion overloading of standard Matlab function.) |
---|
4 | % |
---|
5 | % For a single quaternion, this function computes the inverse, |
---|
6 | % that is the conjugate divided by the modulus. The result will |
---|
7 | % be a NaN if the quaternion has zero modulus. For matrices it |
---|
8 | % computes the matrix inverse using an analytical formula based |
---|
9 | % on partitioning the matrix into sub-matrices. A better method |
---|
10 | % may be substituted in the future. |
---|
11 | |
---|
12 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
13 | % See the file : Copyright.m for further details. |
---|
14 | |
---|
15 | % Reference: Wikipedia article on 'Matrix Inversion', 12 July 2005. |
---|
16 | % Wikipedia article on 'Schur Complement', 13 July 2005. |
---|
17 | % |
---|
18 | % A similar formula is given in: |
---|
19 | % |
---|
20 | % R. A. Horn and C. R. Johnson, 'Matrix Analysis', |
---|
21 | % Cambridge University Press, 1985, §0.7.3, page 18. |
---|
22 | |
---|
23 | [r, c] = size(a); |
---|
24 | |
---|
25 | if r ~= c |
---|
26 | error('Matrix must be square.'); |
---|
27 | end |
---|
28 | |
---|
29 | if r == 1 % c must also be 1, since we have checked that r == c above. |
---|
30 | R = conj(a) ./ modsquared(a); |
---|
31 | elseif r == 2 |
---|
32 | |
---|
33 | % We handle the 2 by 2 case separately since this can be computed |
---|
34 | % using quaternion products, rather than quaternion matrix products. |
---|
35 | % (See note at end on use of substruct function.) |
---|
36 | |
---|
37 | A = subsref(a, substruct('()', {1,1})); % A = a(1,1); |
---|
38 | B = subsref(a, substruct('()', {1,2})); % B = a(1,2); |
---|
39 | C = subsref(a, substruct('()', {2,1})); % C = a(2,1); |
---|
40 | D = subsref(a, substruct('()', {2,2})); % D = a(2,2); |
---|
41 | |
---|
42 | IA = inv(A); % If A is zero, this will cause a NaN and the whole |
---|
43 | % result will be NaNs. Solution unknown as yet. |
---|
44 | |
---|
45 | IAB = IA .* B; % These two subexpressions are needed more than |
---|
46 | CIA = C .* IA; % once each, and are computed once for efficiency. |
---|
47 | |
---|
48 | T = inv(D - CIA .* B); |
---|
49 | |
---|
50 | TCIA = T .* CIA; % This is needed twice, so compute it once and reuse. |
---|
51 | |
---|
52 | R = [[IA + IAB .* TCIA, -IAB .* T];... |
---|
53 | [ -TCIA, T]]; |
---|
54 | else |
---|
55 | |
---|
56 | % For matrices larger than 2 by 2, we partition the matrix into sub |
---|
57 | % blocks of half the height/width of the input matrix. If the input |
---|
58 | % matrix has an even number of rows (and therefore columns, since it |
---|
59 | % is square), then the sub-blocks will be exactly half the height |
---|
60 | % (width), otherwise the left block will have one less element, |
---|
61 | % because we round towards zero. The layout of the sub-blocks is: |
---|
62 | % |
---|
63 | % [ A B ] |
---|
64 | % [ C D ] |
---|
65 | % |
---|
66 | % the same as in the 2 by 2 case above (except that here they are |
---|
67 | % matrices and not scalars). |
---|
68 | |
---|
69 | p = fix(r./2); % Calculate half the number of rows/columns rounded down. |
---|
70 | q = p + 1; % Index of first element of right or bottom block. |
---|
71 | |
---|
72 | A = subsref(a, substruct('()', {1:p, 1:p})); % A = a(1:p, 1:p); |
---|
73 | B = subsref(a, substruct('()', {1:p, q:c})); % B = a(1:p, q:c); |
---|
74 | C = subsref(a, substruct('()', {q:r, 1:p})); % C = a(q:r, 1:p); |
---|
75 | D = subsref(a, substruct('()', {q:r, q:c})); % D = a(q:r, q:c); |
---|
76 | |
---|
77 | IA = inv(A); |
---|
78 | |
---|
79 | IAB = IA * B; % These two subexpressions are needed more than |
---|
80 | CIA = C * IA; % once each, and are computed once for efficiency. |
---|
81 | |
---|
82 | T = inv(D - CIA * B); |
---|
83 | |
---|
84 | TCIA = T * CIA; % This is needed twice, so compute it once and reuse. |
---|
85 | |
---|
86 | R = [[IA + IAB * TCIA, -IAB * T];... |
---|
87 | [ -TCIA, T]]; |
---|
88 | end |
---|
89 | |
---|
90 | % Note on subscripted references: |
---|
91 | % |
---|
92 | % In the code above, we need to use subscripted references, but these do |
---|
93 | % not work inside a class method (which this is). A brief mention of this |
---|
94 | % can be found in the Matlab documentation under Programming/Classes and |
---|
95 | % Objects/Designing User Classes under the title Object Indexing within |
---|
96 | % Methods. The equally obscure solution is to use the substruct function, |
---|
97 | % as above. For each statement where it is used, the comment gives the |
---|
98 | % normal Matlab notation (which won't work here). |
---|