使用 SciPy 的 Python 绑定到 LAPACK 的 LDLt 分解

Question

我正在尝试使用dsysv 例程通过 SciPy 与 LAPACK 的 Python 绑定来获取给定对称矩阵的 LDLt 分解，该例程实际上使用此矩阵分解求解线性系统。

我尝试了以下方法：

import numpy as np
from scipy.linalg.lapack import dsysv

A = np.random.randint(1, 1000, size=(5, 5))
A = (A + A.T)
b = np.random.randn(5)
lult, piv, x, _ = dsysv(A, b, lower=1)

x 将是上述线性系统的解决方案，lult 和 piv 包含有关分解的信息。

如何从中重建 LDLt？有时负值包含在piv 和docs 中，我无法理解它们的含义。

LAPACK 的 sytrf 实际上计算了这个因式分解（没有求解任何线性系统），但它似乎无法通过 SciPy 获得。

有一个示例 here 与我感兴趣的输出（见 eq. 3-23）。

Answer 1

所有必需的信息都可以在documentation of systrf 中找到。但不可否认，这有点冗长。

所以只要给我代码：

import numpy as np
from scipy.linalg.lapack import dsysv

def swapped(i, k, n):
    """identity matrix where ith row and column are swappend with kth row and column"""
    P = np.eye(n)
    P[i, i] = 0
    P[k, k] = 0
    P[i, k] = 1
    P[k, i] = 1
    return P


# example

n = 5

A = np.random.rand(n, n)
A = (A + A.T)
b = np.random.randn(n)
lult, piv, x, _ = dsysv(A, b, lower=1)


# reconstruct L and D

D = np.zeros_like(A, dtype=float)
L = np.eye(n)

k = 0
while k < n:
    i = piv[k]

    if i < 0:
        s = 2
    else:
        s = 1

    if s == 1:
        i = i - 1
        D[k, k] = lult[k, k]  # D(k) overwrites A(k,k)
        Pk = swapped(k, i, n)
        v = lult[k+1:n, k]  # v overwrites A(k+1:n,k)
        Lk = np.eye(n)
        Lk[k+1:n, k] = v
    else:
        m = -i - 1
        D[k:k+2, k:k+2] = lult[k:k+2, k:k+2]  #  the lower triangle of D(k) overwrites A(k,k), A(k+1,k), and A(k+1,k+1)
        D[k, k+1] = D[k+1, k]  # D is symmeric
        Pk = swapped(k+1, m, n)
        v = lult[k+2:n, k:k+2]  # v overwrites A(k+2:n,k:k+1)
        Lk = np.eye(n)
        Lk[k+2:n, k:k+2] = v   

    L = L.dot(Pk).dot(Lk)

    if s == 1:
        k += 1
    else:
        k += 2

print(np.max(np.abs(A - L.dot(D).dot(L.T))))  # should be close to 0

上面的片段从分解中重构 L 和 D（它需要适应从 UDUt 分解中重构 U）。我将在下面尝试解释。首先引用文档：

... 需要额外的行交换来显式恢复 U 或 L（这很少需要）。

重构 L（或 U）需要多次迭代，包括行交换操作和矩阵乘法。这不是很有效（在 Python 中效率较低），但幸运的是这种重建很少需要。因此，请确保您确实必须这样做！

我们从L = P(1)*L(1)* ... *P(k)*L(k)*..., 重构L。 （Fortran 索引从 1 开始）。所以我们需要将k从0迭代到n，每一步得到K和L，然后相乘。

P 是一个置换矩阵，由piv 定义。 piv 的正值是直截了当的 (i = piv[k])。这意味着在执行操作之前，在 A 中交换了第 i 行/第 k 行/列。在这种情况下，lult 的第 k 个对角元素对应于 D 的第 k 个对角元素。 L(k) 包含下对角矩阵的第 k 列 - 交换后。

piv为负值表示D对应的元素是一个2x2的块而不是一个元素，L(k)对应下对角矩阵的两列。

现在对于k 中的每个步骤，我们获得L(k)，应用交换操作P(k)，并将其与现有的L 结合。我们还得到了 D 的 1x1 或 2x2 块，并相应地将k 增加 1 或 2 用于下一步。

---

我不会责怪任何人不理解我的解释。我只是把它写下来，因为我想出来了......希望代码 sn-p、描述和原始文档的组合证明是有用的 :)

Answer 2

dsysv 是线性系统求解器，它在内部完成所有魔法，包括对dsytrf 的调用。因此，对于因式分解，它是不需要的。正如 kazemakase 提到的，这现在可以在 SciPy 中使用（PR 7941 并将在 1.1 版中正式出现），您可以使用scipy.linalg.ldl() 来获取外部因子的分解和排列信息。其实这就是?sytrf和?hetrf添加的原因。

您可以查看它的源代码以了解如何清理ipiv。

在 Windows 10 机器上使用 OpenBlas 构建的 SciPy v.1.1 与使用 mkl 的 matlab 相比，性能如下所示

在其之上添加额外的 JIT 编译器可能会使其达到 matlab 速度。由于 ipiv 处理和分解构造是在纯 numpy/python 中完成的。如果性能是最重要的，或者更好地对其进行 cythonize。

Answer 3

将 scipy 更新到 >= 1.0.0 应该可以解决问题。

9 月中旬，sytrf 的包装器已添加到主分支，就在 1.0.0 Beta release 之前。 你可以在 Github 上找到相关的pull-request 和commit。