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

Question

我正在尝试通过 SciPy 的 Python 绑定到 LAPACK 使用 dsysv 例程获得给定对称矩阵的 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，其中包含我感兴趣的输出（参见等式 3-23）。

Answer 1

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

sytrf 的包装器已于 9 月中旬添加到主分支，就在 1.0.0 Beta release 之前。您可以在 Github.

上找到相关的 pull-request and commit

Answer 2

所有需要的信息都可以在 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）。我将在下面尝试解释。首先引用文档：

... additional row interchanges are required to recover U or L explicitly (which is seldom necessary).

重构 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])。这意味着在执行操作之前，第 i 个和第 k 个 row/column 在 A 中被交换了。在这种情况下，lult 的第 k 个对角线元素对应于 D 的第 k 个对角线元素。 L(k) 包含交换后下对角矩阵的第 k 列。

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

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

我不会责怪任何人不理解我的解释。我只是在想出来的时候把它写下来……希望代码片段、描述和原始文档的组合证明是有用的:)

Answer 3

dsysv 是线性系统求解器，它在内部执行所有的魔术操作，包括对 dsytrf 的调用。所以对于因式分解来说，它是不需要的。正如 kazemakase 提到的，现在可以在 SciPy（PR 7941 and will appear officially in version 1.1) and you can just use the scipy.linalg.ldl() 中获取外部因子的因式分解和排列信息。实际上这就是 ?sytrf 和 ?hetrf 的原因已添加。

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

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

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

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

LDLt factorization using SciPy's Python bindings to LAPACK

scipy

lapack