未确定线性系统的随机解

Question

考虑一个不确定的线性方程组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模块可以让我在没有特定目标的解决方案之间跳转，并返回它们？

Answer 1

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

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

1. 找到a·x =b的任何解x'。你可以使用np.linalg.lstsq来做到这一点，它找到了一个最小化x'的 L2 范数的解决方案。
2. 找到A的空空间。有许多不同的方法可以做到这一点，其中大多数都在上一个问题中介绍。
3. 对随机向量 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向量的空间是无限的 - 你必须选择如何采样这些向量。

Answer 2

这是我评论附带的代码。它使用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，您可以选择彼此"不同"的解决方案，无论您如何定义"不同"。