How to solve Ax=0 with non-zero solution?-CodePudding

I tried to find a non-zero solution for Ax=0 using 'numpy.linalg.solve', but this package only gives me some solutions with zero vector.

So I tried to search related code as below:

import sympy
import numpy as np
 
A = np.array([[0.8, 0.1, 0.1],
              [0.7, 0.2, 0.1],
              [0.6, 0.3, 0.1]]) - np.eye(3)
 
dim = A.shape[0]
x = sympy.symbols([f"x{idx}" for idx in range(dim)])

# Create a list of unknowns
gen_sol = sympy.solve(np.array(x) @ A, *x)

# to obtain the general solution for Ax = 0
cond = sympy.Eq(sum(x), 1)

# Creating normalization conditions
equals = [sympy.Eq(key, value) for key, value in gen_sol.items()]   [cond]

# Create a system of equations with a general solution, with additional conditions
part_sol = sympy.solve(equals, x)

# Obtaining special solutions for specific conditions
assert part_sol, "Ax = 0 No solution under this condition"
result = np.array([part_sol[key] for key in x])
print(result)

the result is [0.766666666666667 0.133333333333333 0.100000000000000], then A.dot(result) is not a zero-vector which is quite strange.

Would someone help me with that? Or is there any another possible solution for the Ax=0 question?

CodePudding user response：

In your example, x=[0,0,0]' is the solution. np.linalg.solve works, try with another system.

>>> A = np.array([[-.2,.1,.1],
                  [.7,-.8,.2],  # changed .1 -> .2 to get rank 3
                  [.6,.3,-.9]])
>>> b = [0,0,1]  # changed to get non-zero solution
>>> x = np.linalg.solve(A, b)
>>> x
array([8.33333333, 9.16666667, 7.5       ])
>>> A @ x
array([-1.1379786e-16,  2.1649349e-16,  1.0000000e 00])

I have following remark. Your matrix has only rank 2, i.e. the rows are linearly dependent. You should really look first at what A is (e.g., rank, invertability, etc.) before blindly applying solve() method. Let me elaborate.

Solving a linear system

An analytical approach to solve a linear system typically consists

transform a general matrix A into certain form, and
apply a fast algorithm solving Ax=0 with certain assumptions about A. The assumptions are fulfilled by the preceding transformation of A.

A classical high-school approach is Gaussian elimination, which transform A into upper triangular form. After that, a back-substitution solves the system.

For A with higher ranks, Gaussian elimination is inefficient. There are other ways - Cholesky decomposition, LT codes, algorithms for sparse A, etc.

The linear systems can be, and often are, solved numerically (e.g., Gauss-Seidl method).

In case you can apriori assume something about your matrix A, e.g., it is a covariance matrix, hence, it is positive definite, then you might avoid calling solve() and just perform specialized method (e.g., Cholesky decomposition).

Online resources, further elaborating on this topic, are

CodePudding user response：

First of all, you can easily see that an obvious solution is [1,1,1]. (The sum of columns of your matrix is 1)

That's not an answer, but something to keep in mind: we want to find this solution.

That being said, your code works fine (Thank you, btw, I wasn't at all familiar with this usage of sympy, that was interesting to read).

The only problem is that you are solving the wrong equation.

You are not solving AX=0. But AᵀX=0. Because the equation you gave was np.array(x) @ A. x is a 1D array. Which means that it is treated as a line or as a column, depending on what makes sense.

A@x is AX (A times column X), with the result in a 1D array. x@A is XᵀA (that is the line X times A), also with result in a 1D array. So, long story short, np.array(x)@A is

             /-0.2  0.1  0.1 \
[x₁ x₂ x₃] × | 0.7 -0.8  0.1 |
             \ 0.6  0.3 -0.9 /

What you wanted to solve is

/-0.2  0.1  0.1 \   / x₁ \
| 0.7 -0.8  0.1 | × | x₂ |
\ 0.6  0.3 -0.9 /   \ x₃ /

Since u×v is (vᵀuᵀ)ᵀ, the equation your are solving begin xᵀA = 0 <=> (Aᵀxᵀᵀ)ᵀ =0 <=> Aᵀx=0

As you can see, ``` A.T.dot([0.766666666666667,0.133333333333333,0.100000000000000])

is indeed 0. So your code works. Just not solving the correct equation.


So, at last, the answer (I am sure you've got it yourself by now, but I conclude anyway); line 
```python
gen_sol = sympy.solve(np.array(x) @ A, *x)

should be

gen_sol = sympy.solve(A.dot(x), *x)

Or,

gen_sol = sympy.solve(np.array(x) @ A.T, *x)

as you wish Then you get the result x0=x1, x1=x2, which combine with your x0 x1 x2=1, gives, with your code, [0.333333333333, 0.333333333333, 0.3333333333333], as expected.

Edit

Note: (but keep in mind I am discovering, through your code, the symbolic computation capabilities of sympy — I've known symbolic computation for a while, but I've never tried it with sympy, even tho, I wanted to try since a while. Reason why I take advantage of your question to play a little bit with it) I think that you're usage might be too convoluted. I mean, you solve the equation Ax=0. And then you solve the result of this equation, with equality system.

Why not do that in 1 step?

eq=A@x # Equations
equals=[sympy.Eq(eq[i], 0) for i in range(dim)]   [cond] 
part_sol = sympy.solve(equals, x)

Second Edit:

You have also a more direct, more efficient (but let interesting from my point of view) way to find the result

import scipy.linalg
import numpy as np
A = np.array([[0.8, 0.1, 0.1],
              [0.7, 0.2, 0.1],
              [0.6, 0.3, 0.1]]) - np.eye(3)
result=scipy.linalg.null_space(A)

it return

array([[-0.57735027],
       [-0.57735027],
       [-0.57735027]])

(1 vector because null space is 1D. It could have been 2 vectors in matrix rank were 1. Or even 3 vectors if A had been 0)

Less fun, tho. But, if you were not aware of null_space, maybe you needed to know :)

Third Edit:

(Just because I wonder what other, more numerical way there could be, without symbolic computation — null_space is not symbolic, but well, it is part of scipy, so maybe it is a little under the hood) You could also compute eigen values and vectors. And the solutions you are looking for are the eigen vectors associated to the 0 eigen value.

In your case

import numpy as np
# See, no import scipy
A = np.array([[0.8, 0.1, 0.1],
              [0.7, 0.2, 0.1],
              [0.6, 0.3, 0.1]]) - np.eye(3)
val,vec=np.linalg.eig(A)
result=[vec[:,i] for i in range(len(A)) if np.abs(val[i])<1e-14]