python - 在python中测试矩阵是否为半正定(PSD)时出错

Question

我正在使用 python 2，我想测试矩阵是否为半正定 (PSD) 矩阵。 我建立了一个随机矩阵 X，我想测试 Q = X^TX 的 SDP 属性。

为此，我调整了一个测试正定属性的函数，以便测试半正定属性。然后我计算矩阵 Q = X^TX。 X^TX 应该是 PSD，因为它是矩阵及其转置的乘积(有关 PSD 保证的更多详细信息，请参见 Is a matrix multiplied with its transpose something special?)。但是，当我测试 Q 是否为 PSD 时，该函数返回 false。

有人知道我的代码哪里错了吗？是我的测试功能不测试 PSD 属性还是其他？

这是我的脚本(比实际更简单):

from scipy.stats import bernoulli
from scipy import linalg
import numpy as np

p = 300
N = 100
np.random.seed(18) 
X = bernoulli.rvs(0.5, size=p*N).reshape((N, p))
X = 2 * X - 1* np.ones_like(X)
Q = np.dot(X.T, X)
def is_semi_pos_def(x):
    return np.all(np.linalg.eigvals(x) >= 0)
is_semi_pos_def(Q)

它返回:

False

非常感谢您对此提供的任何帮助。

Answer 1

如果 N < p，您在脚本中创建的矩阵 Q 没有满秩，一些特征值将为 0。但是正如评论中提到的，np.eigvals 中的数值错误导致了这些变成负数:

P = 5
N = 3
np.random.seed(18) 
X = bernoulli.rvs(0.5, size=N*P).reshape((N, P))
X = 2 * X - 1* np.ones_like(X)
Q = np.dot(X.T, X)
np.linalg.eigvals(Q)

返回

[1.24244289e+01,  -2.96746135e-16,   2.57557110e+00, 4.23704588e-33,   4.25752762e-18]

您可以通过测试特征值是否大于 0 - E 来解释错误，对于一些适当小的 E > 0。但是，如果在实践中您知道您的矩阵将是满秩的，您可以计算 Cholesky decomposition这似乎是使用 numpy 的最快方法。

import numpy as np

P = 300
N = 100
X = bernoulli.rvs(0.5, size=N*P).reshape((N, P))
X = 2 * X - 1* np.ones_like(X)
Q = np.dot(X.T, X)

@timing
def is_semi_pos_def_chol(x):
    try:
        np.linalg.cholesky(x)
        return True
    except np.linalg.linalg.LinAlgError:
        return False

@timing
def is_semi_pos_def_eigsh(x, epsilon=1e-10):
    return np.all(np.linalg.eigvalsh(x) >= -epsilon)

@timing
def is_semi_pos_def_eigs(x, epsilon=1e-10):
    return np.all(np.linalg.eigvals(x) >= -epsilon)

print is_semi_pos_def(Q)
print is_semi_pos_def_eigs(Q)
print is_semi_pos_def_eigsh(Q)

返回:

is_semi_pos_def function took 0.001 s
False
is_semi_pos_def_eigs function took 0.073 s
True
is_semi_pos_def_eigsh function took 0.011 s
True

我们可以使用更快的 eigvalsh 因为矩阵是对称的。 请注意，由于 Q 的非平凡内核，Cholesky 分解错误地返回 False。但是对于 N=P=300，这不是问题:

is_semi_pos_def_chol function took 0.002 s
True
is_semi_pos_def_eigs function took 0.118 s
True
is_semi_pos_def_eigsh function took 0.023 s
True

如果您想要一个分析上合理的答案，您可以使用 sympy 来计算没有数值错误的特征值，或者使用 sympy 来计算 Cholesky 分解，这也可以避免 float 错误。然而，对于大的 N 和 P，计算时间变得令人望而却步:

from sympy.mpmath import mp

P = 30
N = 10
X = bernoulli.rvs(0.5, size=N*P).reshape((N, P))
X = 2 * X - 1* np.ones_like(X)
M = Matrix(np.dot(X.T, X))
Q = np.dot(X.T, X)


@timing
def is_semi_pos_def_symbolic(x):
    try:
        M.cholesky()
        return True
    except ValueError as e:
        print e
        return False

print is_semi_pos_def_symbolic(M)
print is_semi_pos_def_chol(Q)
print is_semi_pos_def_eigsh(Q)

注意分析真相的巨大成本:

is_semi_pos_def_symbolic function took 0.908 s
True
is_semi_pos_def_chol function took 0.000 s
False
is_semi_pos_def_eigsh function took 0.000 s
True