Scipy optimize unable to find the correct results

Question

I am trying to use scipy.optimize.minimize to fit parameters for a multivariate function, however, regardless of how many noise free data points I am providing to the optimizer, the optimizer could not converge to a correct (or close) answer.

I wonder if there is a mistake in the way I am using the optimizer but I have been scratching my head to find the mistake. I would appreciate any advice or guesses, thanks!

import numpy as np
from scipy.optimize import minimize
import math

def get_transform(ai,aj,ak,x,y,z):

    i,j,k = 0, 1, 2
    si, sj, sk = math.sin(ai), math.sin(aj), math.sin(ak)
    ci, cj, ck = math.cos(ai), math.cos(aj), math.cos(ak)
    cc, cs = ci*ck, ci*sk
    sc, ss = si*ck, si*sk
    M = np.identity(4)
    M[i, i] = cj*ck
    M[i, j] = sj*sc-cs
    M[i, k] = sj*cc+ss
    M[j, i] = cj*sk
    M[j, j] = sj*ss+cc
    M[j, k] = sj*cs-sc
    M[k, i] = -sj
    M[k, j] = cj*si
    M[k, k] = cj*ci
    M[0, 3] = x
    M[1, 3] = y
    M[2, 3] = z
    
    return M

def camera_intrinsic(fx, ppx, fy, ppy):
    K = np.zeros((3, 3), dtype='float64')
    K[0, 0], K[0, 2] = fx, ppx
    K[1, 1], K[1, 2] = fy, ppy

    K[2, 2] = 1

    return K

def apply_transform(p, matrix):
    rotation = matrix[0:3,0:3]
  
    T = np.array([matrix[0][3],matrix[1][3],matrix[2][3]])
    transformed = (np.dot(rotation, p.T).T)+T
    return transformed

def project(points_3D,internal_calibration):
    points_3D = points_3D.T
    projections_2d = np.zeros((2, points_3D.shape[1]), dtype='float32')
    camera_projection = (internal_calibration).dot(points_3D)
    projections_2d[0, :] = camera_projection[0, :]/camera_projection[2, :]
    projections_2d[1, :] = camera_projection[1, :]/camera_projection[2, :]

    return projections_2d.T

    

def error(x):
    global points,pixels
    transform = get_transform(x[0],x[1],x[2],x[3],x[4],x[5])
    points_transfered = apply_transform(points, transform)
    internal_calibration = camera_intrinsic(x[6],x[7],x[8],x[9])
    projected = project(points_transfered,internal_calibration)
    # print(((projected-pixels)**2).mean())
    return ((projected-pixels)**2).mean()


def generate(points, x):

    transform = get_transform(x[0],x[1],x[2],x[3],x[4],x[5])
    points_transfered = apply_transform(points, transform)
    internal_calibration = camera_intrinsic(x[6],x[7],x[8],x[9])
    projected = project(points_transfered,internal_calibration)
    return projected


points = np.random.rand(100,3)
x_initial = np.random.rand(10)
pixels = generate(points,x_initial)
x_guess = np.random.rand(10)
results = minimize(error,x_guess, method='nelder-mead', tol = 1e-15)
x = results.x
print(x_initial)
print(x)

The optimizer reaches the limit of maximal function evaluations. What about the results if you increase this limit via `minimize(...., options = {'maxiter' : 10**9, 'maxfev': 10**9})` ? — joni, Feb 08 '21 at 20:22
Is it required to use only `'nelder-mead'` method? Can you another one? Also `tol` argument seems too low, because it is just above precision limit for float (about 1e-16). Try increasing it to something like 1e-9 and check if helps the optimization to converge — fdermishin, Feb 09 '21 at 10:08
`'Powell'` optimizer performs better and produces results that are somewhat similar to x_initial, but still fails find the correct solution. It seems that there is a problem in computation of the cost function `error(x)` — fdermishin, Feb 09 '21 at 10:49

fdermishin · Accepted Answer · 2021-02-14T12:14:54.950

You are solving least squares problem, but trying to optimize it using a solver that minimizes a scalar function. While it can possibly solve the problem, it does so very inefficiently. It can require much more iterations or can fail to converge at all.

The better way is to use least_squares instead of minimize.

For it to work properly you should modify error function by returning 1D numpy array instead of a scalar:

def error(x):
    ...
    return (projected-pixels).flatten()

Then call least_squares:

results = least_squares(error, x_guess)
x = results.x
print(x_initial)
print(x)
print('error:', np.linalg.norm(error(x)))

Also, error(x) currently returns array of float32, because an array of float32 is created in project. It should be replaced by float64, otherwise minimization fails to converge, because most of gradients become zeros when 32 bit precision is used.

def project(points_3D,internal_calibration):
    ...
    projections_2d = np.zeros((2, points_3D.shape[1]), dtype='float64')

With these modifications the solver converges to the solution most of the times, but can sometimes fail to do so. It happens because you generate the problem randomly, so in some cases the problem may be degenerate or make no physical sense. Such cases should be investigated on their own.

It can also help to use a robust loss, such as 'arctan', instead of linear loss:

results = least_squares(error, x_guess, loss='arctan')

Result:

[0.68589904 0.68782115 0.83299068 0.02360941 0.19367124 0.54715374
 0.37609235 0.62190714 0.98824796 0.88385802]
[0.68589904 0.68782115 0.83299068 0.02360941 0.19367124 0.54715374
 0.37609235 0.62190714 0.98824796 0.88385802]
error: 1.2269443642313758e-12

Interesting comment about the float32. I had no trouble using `least_squares` with all data as `float32`, but scratched my head for a long time when a seemingly minor calculation turned all my gradients turned into zeros. That minor calculation had turned some of the data into `float64`. I could get it to work with `diff_step=1e-7`, but the better solution was to avoid mixed the dtypes. — adr, Dec 03 '21 at 16:35

Scipy optimize unable to find the correct results

1 Answers1