|
""" 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 : (m, n) array_like The array to be factored. side : {'left', 'right'}, 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 : (m, n) ndarray 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 : ndarray `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.
References ---------- .. [1] R. A. Horn and C. R. Johnson, "Matrix Analysis", Cambridge University Press, 1985. .. [2] N. J. Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008.
Examples -------- >>> from scipy.linalg import polar >>> a = np.array([[1, -1], [2, 4]]) >>> u, p = polar(a) >>> u array([[ 0.85749293, -0.51449576], [ 0.51449576, 0.85749293]]) >>> p array([[ 1.88648444, 1.2004901 ], [ 1.2004901 , 3.94446746]])
A non-square example, with m < n:
>>> b = np.array([[0.5, 1, 2], [1.5, 3, 4]]) >>> u, p = polar(b) >>> u array([[-0.21196618, -0.42393237, 0.88054056], [ 0.39378971, 0.78757942, 0.4739708 ]]) >>> p array([[ 0.48470147, 0.96940295, 1.15122648], [ 0.96940295, 1.9388059 , 2.30245295], [ 1.15122648, 2.30245295, 3.65696431]]) >>> u.dot(p) # Verify the decomposition. array([[ 0.5, 1. , 2. ], [ 1.5, 3. , 4. ]]) >>> u.dot(u.T) # The rows of u are orthonormal. array([[ 1.00000000e+00, -2.07353665e-17], [ -2.07353665e-17, 1.00000000e+00]])
Another non-square example, with m > n:
>>> c = b.T >>> u, p = polar(c) >>> u array([[-0.21196618, 0.39378971], [-0.42393237, 0.78757942], [ 0.88054056, 0.4739708 ]]) >>> p array([[ 1.23116567, 1.93241587], [ 1.93241587, 4.84930602]]) >>> u.dot(p) # Verify the decomposition. array([[ 0.5, 1.5], [ 1. , 3. ], [ 2. , 4. ]]) >>> u.T.dot(u) # The columns of u are orthonormal. array([[ 1.00000000e+00, -1.26363763e-16], [ -1.26363763e-16, 1.00000000e+00]])
""" if side not in ['right', 'left']: raise ValueError("`side` must be either 'right' or 'left'") a = np.asarray(a) if a.ndim != 2: raise ValueError("`a` must be a 2-D array.")
w, s, vh = svd(a, full_matrices=False) u = w.dot(vh) if side == 'right': # a = up p = (vh.T.conj() * s).dot(vh) else: # a = pu p = (w * s).dot(w.T.conj()) return u, p |