Random solutions of indefinite linear systems
Consider an indefinite linear system of equations Ax=b
.
I would like to find a set of vectors x_1, ..., x_n
such that they decide everything Ax=b
, and they are as different from each other as possible.
The second part is actually less important; I would be happy with an algorithm that returns a random solution Ax=b
every time I call it.
I know that scipy.sparse.linalg.lsqr
and numpy.linalg.lstsq
return a sparse solution (in terms of least squares) to an underdetermined linear system Ax=b
, but I am not interested in the properties of the solution; I just want any solution Ax=b
if I can create many different solutions.
In fact, scipy.sparse.linalg.lsqr
and numpy.linalg.lstsq
must follow an iterative process that moves from solution to solution until it finds a solution that appears to be least squares. Well, is there a python module that allows me to switch between solutions without any specific purpose and return them?
source to share
For underdetermined system A · x = b , you can calculate the empty space of your matrix coefficients A . The empty space, Z , is the set of basis vectors enclosing the subspace A such that A · Z = 0 . In other words, the column Z are vectors orthogonal to all the rows in A . This means that for any solution x ' up to A x = b , then x' + Z c must also be a solution for any arbitrary vector c .
So, if you want to try random solutions for A x = b , you can do the following:
- Find any solution x ' up to A x = b . You can do this with help
np.linalg.lstsq
, which finds a solution that minimizes the L2 x ' norm . - Find the null space of A . There are several ways to do this, most of which are described in the previous question .
- Pick a random vector c and calculate x ' + Z · c . This will be the solution A x = b .
For example:
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))
Output example:
3.33066907388e-15
3.58036167305e-15
4.63775652864e-15
4.67877015036e-15
4.31132637123e-15
The space of possible c vectors is infinite - you will need to make choices about how you want to try them.
source to share
Here is the code accompanying my comment. It uses a module rank_nullspace.py
from the Scipy Cookbook .
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]))
With a bunch of solutions in there, x
you can choose ones that are "different" from each other, however you define "different".
source to share