Numpy's 'linalg.solve' and 'linalg.lstsq' not giving same answer as Matlab's '\' or mldivide
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.
%% 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)
# 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 for i in camera_vectors]) cam_y = np.array([i 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) 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 for i in test_cam_vec]) test_cam_y = np.array([i 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 for i in results] results_y = [i for i in results] actual_x = [i for i in test_scr_vec] actual_y = [i 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
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) 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 : x1 Out: 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 : x2 Out: 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 : A_star.dot(x1) Out: 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 : A_star.dot(x2) Out: 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 : A_star.dot(x1-x2) Out: 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 : cond(A_star) Out: 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 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.
The matrix A_matrix is a 6 by 5 matrix, so A_star is a singular matrix. As a result there is no unique solution, and the result of both programs are correct. This corresponds to the original problem being under-determined, as opposed to over-determined.