scipy - Bivariate Quadratic Polynomial Solution

Question

I am trying to implement the code for the algorithm from this paper: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.123.2361&rep=rep1&type=pdf. I came to the last equation Eq. (5) where I would like to get coefficients a and b. I tried to solve it using scipy in code which you can find below. But I don't understand how these methods work. I just get errors.

Code:

import numpy as np
from scipy import optimize

def solve_for_a_b(coors, error_vec):
     
    cx = coors[:,0,0]
    cy = coors[:,0,1]

    def cost_fun(coeffs):
        a     = coeffs[0]
        b     = coeffs[1]
        c     = coeffs[2]
        c_hat = coeffs[3]

        cost =  np.linalg.norm( a**2*cx**2 + b**2*cy**2 + 2*(a*b)*(cx*cy)
                              + 2*(a*c)*cx + 2*(b*c)*cy + c_hat - re)
    return 4*[cost]

    coeffs = [0,0,0,0]
    sol = optimize.root(cost_fun, coeffs, method = 'lm')

return sol


coors = np.array(
  [[[2269., 1785.]],[[1977., 1787.]],[[1686., 1789.]],[[1393., 1792.]],
   [[1100., 1794.]],[[ 807., 1796.]],[[ 514., 1798.]],[[ 220., 1799.]],
   [[2267., 1492.]],[[1976., 1495.]],[[1684., 1497.]],[[1392., 1499.]],
   [[1099., 1501.]],[[ 806., 1502.]],[[ 512., 1504.]],[[ 219., 1506.]],
   [[2265., 1200.]],[[1974., 1202.]],[[1682., 1204.]],[[1390., 1206.]],
   [[1097., 1207.]],[[ 804., 1209.]],[[ 510., 1211.]],[[ 217., 1212.]],
   [[2263.,  909.]],[[1972.,  910.]],[[1680.,  911.]],[[1388.,  913.]],
   [[1095.,  915.]],[[ 802.,  916.]],[[ 509.,  918.]],[[ 215.,  919.]],
   [[2261.,  617.]],[[1970.,  618.]],[[1678.,  619.]],[[1386.,  620.]],
   [[1093.,  622.]],[[ 800.,  623.]],[[ 507.,  625.]],[[ 214.,  626.]],
   [[2258.,  325.]],[[1967.,  326.]],[[1676.,  327.]],[[1384.,  328.]],
   [[1092.,  329.]],[[ 799.,  330.]],[[ 506.,  332.]],[[ 213.,  333.]]])

error_vec = np.array(
  [-4.76493709e-07, -4.60876991e-07, -4.00291734e-07, -3.44813903e-07,
   -3.42365319e-07, -2.69245444e-07, -1.92865310e-07,  0.00000000e+00,
   -5.18623766e-07, -4.41326967e-07, -4.25394901e-07, -3.65253429e-07,
   -3.11198505e-07, -2.85033873e-07, -2.50725220e-07,  2.08320684e-08,
   -5.09085021e-07, -4.76810581e-07, -3.76267916e-07, -3.41724156e-07,
   -2.78167033e-07, -2.48818089e-07, -1.61401468e-07,  6.93690312e-08,
   -5.68638565e-07, -4.55043519e-07, -4.27832611e-07, -3.47848034e-07,
   -3.13846880e-07, -2.09616720e-07, -1.26412286e-07,  9.22803981e-08,
   -5.82171462e-07, -5.22886081e-07, -4.92479402e-07, -4.32757190e-07,
   -4.23322282e-07, -3.33543040e-07, -2.89759051e-07, -7.35946780e-08,
   -7.25028074e-07, -5.58437533e-07, -5.35705608e-07, -5.02436595e-07,
   -4.56031360e-07, -3.80862038e-07, -3.11127887e-07, -2.20412620e-07])

Answer 1

I solved it using the following:

def find_a_b_coeffs(coors, error_vec):

    def cost_fun(coeffs):
        
        a     = coeffs[0]
        b     = coeffs[1]
        c     = coeffs[2]
        c_hat = coeffs[3]
        
        cx = coors[:,0,0]
        cy = coors[:,0,1]
        ev = error_vec
        
        cost = np.linalg.norm( a**2*cx**2 + b**2*cy**2 + 2*a*b*cx*cy 
                             + 2*c*a*cx + 2*c*b*cy - c_hat - ev)
        
        return cost
    
    coeffs = [0,0,0,0]
    sol = optimize.minimize(cost_fun, coeffs , method='Nelder-Mead', tol=1e-10)
    
    return sol


def plot_sol_3d(coors, e_vec, coeffs ):
  
    cx = coors[:,0,0]
    cy = coors[:,0,1]
    a,b,c,c_hat = coeffs[:]  

    linx = np.linspace(0,np.max(cx),500)
    liny = np.linspace(0,np.max(cy),500)
    xran, yran = np.meshgrid(linx,liny)
    linx = xran.ravel()
    liny = yran.ravel()
    
    z = (a**2*xran**2 + b**2*yran**2 + 2*a*b*xran*yran 
                      + 2*c*a*xran + 2*c*b*yran -c_hat)
    
    fig = plt.figure(dpi = 120)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(cx,cy, error_vec,'.')
    ax.plot_surface(xran, yran,z)
    
    fig = plt.figure(dpi=120)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(xran,yran, a*xran+b*yran)