CODEWITHSUNDEEP
pythonnumpymatrixlinear-algebra
tobias hassebrock
tobias hassebrock

Reputation: 357

Solve non square matrix with python: How to use numpy.linalg.lstsq()?

Requested behaviour
I would like to solve a non-square matrix with python. The matrix has two linearly dependent vectors.

Current State
I tried to use numpy.linalg.solve() first, but that only works for square matrices. Other StackOverflow posts recommended using numpy.linalg.lstsq().

issue
However, I do not understand how to implement numpy.linalg.lstsq() correctly. The function solves the last parameter correctly, but not the other parameters. One post recommendet this solution which I do not understand either.

Do I have to implemt a loop somehow?

Can someone provide me with a code example? How do solve this matrix problem using python?

My current Code

import numpy as np

# defining a linear equation system E=F with
#   | -2 * x1 - 4 * x2 + 1 * x3 -  9 * x4 + 0 * x5 =  +5  |
#   |  3 * x1 + 6 * x2 + 0 * x3 + 12 * x4 + 3 * x5 = +15  |
#   |  1 * x1 + 2 * x2 + 1 * x3 +  3 * x4 + 1 * x5 = -17  |
#   | -5 * x1 - 4 * x2 + 1 * x3 -  9 * x4 + 0 * x5 = +14  |


E=np.array(
    [
        [-2,-4,1,-9,0],
        [3,6,0,12,3],
        [1,2,1,3,1],
        [-5,-10,3,-23,1]
    ]
)

F=np.array(
    [3,15,-17,14]
)

solutionNonSquare = np.linalg.lstsq(E, F)
print('the solution vector is: {x1, x2, x3, x4, x5}=')
print(solutionNonSquare)

Written Matrix solution enter image description here

Upvotes: 2

Views: 5800

Answers (1)

Subhaneil Lahiri
Subhaneil Lahiri

Reputation: 458

This is an underdetermined system of equations. This means that there are many solutions, and there is no such thing as "the" solution. The fact that Gaussian elimination and lstsq give different solutions does not mean that anything is wrong.

Let's generate and check various solutions:

import scipy.linalg as sla

E_null = sla.null_space(E)

def check_solution(coeffs):
    x = solutionNonSquare[0] + E_null @ coeffs
    check = E @ x - F
    with np.printoptions(precision=2, suppress=True):
        print('x = {}'.format(x))
    with np.printoptions(precision=5, suppress=True):
        print('E . x - F = {}'.format(check))
    print('|x| = {}'.format(np.linalg.norm(x)))

We can check the minimum norm solution produced by lstsq:

>>> check_solution([0, 0])
x = [ -4.35  -8.69 -19.69   2.31  17.5 ]
E . x - F = [ 0. -0. -0.  0.]
|x| = 28.174593028253167

We can generate and test many other solutions

>>> check_solution(100 * np.random.randn(2))
x = [ -88.93 -139.06   66.64   88.64   17.5 ]
E . x - F = [ 0.  0. -0.  0.]
|x| = 199.62363490542995
>>> check_solution(100 * np.random.randn(2))
x = [-25.2  -26.99  -5.33  16.67  17.5 ]
E . x - F = [ 0. -0. -0.  0.]
|x| = 44.455362582961335
>>> check_solution(100 * np.random.randn(2))
x = [ 93.34  14.57 -55.74 -33.74  17.5 ]
E . x - F = [ 0. -0. -0. -0.]
|x| = 116.09338153741933

We can even look at your solution:

>>> my_favourite_solution = np.array([-12.5, 0, -22, 0, 17.5 ])
>>> my_favourite_coeffs = my_favourite_solution @ E_null
>>> check_solution(my_favourite_coeffs)
x = [-12.5   0.  -22.   -0.   17.5]
E . x - F = [ 0. -0. -0.  0.]
|x| = 30.765240125830324

Upvotes: 3

Related Questions