複數2×2矩陣的奇異值分解

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