1 files changed, 56 insertions, 0 deletions
diff --git a/polardecomp.m b/polardecomp.m
new file mode 100644
index 0000000..698bd90
--- /dev/null
+++ b/polardecomp.m
@@ -0,0 +1,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)