用numpy找出矩阵是否是正定的

Question

我需要确定矩阵是否为positive definite。我的矩阵是 numpy 矩阵。我期待在 numpy 库中找到任何相关的方法，但没有成功。 感谢您的帮助。

Answer 1

您可以尝试计算 Cholesky 分解 (numpy.linalg.cholesky)。如果矩阵不是正定矩阵，这将引发LinAlgError。

Answer 2

你还可以检查矩阵的所有特征值是否为正，如果是则矩阵是正定的：

import numpy as np

def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 0)

Answer 3

我不知道为什么 NPE 的解决方案被低估了。这是最好的方法。我在Wkipedia 上发现复杂度是立方的。

此外，据说它在数值上比 Lu 分解更稳定。并且Lu分解比寻找所有特征值的方法更稳定。

而且，这是一个非常优雅的解决方案，因为这是事实：

矩阵具有 Cholesky 分解当且仅当它是对称正的。

那么为什么不使用数学呢？也许有些人害怕引发异常，但这也是一个事实，使用异常进行编程非常有用。

Answer 4

对于一个实数矩阵 $A$，我们有 $x​​^TAx=\frac{1}{2}(x^T(A+A^T)x)$，$A+A^T$ 是对称实矩阵。所以 $A$ 是正定的，如果 $A+A^T$ 是正定的，如果 $A+A^T$ 的所有特征值都是正的。

import numpy as np

def is_pos_def(A):
    M = np.matrix(A)
    return np.all(np.linalg.eigvals(M+M.transpose()) > 0)

Answer 5

用一些现成的代码来说明@NPE的答案：

import numpy as np

def is_pd(K):
    try:
        np.linalg.cholesky(K)
        return 1 
    except np.linalg.linalg.LinAlgError as err:
        if 'Matrix is not positive definite' in err.message:
            return 0
        else:
            raise 

Answer 6

以上所有答案似乎都存在一些小混乱（至少与问题有关）。

对于实矩阵，np.linalg.cholesky 中的正特征值和正前导项的检验仅适用于对称矩阵。所以首先需要测试矩阵是否对称，然后应用其中一种方法（正特征值或 Cholesky 分解）。

例如：

import numpy as np

#A nonsymmetric matrix
A = np.array([[9,7],[6,14]])

#check that all eigenvalues are positive:
np.all(np.linalg.eigvals(A) > 0)

#take a 'Cholesky' decomposition:
chol_A = np.linalg.cholesky(A)

矩阵 A 不是对称的，但特征值是正的，Numpy 返回的 Cholesky 分解是错误的。你可以检查一下：

chol_A.dot(chol_A.T)

与 A 不同。

您还可以检查上面的所有 python 函数是否会为“正定性”测试为阳性。如果您尝试使用 Cholesky 分解来计算逆，这可能是一个严重的问题，因为：

>np.linalg.inv(A)
array([[ 0.16666667, -0.08333333],
   [-0.07142857,  0.10714286]])

>np.linalg.inv(chol_A.T).dot(np.linalg.inv(chol_A))
array([[ 0.15555556, -0.06666667],
   [-0.06666667,  0.1       ]])

不同。

总之，我建议在上面的任何函数中添加一行来检查矩阵是否对称，例如：

def is_pos_def(A):
    if np.array_equal(A, A.T):
        try:
            np.linalg.cholesky(A)
            return True
        except np.linalg.LinAlgError:
            return False
    else:
        return False

您可能希望将上述函数中的 np.array_equal(A, A.T) 替换为 np.allclose(A, A.T) 以避免因浮点错误导致的差异。

Answer 7

如果您特别想要对称（厄米特，如果复杂）正半定矩阵，则可以使用以下方法。如果您不关心对称性（厄米特，如果复杂），请删除检查它的“if”状态。如果您想要正定而不是正半定而不是删除正则化行（并将传递给 'np.lingalg.cholesky()' 的值从 'regularized_X' 更改为 'X'）。下面

import numpy as np

def is_hermitian_positive_semidefinite(X):
    if X.shape[0] != X.shape[1]: # must be a square matrix
        return False

    if not np.all( X - X.H == 0 ): # must be a symmetric or hermitian matrix
        return False

    try: # Cholesky decomposition fails for matrices that are NOT positive definite.

        # But since the matrix may be positive SEMI-definite due to rank deficiency
        # we must regularize.
        regularized_X = X + np.eye(X.shape[0]) * 1e-14

        np.linalg.cholesky(regularized_X)
    except np.linalg.LinAlgError:
        return False

    return True

Answer 8

对于非对称矩阵，您可以使用 Principal Minor Test ：

def isPD(Y):
  row = X.shape [0]
  i = 0
  j = 0
  for i in range(row+1) : 
    Step = Y[:i,:j]
    j+=1
    i+=1
    det = np.linalg.det(Step)
    if det > 0 : 
        continue 
    else :
        return ("Not Positive Definite, Test Principal minor failed")

  return ("Positive Definite")

Answer 9

正定

numpy.linalg.cholesky(x) # just handle the error LinAlgError

正半定

np.all(np.linalg.eigvals(x) >= 0)

注意：大多数情况下np.all(np.linalg.eigvals(x) > 0) 也会给你，如果你的矩阵是PSD，即使你看到> 而不仅仅是>=，我几天前就遇到了这个问题。我认为这应该与舍入误差有关，因为我们的特征值非常小，甚至 Cholesky 分解也可能产生错误。

注意

为了测试，您可能需要创建一些正半定矩阵和一些正定矩阵：

n_size=4
a = np.random.rand(n_size)
A_PSD = np.outer(a,a)  # the outer product of any vector generates a PSD matrix
A_PD = A_PSD+np.identity(n_size) # little trick I found for PS matrix