待定线性系统的随机解

Question

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

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

第二部分其实没那么重要；我会很高兴有一种算法，每次我调用它时 return 都是 Ax=b 的随机解。

我知道 scipy.sparse.linalg.lsqr 和 numpy.linalg.lstsq return 是欠定线性系统 Ax=b 的稀疏解（根据最小二乘法），但我不知道关心解决方案的属性；我只要Ax=b的任意解，只要能生成一堆不同的解即可

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

Answer 1

这是我的评论附带的代码。它使用 Scipy Cookbook.

中的 rank_nullspace.py 模块

import numpy as np
from numpy.linalg import lstsq

from rank_nullspace import nullspace
# rank_nullspace from
# http://scipy-cookbook.readthedocs.io/items/RankNullspace.html


def randsol(A, b, num=1, check=False):
    xLS, *_ = lstsq(A, b)

    colsOfNullspace = nullspace(A)
    nullrank = colsOfNullspace.shape[1]
    if check:
        assert(np.allclose(np.dot(A, xLS), b))
        assert(np.allclose(np.dot(A, xLS + np.dot(colsOfNullspace,
                                                  np.random.randn(nullrank))),
                           b))

    sols = xLS[:, np.newaxis] + np.dot(colsOfNullspace,
                                       np.random.randn(nullrank, num))
    return sols


A = np.random.randn(2, 10)
b = np.random.randn(2)
x = randsol(A, b, num=50, check=True)
assert(np.allclose(np.dot(A, x), b[:, np.newaxis]))

在 x 中有一堆解决方案，您可以选择彼此“不同”的解决方案，但是您定义“不同”。

Answer 2

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

因此，如果您想对 A·x = b[ 的随机解进行采样=84=] 那么你可以这样做：

求 x' 到 A·x = b。您可以使用 np.linalg.lstsq 来执行此操作，它会找到一个解决方案，使 x'.
的 L2 范数最小化
找到 A 的 null space。有许多不同的方法可以做到这一点，其中大部分都在 this previous question.
中介绍。
采样一个随机向量c，然后计算x' +Z·c。这将是A·x = b[=85=的解].

例如：

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 向量的 space 是无限的 - 你必须做出一些选择想尝尝这些。

待定线性系统的随机解

Random solutions of undetermined linear systems

python

numpy

linear-algebra

scipy