Use the pseudoinverse to find the conic section of best fit to the data

Question

I am working on a group project and none of us can figure out how to find the answer. Our professor insists that all of our work be done in maple.

The problem is: Use the pseudoinverse to find the conic section of best fit to the data {(-7, -11), (-4, -3), (-1, 5), (2, 8), (6, 1), (9, -4), (13, -10)} and plot it with the data.

He wants this done using both the pseudoinverse method and then part b says:

Now verify that the conic section of part a agrees with the Gauss' regression conic section of best fit to this data.

Any help is appreciated guys!!!

Answer 1

Part a) is somewhat similar to your previous linear least squares problem.

restart:

X := Vector([-7, -4, -1, 2, 6, 9, 13]):
Y := Vector([-11, -3, 5, 8, 1, -4, -10]):

r := Vector[row]([x^2,x*y,y^2,x,y]):
u := Vector([a,b,c,d,e]):

We will use the 7 given data pairs, and p at each pair.

p := r . u = 1;

                          2              2                
                  p := a x  + b x y + c y  + d x + e y = 1

Now, set up the 7 equations in 5 unknowns as M . u = B

M := Matrix(7,5, (i,j) -> eval(r[j],[x=X[i],y=Y[i]]) ):
B := Vector(7,1):

Solve the least squares problem.

C := LinearAlgebra:-LeastSquares( evalf(M), B );

                          [ 0.0535681724900081]
                          [                   ]
                          [ 0.0144560334783153]
                          [                   ]
                     C := [-0.0134458111875932]
                          [                   ]
                          [ -0.195957112862695]
                          [                   ]
                          [  0.226160844326909]

One alternate way,

C := LinearAlgebra:-MatrixInverse( evalf(M) ) . B;

                          [ 0.0535681724900080]
                          [                   ]
                          [ 0.0144560334783152]
                          [                   ]
                     C := [-0.0134458111875931]
                          [                   ]
                          [ -0.195957112862694]
                          [                   ]
                          [  0.226160844326908]

Construct the final answer.

sol := r . C - 1;

                               2                                                2
    sol := 0.0535681724900080 x  + 0.0144560334783152 x y - 0.0134458111875931 y 

       - 0.195957112862694 x + 0.226160844326908 y - 1


with(plots):
display( implicitplot( sol, x=3*min(X)..3*max(X),
                            y=3*min(Y)..3*max(Y), gridrefine=2 ),
         pointplot( <X|Y>, color=red,symbolsize=20 ) );

Another way, using stock command,

Statistics:-LinearFit( r . u, <X|Y>, B, [x,y] );

                       2                                                2
   0.0535681724900080 x  + 0.0144560334783152 x y - 0.0134458111875932 y 

      - 0.195957112862694 x + 0.226160844326909 y

[addendum: following up on the attempt in your comment]

restart:

p5Data:=[[-7, -11], [-4, -3], [-1, 5],
         [2, 8], [6, 1], [9, -4], [13, -10]]:

L := nops(p5Data):
x := [seq(p5Data[j][2], j=1..L)]:
y := [seq(p5Data[j][2], j=1..L)]:

G := add((y[j]-A*x[j]^2-B*x[j]*y[j]-C*y[j]^2-d*x[j]-F)^2, j=1..L):
dGdA := diff(G, A):
dGdB := diff(G, B):
dGdC := diff(G, C):
dGdd := diff(G, d):
dGdF := diff(G, F):

And now, solving,

Soln := solve({dGdA=0,dGdB=0,dGdC=0,dGdF=0,dGdd=0}, {A,B,C,F,d}):

eqnConic := Y = A*X^2 + B*X*Y + C*Y^2 + d*X + F;

                  2              2          
           Y = A X  + B X Y + C Y  + X d + F

# ...which is equivalent to,
eqnConicAlt := ( 0 = rhs(eqnConic) - lhs(eqnConic) )/F;

                2              2              
             A X  + B X Y + C Y  + X d + F - Y
         0 = ---------------------------------
                             F                

eqnConicfloat := evalf( subs(Soln, eqnConicAlt) );

                          2                                      2
     0. = -0.05373994787 X  - 0.01447136784 X Y + 0.01351475167 Y 

    + 0.1969003916 X - 0.2269982688 Y + 1.

Now, plotting,

with(plots):
display( implicitplot( eqnConicfloat, X=3*min(x)..3*max(x),
                                      Y=3*min(y)..3*max(y), gridrefine=2 ),
         pointplot( <<x>|<y>>, color=red,symbolsize=20 ) );