1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
function [u, p] = polardecomp(a, side)
% Compute the polar decomposition.
%
% Returns the factors of the polar decomposition [1] `u` and `p` such
% that ``a = up`` (if `side` is "right") or ``a = pu`` (if `side` is
% "left"), where `p` is positive semidefinite. Depending on the shape
% of `a`, either the rows or columns of `u` are orthonormal. When `a`
% is a square array, `u` is a square unitary array. When `a` is not
% square, the "canonical polar decomposition" [2] is computed.
%
% Parameters:
% a: The array to be factored.
% side: (optional) Determines whether a right or left polar decomposition
% is computed. If `side` is "right", then ``a = up``. If `side` is
% "left", then ``a = pu``. The default is "right".
%
% Returns:
% u: If `a` is square, then `u` is unitary. If m > n, then the columns
% of `a` are orthonormal, and if m < n, then the rows of `u` are
% orthonormal.
% p: `p` is Hermitian positive semidefinite. If `a` is nonsingular, `p`
% is positive definite. The shape of `p` is (n, n) or (m, m), depending
% on whether `side` is "right" or "left", respectively.
if nargin < 2
side = 'right';
end
if ~strcmp(side, 'right') && ~strcmp(side, 'left')
error('`side` must be either ''right'' or ''left''');
end
if ~ismatrix(a)
error('`a` must be a 2-D array.');
end
[w, s, vh] = svd(a, 'econ');
u = w * vh';
if strcmp(side, 'right')
% a = up
p = (vh' * diag(diag(s))) * vh;
else
% a = pu
p = (w * diag(diag(s))) * w';
end
end
%!test
%! a = [1,-1; 2, 4];
%! [u, p] = polardecomp(a);
%! u0 = [0.85749293, -0.51449576 ; 0.51449576, 0.85749293];
%! p0 = [1.88648444, 1.2004901; 1.2004901, 3.94446746];
%! assert(u,u0, 1e-8)
%! assert(p,p0, 1e-8)
%! assert(u*u', eye(2), 1e-15)
|
