NumPy SVD与R实现不一致

kqlmhetl 于 2023-04-21 发布在其他

关注(0)|答案(1)|浏览(143)

我在Stack Overflow上看到一个question about inverting a singular matrix使用NumPy。我想看看NumPy SVD是否可以提供一个可接受的答案。
我已经演示了在R中使用SVD来获得另一个Stack Overflow答案。我使用已知的解决方案来确保我的NumPy代码在应用于新问题之前能够正确工作。
我很惊讶地发现NumPy的解与R的答案不匹配。当我把NumPy的解代入方程时，我没有得到一个恒等式。
来自R和NumPy的U矩阵具有相同的形状（3x 3），并且值相同，但是符号不同。下面是我从NumPy得到的U矩阵：

对于R和NumPy，D矩阵是相同的。这里是大对角元素被归零之后的D：

我从NumPy得到的V矩阵具有3x 4的形状;R给了我一个4x 3的矩阵。值是相似的，但是符号不同，就像U一样。下面是我从NumPy得到的V矩阵：

R解向量为：

x = [2.41176,-2.28235,2.15294,-3.47059]

当我把它代入原始方程A*x = b时，我从R解中得到RHS向量：

b = [-17.00000,28.00000,11.00000]

NumPy给了我这个解决方案向量：

x = [2.55645,-2.27029,1.98412,-3.23182]

当我将NumPy解代入原始方程A*x = b时，我得到了这个结果：

b = [-15.93399,28.04088,12.10690]

接近，但不正确。
我用NumPy np.linalg.pinv伪逆方法重复了实验，结果与R解一致。
以下是我的完整Python脚本：

# https://stackoverflow.com/questions/75998775/python-vs-matlab-why-my-matrix-is-singular-in-python

import numpy as np

def pseudo_inverse_solver(A, b):
    A_inv = np.linalg.pinv(A)
    x = np.matmul(A_inv, b)
    error = np.matmul(A, x) - b
    return x, error, A_inv

def svd_solver(A, b):
    U, D, V = np.linalg.svd(A, full_matrices=False)
    D_diag = np.diag(np.diag(np.reciprocal(D)))
    D_zero = np.array(D_diag)
    D_zero[D_zero >= 1.0e15] = 0.0
    D_zero = np.diag(D_zero)
    A_inv = np.matmul(np.matmul(np.transpose(V), D_zero), U)
    x = np.matmul(A_inv, b)
    error = np.matmul(A, x) - b
    return x, error, A_inv

if __name__ == '__main__':
    """
    Solution from my SO answer
    https://stackoverflow.com/questions/19763698/solving-non-square-linear-system-with-r/19767525#19767525
    Example showing how to use NumPy SVD
    https://stackoverflow.com/questions/24913232/using-numpy-np-linalg-svd-for-singular-value-decomposition
    """
    np.set_printoptions(20)
    A = np.array([
        [0.0, 1.0, -2.0, 3.0],
        [5.0, -3.0, 1.0, -2.0],
        [5.0, -2.0, -1.0, 1.0]
    ])
    b = np.array([-17.0, 28.0, 11.0]).T

    x_svd, error_svd, A_inv_svd = svd_solver(A, b)
    error_svd_L2 = np.linalg.norm(error_svd)
    x_pseudo, error_pseudo, A_inv_pseudo = pseudo_inverse_solver(A, b)
    error_pseudo_L2 = np.linalg.norm(error_pseudo)

有什么建议可以告诉我NumPy SVD有什么遗漏吗？我在这一行做错了吗？

A_inv = np.matmul(np.matmul(np.transpose(V), D_zero), U)

更新：Chrysophylaxs指出了我的错误：我需要转置U：

A_inv = np.matmul(np.matmul(np.transpose(V), D_zero), np.transpose(U))

这个变化解决了问题。非常感谢！

numpy

来源：https://stackoverflow.com/questions/76005798/numpy-svd-does-not-agree-with-r-implementation

1条答案

按热度按时间

j8ag8udp1#

感谢Chrysophylaxs，下面是现在可以正常工作的代码：

# https://stackoverflow.com/questions/75998775/python-vs-matlab-why-my-matrix-is-singular-in-python

import numpy as np

def pseudo_inverse_solver(A, b):
    A_inv = np.linalg.pinv(A)
    x = np.matmul(A_inv, b)
    error = np.matmul(A, x) - b
    return x, error, A_inv

def svd_solver_so(A, b, max_diag=1.0e15):
    """
    see https://stackoverflow.com/questions/24913232/using-numpy-np-linalg-svd-for-singular-value-decomposition
    see https://stackoverflow.com/questions/59292279/solving-linear-systems-of-equations-with-svd-decomposition
    :param A: Matrix in equation A*x = b
    :param b: RHS vector in equation A*x = b
    :param max_diag: max value of diagonal for setting to zero.
    :return: x solution, error vector
    """
    U, D, V = np.linalg.svd(A, full_matrices=False)
    D_diag = np.diag(np.diag(np.reciprocal(D)))
    D_zero = np.array(D_diag)
    D_zero[D_zero >= max_diag] = 0.0
    D_zero = np.diag(D_zero)
    A_inv = V.T @ D_zero @ U.T
    c = U.T @ b
    w = D_zero @ c
    x = V.T @ w
    error = np.matmul(A, x) - b
    return x, error, A_inv

def svd_solver(A, b, max_diag=1.0e15):
    U, D, V = np.linalg.svd(A, full_matrices=False)
    D_diag = np.diag(np.diag(np.reciprocal(D)))
    D_zero = np.array(D_diag)
    D_zero[D_zero >= max_diag] = 0.0
    D_zero = np.diag(D_zero)
    A_inv = np.matmul(np.matmul(np.transpose(V), D_zero), np.transpose(U))
    x = np.matmul(A_inv, b)
    error = np.matmul(A, x) - b
    return x, error, A_inv

if __name__ == '__main__':
    """
    Solution from my SO answer
    https://stackoverflow.com/questions/19763698/solving-non-square-linear-system-with-r/19767525#19767525
    Example showing how to use NumPy SVD
    https://stackoverflow.com/questions/24913232/using-numpy-np-linalg-svd-for-singular-value-decomposition
    """
    np.set_printoptions(20)
    A = np.array([
        [0.0, 1.0, -2.0, 3.0],
        [5.0, -3.0, 1.0, -2.0],
        [5.0, -2.0, -1.0, 1.0]
    ])
    b = np.array([-17.0, 28.0, 11.0]).T

    x_svd, error_svd, A_inv_svd = svd_solver(A, b)
    error_svd_L2 = np.linalg.norm(error_svd)
    x_pseudo, error_pseudo, A_inv_pseudo = pseudo_inverse_solver(A, b)
    error_pseudo_L2 = np.linalg.norm(error_pseudo)
    x_so, error_so, A_inv_so = svd_solver_so(A, b)
    error_so_L2 = np.linalg.norm(error_so)

赞(0）回复(0）举报 2023-04-21

我来回答

NumPy SVD与R实现不一致

1条答案

相关问题

热门标签

最新问答