Fit plane to a set of points in 3D: scipy.optimize.minimize vs scipy.linalg.lstsq

by

Least squares (scipy.linalg.lstsq) is guaranteed to converge. In fact, there is a closed form analytical solution (given by (A^T A)^-1 A^Tb (where ^T is matrix transpose and ^-1 is matrix inversion)

The standard optimization problem, however, is not generally solvable – we are not guaranteed to find a minimizing value. However, for the given equation, find some a, b, c such that z = a*x + b*y + c, we have a linear optimization problem (the constraints and objective are linear in the variables that we’re trying to optimize for). Linear optimization problems are generally solvable, so scipy.optimize.minimize should converge to the optimal value.

Note: this is linear in our constraints even if we do z = a*x + b*y + d*x^2 + e*y^2 + f — we don’t have to restrict ourselves to a linear space of (x,y), since we will have these points (x, y, x^2, y^2) already. To our algorithm, these just look like points in the matrix A. So we can actually get a higher order polynomial using least squares!

A brief aside: In reality, all of the solvers which don’t use an exact analytical solution generally stop within some acceptable range of the actual answer, so it is rarely the case that we get the exact solution, but it tends to be so close that we accept it as exact in practice. Additionally, even least-squares solvers rarely use the analytical solution and instead resort to something quicker like Newton’s Method.

If you were to change the optimization problem, this would not be true. There are certain classes of problems for which we can generally find an optimal value (the largest class of these are called convex optimization problems — although there are many non-convex problems which we can find an optimal value for under certain conditions).

If you’re interested in learning more, take a look at Convex Optimization by Boyd and Vandenberghe. The first chapter doesn’t require much mathematical background and it overviews the general optimization problem and how it relates to solvable optimization problems like linear and convex programming.

More Related Contents:

  1. How to display progress of scipy.optimize function?
  2. scipy.optimize.leastsq with bound constraints
  3. Histogram Matplotlib
  4. Does Python SciPy need BLAS?
  5. Plotting a fast Fourier transform in Python
  6. How to check BLAS/LAPACK linkage in NumPy and SciPy?
  7. How to make scipy.interpolate give an extrapolated result beyond the input range?
  8. Inverse Distance Weighted (IDW) Interpolation with Python
  9. How to create a density plot in matplotlib?
  10. How to do n-D distance and nearest neighbor calculations on numpy arrays
  11. Getting standard errors on fitted parameters using the optimize.leastsq method in python
  12. How to display a 3D plot of a 3D array isosurface in matplotlib mplot3D or similar?
  13. How do I install SciPy on 64 bit Windows?
  14. Multivariate spline interpolation in python/scipy?
  15. Removing duplicate columns and rows from a NumPy 2D array
  16. scipy curve_fit doesn’t like math module
  17. Numpy Resize/Rescale Image
  18. scipy.io.loadmat nested structures (i.e. dictionaries)
  19. Storing numpy sparse matrix in HDF5 (PyTables)
  20. How to calculate the statistics “t-test” with numpy
  21. How do I use numpy.where()? What should I pass, and what does the result mean? [closed]
  22. Creating a Confidence Ellipsis in a scatterplot using matplotlib
  23. Load just part of an image in python
  24. Getting spline equation from UnivariateSpline object
  25. Numpy: Divide each row by a vector element
  26. multidimensional confidence intervals [closed]
  27. Consecutive, Overlapping Subsets of Array (NumPy, Python)
  28. How to display a 3D plot of a 3D array isosurface with mplot3D or similar
  29. Computing cross-correlation function?
  30. Understanding scipy deconvolve

Leave a Comment