import numpy as np
import scipy.linalg


$$\newcommand{\trace}{\operatorname{tr}}$$

For a given matrix $$\hat{A}$$, we want to find the closest unitary matrix $$\hat{U}$$, in the sense that the operator norm (aka 2-norm) of their difference should be minimal.

$$\hat{U} = \arg \min_{\tilde{U}} \Vert \hat{A}-\tilde{U} \Vert$$

For example, let’s consider the following matrix:

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]
])


In cases where $$\hat{A}$$ of dimension $$n$$ is obtained from projecting down a unitary matrix $$\hat{U}_m$$ of a larger dimension ($$\hat{A} = \hat{P} \hat{U}_m \hat{P}$$, where $$\hat{P}$$ is the projector from dimension $$m$$ to dimension $$n$$) , one possible way to quantify the “distance from unitarity” is

$$d_u(\hat{A}) = 1 - \frac{1}{n} \trace[\hat{A}^\dagger \hat{A}].$$

This situation is common in quantum information, where $$\hat{A}$$ is the projection of the unitary evolution of a large Hilbert space of dimension $$m$$ into a small logical subspace of dimension $$n$$. The quantity $$d_u(\hat{A})$$ is then simply the population lost from logical subspace.

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]


For the given matrix, we lose about 80% of the population in the logical subspace:

>>> delta_uni(A)
0.80861752310942581


A $$\hat{U}$$ that minimizes $$\Vert \hat{A} - \hat{U}\Vert$$ can be calculated via a singular value decomposition (SVD) [Reich]:

$$\hat{A} = \hat{V} \hat{\Sigma} \hat{W}^\dagger,$$
$$\hat{U} = \hat{V} \hat{W}^\dagger.$$
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)


The SVD also allows to calculate the distance that $$\hat{A}$$ has from $$\hat{U}$$.

$$d(\hat{A}, \hat{U}) = \max_i \vert \sigma_i - 1\vert,$$

where $$\sigma_i$$ are the diagonal entries of $$\hat{\Sigma}$$ from the SVD. This is a more general measure of “distance from unitarity” than delta_uni

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


For matrices obtained from projecting down from a larger Hilbert space, the maximum distance is 1. For general matrices, it can be larger.

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

>>> deltaU(U) # should be zero, within machine precision
6.6613381477509392e-16


We can double check this with the implemenation of the 2-norm in SciPy:

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


## Reference

[Reich] D.M.Reich. “Characterisation and Identification of Unitary Dynamics Maps in Terms of Their Action on Density Matrices” (unpublished)