python - 未定线性系统的随机解

考虑 欠定 线性方程组 Ax=b 。

我想找到一组向量 x_1, ..., x_n 以便它们都解决 Ax=b 并且它们彼此之间尽可能不同。

第二部分实际上不那么重要；我会对每次调用它时返回 Ax=b 的随机解的算法感到满意。

我知道 scipy.sparse.linalg.lsqr 和 numpy.linalg.lstsq 返回一个欠定线性系统 Ax=b 的稀疏解(就最小二乘而言)，但我不关心解的属性；我只想要 Ax=b 的任何解决方案，只要我可以生成一堆不同的解决方案。

事实上，scipy.sparse.linalg.lsqr 和 numpy.linalg.lstsq 应该遵循一个迭代过程，从一个解决方案跳到另一个解决方案，直到他们找到一个在最小二乘方面似乎是最小值的解决方案。那么，是否有一个 python 模块可以让我在没有特定目标的解决方案之间跳转，并返回它们？

最佳答案

对于欠定系统 A·x = b 您可以计算系数矩阵 A 的 null space 。零空间 Z 是一组跨越 A 子空间的基向量，使得 A·Z = 0 。换句话说， Z 的列是与 A 中的所有行正交的向量。这意味着对于 x' 到 A·x = b 的任何解，那么 x' + Z·c 也必须是任意向量 的解。

因此，如果您想对 A·x = b 的随机解进行采样，那么您可以执行以下操作:

找到 x' 到 A·x = b 的任何解。您可以使用 np.linalg.lstsq 来做到这一点，它会找到一个最小化 x' 的 L2 范数的解决方案。

找到 A 的零空间。有许多不同的方法可以做到这一点，其中大部分都包含在 this previous question 中。

采样一个随机向量 c ，并计算 x' + Z·c 。这将是 A·x = b 的解。

例如:

import numpy as np
from scipy.linalg import qr


def qr_null(A, tol=None):
    """Computes the null space of A using a rank-revealing QR decomposition"""
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    tol = np.finfo(R.dtype).eps if tol is None else tol
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol)
    return Q[:, rnk:].conj()


# An underdetermined system with nullity 2
A = np.array([[1, 4, 9, 6, 9, 2, 7],
              [6, 3, 8, 5, 2, 7, 6],
              [7, 4, 5, 7, 6, 3, 2],
              [5, 2, 7, 4, 7, 5, 4],
              [9, 3, 8, 6, 7, 3, 1]])
b = np.array([0, 4, 1, 3, 2])

# Find an initial solution using `np.linalg.lstsq`
x_lstsq = np.linalg.lstsq(A, b)[0]

# Compute the null space of `A`
Z = qr_null(A)
nullity = Z.shape[1]

# Sample some random solutions
for _ in range(5):
    x_rand = x_lstsq + Z.dot(np.random.rand(nullity))
    # If `x_rand` is a solution then `||A·x_rand - b||` should be very small
    print(np.linalg.norm(A.dot(x_rand) - b))

示例输出:

3.33066907388e-15
3.58036167305e-15
4.63775652864e-15
4.67877015036e-15
4.31132637123e-15

可能的 c 向量的空间是无限的——您必须对如何对这些向量进行采样做出一些选择。

关于python - 未定线性系统的随机解，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/43423346/