import numpy as np
import scipy.linalg

A = np.matrix([
[-1.75900766E-02-1.15354406E-01j, 6.10816904E-03+9.49971160E-03j, 1.79090787E-02+1.33311069E-02j,-1.82163102E-03-8.77682357E-04j],
[-9.77987203E-03+5.01950535E-03j, 8.74085180E-04+3.25580543E-04j,-6.74874670E-03-5.82800747E-03j,-1.95106265E-03+9.84138284E-04j],
[-6.11175534E-03+2.26761191E-02j,-5.04355339E-03-2.57661178E-02j, 2.15674643E-01+8.36337993E-01j, 1.76098908E-02+1.74391779E-02j],
[ 1.51473418E-03+1.07178057E-03j, 6.40793740E-04-1.94176372E-03j,-1.28408900E-02+2.66263921E-02j, 4.84726807E-02-3.84341687E-03j]
])

def delta_uni(A):
    """ Assuming that A is a matrix obtained from projecting a unitary matrix
        to a smaller subspace, return the loss of population of subspace, as a
        distance measure of A from unitarity.
        
        Result is in [0,1], with 0 if A is unitary.
    """
    return 1.0 - (A.H * A).trace()[0,0].real / A.shape[0]

delta_uni(A)

def closest_unitary(A):
    """ Calculate the unitary matrix U that is closest with respect to the
        operator norm distance to the general matrix A.

        Return U as a numpy matrix.
    """
    V, __, Wh = scipy.linalg.svd(A)
    U = np.matrix(V.dot(Wh))
    return U


U = closest_unitary(A)

def deltaU(A):
    """ Calculate the operator norm distance \Vert\hat{A} - \hat{U}\Vert between
        an arbitrary matrix $\hat{A}$ and the closest unitary matrix $\hat{U}$
    """
    __, S, __ = scipy.linalg.svd(A)
    d = 0.0
    for s in S:
        if abs(s - 1.0) > d:
            d = abs(s-1.0)
    return d

deltaU(A) # should be close to 1, we are *very* non-unitary

deltaU(U) # should be zero, within machine precision

scipy.linalg.norm(A-U, ord=2) # should be equal to deltaU(A)