Numpy 'linalg.solve' and 'linalg.lstsq' do not give the same answer as Matlab '' or mldivide

Question

I am trying to implement a least squares curve approximation algorithm in Python, already writing it in Matlab. However, it’s hard for me to get the right transformation matrix, and the problem seems to be happening at the solution stage. (Edit: My transformation matrix is ​​incredibly accurate with Matlab, but completely off with Python.)

I have looked at many sources on the Internet, and all of them indicate that to translate Matlab "mldivide" you need to use "np.linalg.solve" if the matrix is ​​square and non-singular, and "np.linalg.lstsq 'otherwise. But mine Results do not match.

What is the problem? If this is related to the implementation of functions, what is the correct translation of mldivide to numpy?

I have attached both versions of the code below. They are essentially the same implementation, with the exception of the crucial 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() 

Change (as proposed):

 # 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

Interestingly, you will get completely 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, they give two completely 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 exact (up to the accuracy of calculation) solution of a group of linear equations, and for a nonsingular matrix, exactly one solution.

Something must be wrong. Let's see if these two candidates can 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 essentially coincides with B_star as it should be. This leads us to an 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 indeed it is. Thus, it seems that for the equation Ax = 0 there is a nontrivial solution, x = x1-x2, which means that the matrix is ​​singular and, therefore, there are an infinite number of different solutions for Ax = B.

Thus, the problem is not in NumPy or Matlab, but in the matrix itself.

---

However, in the case of this matrix, the situation is a bit complicated. A_star seems to be singular by definition above (Ax = 0 for x & lt> 0). On the other hand, its determinant is nonzero and not singular.

In this case, A_star is an example of a matrix that is numerically unstable rather than singular. The solve method solves it using simple multiplication by the inverse. In this case, this is a bad choice, since the inversion contains very large and very small values. This makes the animation prone to rounding 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, and the matrix is ​​poorly conditioned.

In this case, using a reverse solution to the problem is a bad approach. As you can see, the least squares approach gives the best results. However, the best 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 can consider is to try to use the NumPy indexing capabilities. For example:

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

equivalent to:

 cam_x = camera_vectors[:,0] 

and you can build your array A as follows:

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

There is no need to create lists of lists or any loops.

Answer 2

Matrix A_matrix is ​​a 6 × 5 matrix, so A_star is a special matrix. As a result, there is no single solution, and the result of both programs is correct. This corresponds to the outstanding original task, and not overridden.