复 2x2 矩阵的奇异值分解

Question

我正在寻找示例代码，展示如何计算可以包含复数值的 2x2 矩阵的 singular value decomposition。

例如，这对于将用户输入的矩阵“修复”为单一矩阵很有用。您只需使用u, s, v = svd(m)，然后从产品中省略s 部分：repaired = u * v。

Answer 1

这里有一些 python 代码可以解决问题。它基本上只是提取复杂的部分，然后委托给the solution from this answer for real 2x2 matrices。

我已经使用 numpy 在 python 中编写了代码。这有点讽刺，因为如果你有 numpy 你应该只使用np.linalg.svd。显然，这旨在作为示例代码，适合在紧要关头学习或翻译成其他语言。

我也不是数值稳定性方面的专家，所以……买家要小心。

import numpy as np
import math

# Note: in practice in python just use np.linalg.svd instead

def singular_value_decomposition_complex_2x2(m):
    """
    Returns a singular value decomposition of the given 2x2 complex numpy
    matrix.

    :param m: A 2x2 numpy matrix with complex values.
    :returns: A tuple (U, S, V) where U*S*V ~= m, where U and V are complex
    2x2 unitary matrices, and where S is a 2x2 diagonal matrix with
    non-negative real values.
    """

    # Make top row non-imaginary and non-negative by column phasing.
    # m2 = m p = | >     >    |
    #            | ?+?i  ?+?i |
    p = phase_cancel_matrix(m[0, 0], m[0, 1])
    m2 = m * p

    # Cancel top-right value by rotation.
    # m3 = m p r = | ?+?i  0    |
    #              | ?+?i  ?+?i |
    r = rotation_matrix(math.atan2(m2[0, 1].real, m2[0, 0].real))
    m3 = m2 * r

    # Make bottom row non-imaginary and non-negative by column phasing.
    # m4 = m p r q = | ?+?i  0 |
    #                | >     > |
    q = phase_cancel_matrix(m3[1, 0], m3[1, 1])
    m4 = m3 * q

    # Cancel imaginary part of top left value by row phasing.
    # m5 = t m p r q = | > 0 |
    #                  | > > |
    t = phase_cancel_matrix(m4[0, 0], 1)
    m5 = t * m4

    # All values are now real (also the top-right is zero), so delegate to a
    # singular value decomposition that works for real matrices.
    # t m p r q = u s v
    u, s, v = singular_value_decomposition_real_2x2(np.real(m5))

    # m = (t* u) s (v q* r* p*)
    return adjoint(t) * u, s, v * adjoint(q) * adjoint(r) * adjoint(p)


def singular_value_decomposition_real_2x2(m):
    """
    Returns a singular value decomposition of the given 2x2 real numpy matrix.

    :param m: A 2x2 numpy matrix with real values.
    :returns: A tuple (U, S, V) where U*S*V ~= m, where U and V are 2x2
    rotation matrices, and where S is a 2x2 diagonal matrix with
    non-negative real values.
    """
    a = m[0, 0]
    b = m[0, 1]
    c = m[1, 0]
    d = m[1, 1]

    t = a + d
    x = b + c
    y = b - c
    z = a - d

    theta_0 = math.atan2(x, t) / 2.0
    theta_d = math.atan2(y, z) / 2.0

    s_0 = math.sqrt(t**2 + x**2) / 2.0
    s_d = math.sqrt(z**2 + y**2) / 2.0

    return \
        rotation_matrix(theta_0 - theta_d), \
        np.mat([[s_0 + s_d, 0], [0, s_0 - s_d]]), \
        rotation_matrix(theta_0 + theta_d)


def adjoint(m):
    """
    Returns the adjoint, i.e. the conjugate transpose, of the given matrix.
    When the matrix is unitary, the adjoint is also its inverse.

    :param m: A numpy matrix to transpose and conjugate.
    :return: A numpy matrix.
    """
    return m.conjugate().transpose()


def rotation_matrix(theta):
    """
    Returns a 2x2 unitary matrix corresponding to a 2d rotation by the given angle.

    :param theta: The angle, in radians, that the matrix should rotate by.
    :return: A 2x2 orthogonal matrix.
    """
    c, s = math.cos(theta), math.sin(theta)
    return np.mat([[c, -s],
                   [s, c]])


def phase_cancel_complex(c):
    """
    Returns a unit complex number p that cancels the phase of the given complex
    number c. That is, c * p will be real and non-negative (approximately).

    :param c: A complex number.
    :return: A complex number on the complex unit circle.
    """
    m = abs(c)

    # For small values, where the division is in danger of exploding small
    # errors, use trig functions instead.
    if m < 0.0001:
        theta = math.atan2(c.imag, c.real)
        return math.cos(theta) - math.sin(theta) * 1j

    return (c / float(m)).conjugate()


def phase_cancel_matrix(p, q):
    """
    Returns a 2x2 unitary matrix M such that M cancels out the phases in the
    column {{p}, {q}} so that the result of M * {{p}, {q}} should be a vector
    with non-negative real values.

    :param p: A complex number.
    :param q: A complex number.
    :return: A 2x2 diagonal unitary matrix.
    """
    return np.mat([[phase_cancel_complex(p), 0],
                   [0, phase_cancel_complex(q)]])

我测试了上面的代码，用 [-10, 10] + [-10, 10]i 中的随机值填充的矩阵进行模糊测试，并检查分解后的因子是否具有正确的属性（即酉、对角线、实数） ，视情况而定）并且他们的产品（大约）等于输入。

但这里有一个简单的冒烟测试：

m = np.mat([[5, 10], [1j, -1]])
u, s, v = singular_value_decomposition_complex_2x2(m)

np.set_printoptions(precision=5, suppress=True)

print "M:\n", m
print "U*S*V:\n", u*s*v
print "U:\n", u
print "S:\n", s
print "V:\n", v
print "M ~= U*S*V:", np.all(np.abs(m - u*s*v) < 0.1**14)

输出以下内容。您可以确认分解后的 S 与 svd from wolfram alpha 匹配，当然 U 和 V 可以（并且是）不同。

M:
[[  5.+0.j  10.+0.j]
 [  0.+1.j  -1.+0.j]]
U*S*V:
[[ 5.+0.j 10.+0.j]
 [ 0.+1.j -1.-0.j]]
U:
[[-0.89081-0.44541j  0.08031+0.04016j]
 [ 0.08979+0.j       0.99596+0.j     ]]
S:
[[ 11.22533   0.     ]
 [  0.        0.99599]]
V:
[[-0.39679+0.20639j -0.80157+0.39679j]
 [ 0.40319+0.79837j -0.19359-0.40319j]]
M ~= U*S*V: True