Numpy's 'linalg.solve' and 'linalg.lstsq' not giving same answer as Matlab's '\' or mldivide

Question

I'm trying to implement the least squares curve fitting algorithm on Python, having already written it on Matlab. However, I'm having trouble getting the right transform matrix, and the problem seems to be happening at the solve step. (Edit: My transform matrix is incredibly accurate with Matlab, but completely off with Python.)

I've looked at numerous sources online, and they all indicate that to translate Matlab's 'mldivide', you have to use 'np.linalg.solve' if the matrix is square and nonsingular, and 'np.linalg.lstsq' otherwise. Yet my results are not matching up.

What is the problem? If it has to do with the implementation of the functions, what is the correct translation of mldivide in numpy?

I have attached both versions of the code below. They are essentially the exact same implementation, with exception to the solve part.

Matlab code:

%% Least Squares Fit

clear, clc, close all

% Calibration Data
scr_calib_pts = [0,0; 900,900; -900,900; 900,-900; -900,-900];
cam_calib_pts = [-1,-1; 44,44; -46,44; 44,-46; -46,-46];
cam_x = cam_calib_pts(:,1);
cam_y = cam_calib_pts(:,2);

% Least Squares Fitting
A_matrix = [];
for i = 1:length(cam_x)
    A_matrix = [A_matrix;1, cam_x(i), cam_y(i), ...
        cam_x(i)*cam_y(i), cam_x(i)^2, cam_y(i)^2];
end
A_star = A_matrix'*A_matrix
B_star = A_matrix'*scr_calib_pts
transform_matrix = mldivide(A_star,B_star)

% Testing Data
test_scr_vec = [200,400; 1600,400; -1520,1740; 1300,-1800; -20,-1600];
test_cam_vec = [10,20; 80,20; -76,87; 65,-90; -1,-80];
test_cam_x = test_cam_vec(:,1);
test_cam_y = test_cam_vec(:,2);

% Coefficients for Transform
coefficients = [];
for i = 1:length(test_cam_x)
    coefficients = [coefficients;1, test_cam_x(i), test_cam_y(i), ...
        test_cam_x(i)*test_cam_y(i), test_cam_x(i)^2, test_cam_y(i)^2];
end

% Mapped Points
results = coefficients*transform_matrix;

% Plotting
test_scr_x = test_scr_vec(:,1)';
test_scr_y = test_scr_vec(:,2)';
results_x = results(:,1)';
results_y = results(:,2)';

figure
hold on
load seamount
s = 50;
scatter(test_scr_x, test_scr_y, s, 'r')
scatter(results_x, results_y, s)

Python code:

# Least Squares fit

import numpy as np
import matplotlib.pyplot as plt

# Calibration data
camera_vectors = np.array([[-1,-1], [44,44], [-46,44], [44,-46], [-46,-46]])
screen_vectors = np.array([[0,0], [900,900], [-900,900], [900,-900], [-900,-900]])

# Separate axes
cam_x = np.array([i[0] for i in camera_vectors])
cam_y = np.array([i[1] for i in camera_vectors])

# Initiate least squares implementation
A_matrix = []
for i in range(len(cam_x)):
    new_row = [1, cam_x[i], cam_y[i], \
        cam_x[i]*cam_y[i], cam_x[i]**2, cam_y[i]**2]
    A_matrix.append(new_row)
A_matrix = np.array(A_matrix)
A_star = np.transpose(A_matrix).dot(A_matrix)

B_star = np.transpose(A_matrix).dot(screen_vectors)

print A_star
print B_star

try:
    # Solve version (Implemented)
    transform_matrix = np.linalg.solve(A_star,B_star)
    print "Solve version"
    print transform_matrix
except:
    # Least squares version (implemented)
    transform_matrix = np.linalg.lstsq(A_star,B_star)[0]
    print "Least Squares Version"
    print transform_matrix


# Test data
test_cam_vec = np.array([[10,20], [80,20], [-76,87], [65,-90], [-1,-80]])
test_scr_vec = np.array([[200,400], [1600,400], [-1520,1740], [1300,-1800], [-20,-1600]])

# Coefficients of quadratic equation
test_cam_x = np.array([i[0] for i in test_cam_vec])
test_cam_y = np.array([i[1] for i in test_cam_vec])    
coefficients = []
for i in range(len(test_cam_x)):
    new_row = [1, test_cam_x[i], test_cam_y[i], \
        test_cam_x[i]*test_cam_y[i], test_cam_x[i]**2, test_cam_y[i]**2]
    coefficients.append(new_row)
coefficients = np.array(coefficients)

# Transform camera coordinates to screen coordinates
results = coefficients.dot(transform_matrix)

# Plot points
results_x = [i[0] for i in results]
results_y = [i[1] for i in results]
actual_x = [i[0] for i in test_scr_vec]
actual_y = [i[1] for i in test_scr_vec]

plt.plot(results_x, results_y, 'gs', actual_x, actual_y, 'ro')
plt.show()

Edit (in accordance with a suggestion):

# Transform matrix with linalg.solve

[[  2.00000000e+01   2.00000000e+01]
 [ -5.32857143e+01   7.31428571e+01]
 [  7.32857143e+01  -5.31428571e+01]
 [ -1.15404203e-17   9.76497106e-18]
 [ -3.66428571e+01   3.65714286e+01]
 [  3.66428571e+01  -3.65714286e+01]]

# Transform matrix with linalg.lstsq:

[[  2.00000000e+01   2.00000000e+01]
 [  1.20000000e+01   8.00000000e+00]
 [  8.00000000e+00   1.20000000e+01]
 [  1.79196935e-15   2.33146835e-15]
 [ -4.00000000e+00   4.00000000e+00]
 [  4.00000000e+00  -4.00000000e+00]]

% Transform matrix with mldivide:

   20.0000   20.0000
   19.9998    0.0002
    0.0002   19.9998
         0         0
   -0.0001    0.0001
    0.0001   -0.0001

Answer 1

The interesting thing is that you will get quite different results with np.linalg.lstsq and np.linalg.solve.

x1 = np.linalg.lstsq(A_star, B_star)[0]
x2 = np.linalg.solve(A_star, B_star)

Both should offer a solution for the equation Ax = B. However, these give two quite different arrays:

In [37]: x1    
Out[37]: 
array([[  2.00000000e+01,   2.00000000e+01],
       [  1.20000000e+01,   7.99999999e+00],
       [  8.00000001e+00,   1.20000000e+01],
       [ -1.15359111e-15,   7.94503352e-16],
       [ -4.00000001e+00,   3.99999999e+00],
       [  4.00000001e+00,  -3.99999999e+00]]


In [39]: x2
Out[39]: 
array([[  2.00000000e+01,   2.00000000e+01],
       [ -4.42857143e+00,   2.43809524e+01],
       [  2.44285714e+01,  -4.38095238e+00],
       [ -2.88620104e-18,   1.33158696e-18],
       [ -1.22142857e+01,   1.21904762e+01],
       [  1.22142857e+01,  -1.21904762e+01]])

Both should give an accurate (down to the calculation precision) solution to the group of linear equations, and for a non-singular matrix there is exactly one solution.

Something must be then wrong. Let us see if this both candidates could be solutions to the original equation:

In [41]: A_star.dot(x1)
Out[41]: 
array([[ -1.11249392e-08,   9.86256055e-09],
       [  1.62000000e+05,  -1.65891834e-09],
       [  0.00000000e+00,   1.62000000e+05],
       [ -1.62000000e+05,  -1.62000000e+05],
       [ -3.24000000e+05,   4.47034836e-08],
       [  5.21540642e-08,  -3.24000000e+05]])

In [42]: A_star.dot(x2)
Out[42]: 
array([[ -1.45519152e-11,   1.45519152e-11],
       [  1.62000000e+05,  -1.45519152e-11],
       [  0.00000000e+00,   1.62000000e+05],
       [ -1.62000000e+05,  -1.62000000e+05],
       [ -3.24000000e+05,   0.00000000e+00],
       [  2.98023224e-08,  -3.24000000e+05]])

They seem to give the same solution, which is essentially the same as B_star as it should be. This leads us towards the explanation. With simple linear algebra we could predict that A . (x1-x2) should be very close to zero:

In [43]: A_star.dot(x1-x2)
Out[43]: 
array([[ -1.11176632e-08,   9.85164661e-09],
       [ -1.06228981e-09,  -1.60071068e-09],
       [  0.00000000e+00,  -2.03726813e-10],
       [ -6.72298484e-09,   4.94765118e-09],
       [  5.96046448e-08,   5.96046448e-08],
       [  2.98023224e-08,   5.96046448e-08]])

And it indeed is. So, it seems that there is a non-trivial solution for the equation Ax = 0, the solution being x = x1-x2, which means that the matrix is singular and there are thus an infinite number of different solutions for Ax=B.

The problem is thus not in NumPy or Matlab, it is in the matrix itself.

---

However, in the case of this matrix the situation is a bit tricky. A_star seems to be singular by the definition above (Ax=0 for x<>0). On the other hand its determinant is non-zero, and it is not singular.

In this case A_star is an example of a matrix which is numerically unstable while not singular. The solve method solves it by using the simple multiplication-by-inverse. This is a bad choice in this case, as the inverse contains very large and very small values. This makes the multiplicaion prone to round-off errors. This can be seen by looking at the condition number of the matrix:

In [65]: cond(A_star)
Out[65]: 1.3817810855559592e+17

This is a very high condition number, and the matrix is ill-conditioned.

In this case the use of an inverse to solve the problem is a bad approach. The least squares approach gives better results, as you may see. However, the better solution is to rescale the input values so that x and x^2 are in the same range. One very good scaling is to scale everything between -1 and 1.

---

One thing you might consider is to try to use NumPy's indexing capabilities. For example:

cam_x = np.array([i[0] for i in camera_vectors])

is equivalent to:

cam_x = camera_vectors[:,0]

and you may construct your array A this way:

A_matrix = np.column_stack((np.ones_like(cam_x), cam_x, cam_y, cam_x*cam_y, cam_x**2, cam_y**2))

No need to create lists of lists or any loops.