使用NumPy判断矩阵是否正定。

Question

如何判断一个矩阵是否为正定矩阵？

我的矩阵是一个NumPy矩阵。我本来期望在NumPy库中能找到相关的方法，但是没有成功。

- Zygimantas Gatelis

Answer 1

你还可以检查矩阵的所有特征值是否为正。如果是的话，那么这个矩阵就是正定的。

import numpy as np

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

- Akavall

---

你可以使用np.linalg.eigvals代替，它只计算特征值。即使这样，它比@NPE的方法慢得多（对于10x10矩阵来说是3倍，对于1000x1000矩阵来说是40倍）。 - jorgeca

一般来说，所有正特征值都意味着正定性并不正确，除非你知道矩阵是对称的（实数情况）或共轭转置的（复数情况）。例如，A = array([[1, -100],[0, 2]])就不是正定的。有些人可能会将对称或共轭转置作为“正定”定义的一部分，但这并不普遍。 - Warren Weckesser

@WarrenWeckesser 哎呀，对了，不是过于追求细节！实际上，如果使用 np.linalg.cholesky，检查对称性也是必要的（它不会检查并可能返回错误结果，就像您的示例一样）。我想知道如何通过数值方法检查非对称矩阵是否正定... - jorgeca

你可以使用 np.all(x-x.T==0) 来检查对称性。 - shinjin

这非常低效！对于行或列大于6或7的矩阵，请按NPE所指示使用cholesky。虽然cholesky方法可能不太方便（需要捕获异常等），但它要节约得多。 - Alex Flint

如果你正在处理厄米或对称矩阵，建议使用np.linalg.eigvalsh代替；它的计算速度更快。 - Joren

Answer 2

您可以尝试计算Cholesky分解（numpy.linalg.cholesky）。如果矩阵不是正定的，它将引发LinAlgError错误。

- NPE

---

这个方法应该比特征值解法更加高效。 - MRocklin

请注意，在半正定的情况下，从数值上讲，可以向矩阵中添加一个小的单位矩阵（从而将所有特征值移动一小部分，例如几倍机器精度），然后像往常一样使用Cholesky方法。 - jawknee

一个快速的警告：从 numpy.linalg.cholesky 的文档中可以看到：没有执行检查以验证 a 是否是 Hermitian。此外，只使用 a 的下三角和对角线元素。换句话说，如果你测试一个非对称或非 Hermitian 矩阵，即使该矩阵不是正定的，cholesky 也可能不会引发 LinAlgError。请参考 @DanielGarza 的答案，了解潜在的修复方法。 - Praveen

@MRocklin，您所说的"efficient"是指什么？指较低的时间复杂度吗？ - Joren

@MRocklin 你说的"efficient"是指什么？是指更低的时间复杂度吗？ - Joren

Answer 3

以上所有回答中似乎存在一些小混淆（至少与问题相关）。

对于实矩阵，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)，以避免由于浮点误差而导致的差异。

- Daniel Garza

---

挑刺一下：你应该检查 numpy.linalg.LinAlgError，除非你使用 from numpy.linalg import LinAlgError 导入它，但如果我在代码中只捕获这个特定的异常一两次，我不想这样做。 - Imperishable Night

如果使用复杂矩阵，这可能会导致错误（即如果A是复正定的，因此具有严格正的特征值的共轭转置矩阵，则Cholesky技巧仍然是正确的，但它不会通过np.allclose(A, A.T)的第一个“if”语句）。使用np.allclose(A, A.H)将解决此问题（OP表示他正在使用numpy矩阵，如果使用ndarray，请使用A.conj().T而不是.H）。 - H. Rev.

Answer 4

为了用一些现成的代码来说明 @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 

- MarcoMag

Answer 5

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

此外，据说它比Lu分解更具数值稳定性。而Lu分解比找到所有特征值的方法更稳定。

而且，这是一个非常优雅的解决方案，因为事实上： 矩阵具有Cholesky分解，当且仅当它是对称正定的。 那么为什么不使用数学呢？也许有些人害怕异常的出现，但这也是事实，使用异常编程非常有用。

- InfiniteLooper

---

根据相同的维基百科页面，似乎您的说法是不正确的。该页面表明：“如果矩阵A是Hermitian和正半定，则即使允许L的对角线条目为零，它仍然具有形式为A = LL*的分解。” 因此，具有Cholesky分解的矩阵并不意味着矩阵是对称正定的，因为它可能只是半定的。我理解错了吗？此外，似乎您只是在这个蕴涵中添加了“对称”。例如，应该是每个Hermitian正定矩阵都有唯一的Cholesky分解吗？ - user3731622

@user3731622 正确的说，Hermitian半正定矩阵也有Cholesky分解。但是与Hermitian正定情况不同，它并没有唯一的分解。因此，NumPy选择给出错误提示，而不是随机返回可能有无限多种分解的结果。 - CognizantApe

Answer 6

如果你特别需要对称的（如果是复数，则为Hermitian）半正定矩阵，那么下面的代码将可以实现。如果你不关心对称性（如果是复数，则为Hermitian），可以去掉检查它的'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

- CognizantApe

---

@DeepRazi Numpy的Cholesky分解实现适用于复数（即复杂的np.dtype）。因此，在这个意义上它是有效的。但是，我的上面的代码最初检查的是转置是否等于自身，而不是共轭转置，这使得整个函数对于复数无效。我现在已经将转置更改为共轭转置，现在它对于复数是有效的。 - CognizantApe

np.all( X - X.H == 0 ) 这个判断条件太严格了，我认为 np.allclose(X, X.H) 更好。 - syockit

@syockit 如果你只关心矩阵是否近似厄米矩阵而不是实际上的厄米矩阵，那么是的。但要明确的是，因为我们从X中减去元素本身，所以在计算X - X.H时没有浮点数或其他数字误差。因此，如果您以可靠的厄米方式构建X，则此算法是正确的选择。如果您正在逼近一个厄米矩阵，那么我同意使用allclose更好。 - CognizantApe

Answer 7

对于一个实矩阵 $A$，我们有 $x^TAx=\frac{1}{2}(x^T(A+A^T)x)$，而 $A+A^T$ 是对称实矩阵。因此，当且仅当 $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)

- Martin Wang

---

这是唯一一个正确回答OP问题的答案：“如何确定矩阵是否为DP”。所有其他答案都令人困惑地假设对称性是矩阵为正定所必需的，但事实并非如此。 - Guillaume

Answer 8

正定

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+1e-5*np.identity(n_size) # little trick I found for PS matrix

- silgon

---

使用eigvals检查正定性时，您可能需要将其与某个epsilon进行比较。例如，np.linalg.matrix_rank使用机器epsilon乘以最大特征值（和max(x.shape)），而np.linalg.lstsq以类似的方式使用机器epsilon。实际上，如果存在低于机器epsilon但带有负号的特征值，则上述PSD测试可能会失败。 - syockit

你的评论很到位（我还注意到代码中有一个错误）。当你处理这些情况时，你会经常遇到你所谈论的情况，就像我的情况是由于一些求解器的容差率。大多数情况下，我的解决方法是将一个单位矩阵乘以一个小值相加，就像答案中的最后一行代码所示，但这当然取决于你的公差，在我的情况下没问题。 - silgon

Answer 9

对于非对称矩阵，您可以使用主子式测试：

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")

- Pietro Bonazzi

Answer 10

我可以建议这个解决方案，适用于非对称矩阵。基本上，它尝试找到一个非零向量z，使得z^T · M · z的结果最小化。正如您所看到的，它可以收敛于最小值(minV['success'])，或达到最大迭代次数(minV['status'] == 2)。如果在这一点上结果仍然是正的，我们可以认为该矩阵是正定的。

我想可能有分析方法来确定它，但我看到这里关于对称矩阵的一些混淆(不是正定的前提!!)。

from scipy.optimize import minimize
import numpy as np

def is_pos_def(m,z0='rand', maxiter=100):
    #tells if matrix is positive definite
    if m.shape[0] != m.shape[1]:
        raise Exception("Matrix is not square") 
    elif (m==m.T).all(): #symmetry testing
        return np.all(np.linalg.eigvals(m) > 0)
    else:
        def f(z):
            z=np.array(list(z))[:,np.newaxis]
            return np.dot(np.dot(z.T, m),z)
        if z0=='rand':
            z0 = list(np.random.rand(m.shape[0]))
        #constraints for a non-zero vector solution
        cons = ({'type': 'ineq', 'fun': lambda z:  np.sum(np.abs(z))})
        minV = minimize(f, z0, method='COBYLA', options={'maxiter' : maxiter},constraints=cons);

        if minV['success'] or minV['status'] == 2:
            return minV['fun']+0 > 0 
        else:        
            return minV

这种方法适用于对称和非对称的矩阵，你可以使用以下矩阵进行测试（也可以在wolfram alpha上验证）

m=np.array([[3, 0, 0],
            [0, 2, 0],
            [4, 3, 3]])

- porygon-tech