CODEWITHSUNDEEP
pythonnumpy
abeagomez
abeagomez

Reputation: 602

Using numpy and lstsq to solve a 3 dimensions system

I´m trying to transform x,y,z real world coordinates to my own x,y,z virtual world coordinates. As there is noise while getting the real world coordinates I need to use least squares. I have 3 variables as input: r_x,r_y,r_z and I need to have a 3 variables output v_x, v_y, v_z. I know that what I need is just a matrix A so [r_x, r_y, r_z]A =[v_x, v_y, v_z].

Here is my code:

>>> import numpy as np
>>> x = np.array([[1,2,3],[4,5,6],[7,8,9]])
>>> y = np.array([[1,2,3],[4,5,6],[7,8,9]])
>>> A = np.vstack([x, np.ones(len(x))]).T
>>> c, resid, rank,  sigma =np.linalg.lstsq(A,y)
>>> c,m = c[0:3], c[-1]
>>> x1 = np.array([1,2,3])
>>> np.dot(x1,c) + m
array([ 1.94736842,  2.31578947,  2.68421053])

As it can be seen, the answer is wrong, as the initial system has a unique solution and the least squares result is a really bad approximation. Does anybody knows what am I doing wrong?

Upvotes: 1

Views: 1453

Answers (1)

Stephen Rauch
Stephen Rauch

Reputation: 49812

The primary problem is that the rank of your input data was not sufficient to allow accurate inversion of your A matrix. Consider the following:

Test Code:

import numpy as np

def build_a(x_data):
    return np.column_stack((x_data, np.ones(len(x_data))))

def lstsq(x_data, y_data):
    return np.linalg.lstsq(build_a(x_data), y_data)

def show_rank(x_data):
    print(lstsq(x_data, x_data)[2])

show_rank([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
show_rank([[1, 2, 3], [4, 6, 6], [7, 8, 9]])
show_rank([[1, 2, 3], [4, 6, 6], [7, 7, 9]])
show_rank([[1, 2, 3], [4, 6, 6], [7, 7, 9], [7, 8, 10]])
show_rank(np.random.rand(10, 3))

Matrix Ranks:

2
3
3
4
4

This issue is that [1, 2, 3] and [4, 5, 6] are on the same vector, and do not provide anymore information, so by mixing up the pairs slightly the rank can be increased. But we need rank 4 since we are also getting the offsets, so we need at least 4 samples to get a rank 4 matrix.

Example

x = np.array([[1, 2, 3], [4, 5, 16], [17, 8, 9], [7, 8, 19]])
C, resid, rank, sigma = lstsq(x, x + 1)

c, m = C[0:3], C[-1]
x1 = np.array([1, 2, 4])
print(np.dot(x1, c) + m)
print(rank)

Results:

[ 2.  3.  5.]
4

Upvotes: 1

Related Questions