np.reciprocal计算异常：无关操作或打印影响输出问题咨询

问题现象 #

在实现SVD伪逆时，调用np.reciprocal(S, where=(S != 0))会出现异常：当函数中存在无关的U矩阵计算时，伪逆矩阵的非对角元素被意外填充了随机值；但移除U的计算或在求逆前打印S时，结果又恢复正常。

根本原因 #

这不是numpy的bug，而是np.reciprocal结合where参数的行为细节导致的：

• 当使用where参数且未指定out时，numpy会创建一个新数组，但该数组的内存可能复用了之前已释放的对象内存（比如ComputeSVD函数中临时创建的U数组，函数返回后U被销毁，内存被放回numpy的内存池）。
• 未满足where条件的位置（即S中的0元素）不会被np.reciprocal修改，因此会保留内存池中的旧值（也就是U数组的残留数据）。
• 移除U的计算或打印S时，内存布局被改变，新数组不会复用U的残留内存，因此非对角元素保持为0。

解决方案 #

避免依赖np.reciprocal的where参数处理稀疏对角矩阵，改用显式构造伪逆矩阵的方式，确保非对角元素始终为0：

方案1：直接操作对角元素 #

import numpy as np
from numpy.linalg import eig

def ComputeSVD(A):
    m, n = A.shape
    eigenvalues, eigenvectors = eig(A.T @ A)

    Sigma = np.zeros((m,n))
    np.fill_diagonal(Sigma, np.sqrt(eigenvalues))
    print(Sigma)

    U = np.zeros((m,m))
    for i in range(len(eigenvalues)):
         U[:,i] =  A @ (eigenvectors[:,i]) 

    return Sigma


def PseudoInverse(A):
    S = ComputeSVD(A)
    S_inv = np.zeros_like(S)
    # 获取对角元素的索引
    diag_idx = np.diag_indices(min(S.shape[0], S.shape[1]))
    # 仅对对角元素求逆
    S_inv[diag_idx] = np.reciprocal(S[diag_idx])
    return S_inv

A = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 3.0], [1.0, 4.0, 9.0]])
print(PseudoInverse(A))

方案2：重构对角矩阵 #

def PseudoInverse(A):
    S = ComputeSVD(A)
    # 提取对角元素并求逆
    s_diag = np.diag(S)
    s_diag_inv = np.reciprocal(s_diag, where=s_diag != 0)
    # 用逆后的对角元素重构矩阵
    S_inv = np.diag(s_diag_inv)
    # 若原矩阵不是方阵，补0
    if S.shape != s_diag_inv.shape:
        S_inv = np.pad(S_inv, ((0, S.shape[0]-s_diag_inv.shape[0]), (0, S.shape[1]-s_diag_inv.shape[1])), mode='constant')
    return S_inv

两种方案都能确保伪逆矩阵的非对角元素始终为0，不受临时变量内存复用的影响。

内容的提问来源于stack exchange，提问作者Thomas Beck