Scipy LDL 分解返回意外结果

Question

我生成了一个随机的 5*5 矩阵 x，如下所示：

>>> x = np.random.randn(5,5)

然后像这样使用 scipy.linalg.ldl 分解来分解它：

>>> l, d, p = la.ldl(x)

使用 l、d 和 p 我想 return x。我以为我可以做到以下几点：

>>> l[p,:] @ d @ l[p,:].transpose() - x

但这并没有像我预期的那样给我零。谁能解释我哪里出错了？

我的目标是获得下对角矩阵L使得x = LDL^T不需要行置换矩阵p，但我很困惑scipy 作为输出给出。

Answer 1

LDL 分解算法仅适用于Hermitian/symmetric 矩阵。您正在向它传递一个具有随机值的矩阵，该矩阵不太可能是对称的。另外，应该在不将置换矩阵应用于下三角矩阵的情况下执行矩阵乘法。

将非对称矩阵传递给scipy.linalg.ldl时，仅引用矩阵的下三角或上三角部分，具体取决于lower关键字参数的值，默认为True。我们可以通过 np.isclose():

来查看效果

>>> x = np.random.randn(5,5)
>>> l, d, p = la.ldl(x)
>>> np.isclose(l.dot(d).dot(l.T) - x, 0)
[[ True False False False False]
 [ True  True False False False]
 [ True  True  True False False]
 [ True  True  True  True False]
 [ True  True  True  True  True]]

这里我们看到矩阵的上三角部分被假定为对称的，因此算法返回的值如果是这种情况将是正确的。

下面，我们传la.ldl一个实际的对称矩阵，得到预期的结果。

>>> x = np.array([[1, 2, 3],
                  [2, 4, 5],
                  [3, 5, 6]])
>>> l, d, p = la.ldl(x)
>>> print(np.isclose(l.dot(d).dot(l.T) - x, 0))
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]

如果您正在寻找一般的 LDL^T 分解，没有排列，这会进一步减少矩阵的域。您的矩阵也需要 positive definite.

这是一个这样的矩阵示例：

>>> x = np.array([[2, -1, 0],
                  [-1, 3, -1],
                  [0, -1, 4]])
>>> l, d, p = la.ldl(x)
>>> l
array([[ 1. ,  0. ,  0. ],
       [-0.5,  1. ,  0. ],
       [ 0. , -0.4,  1. ]])
>>> d
array([[2. , 0. , 0. ],
       [0. , 2.5, 0. ],
       [0. , 0. , 3.6]])
>>> p
array([0, 1, 2], dtype=int64)

如您所见，排列 p 只是 [0, 1, 2]，而 l 已经是下三角排列。

Scipy LDL 分解返回意外结果

Scipy LDL decomposition returning unexpected result

python

math

numpy

linear-algebra

scipy