draft of bloch_messiah

1 files changed, 105 insertions, 0 deletions

diff --git a/python_originals/polar.py b/python_originals/polar.py
new file mode 100644
index 0000000..69b6a35
--- /dev/null
+++ b/python_originals/polar.py
@@ -0,0 +1,105 @@
+def polar(a, side="right"):
+    """
+    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
+    --------
+    >>> import numpy as np
+    >>> 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
+