Using the pseudoinverse to find the linear combination of functions?

Question

I'm working out this problem with a friend of mine on a group project and we are both stuck

Our professor insists that we do all of our work in Maple. I like Maple, but it's not as great as Mathematica.

a) A stock on the Nasdaq stockmarket has the following 9 values at the end of 9 consecutive business days: {(1, $18.35), (2, $25.25), (3, $22.85), (4, $23.65), (5, $22.05), (6, $25.75), (7, $26.15), (8, $26.25), (9, $25.85)}. Use the pseudoinverse method to find the linear combination of the functions

{1, x-5, (x-5)^2, (x-5)^3, (x-5)^4, (x-5)^5, (x-5)^6}

which best fits this data.

I have used the pseudo inverse method and found the linear equation y = .695x + 20.5416666666667

But when plotted with the data, it does not fit the graph very well.

He hinted that we should solve for the letters a through g of this equation

1*a + b*(x-5) + c*(x-5)^2 + d*(x-5)^3 + e*(x-5)^4 + f*(x-5)^5 + g*(x-5)^6

But I am at a loss.

Answer 1

I suspect that the goal is something like the following,

restart:

X := Vector([1,2,3,4,5,6,7,8,9]):
Y := Vector([18.35,25.25,22.85,23.65,22.05,25.75,26.15,26.25,25.85]):

r := Vector[row]([seq((x-5)^i,i=0..6)]):
u := Vector([a,b,c,d,e,f,g]):
p := r . u;

                               2            3            4            5  
   p := a + (x - 5) b + (x - 5)  c + (x - 5)  d + (x - 5)  e + (x - 5)  f

               6  
      + (x - 5)  g

M := Matrix( 9, 7, (i,j) -> eval(r[j],x=X[i]) ):

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

                          [    23.2395104895108]
                          [                    ]
                          [    1.63841491841497]
                          [                    ]
                          [   0.529261460761524]
                          [                    ]
                     C := [  -0.306191724941740]
                          [                    ]
                          [ -0.0192414529914619]
                          [                    ]
                          [  0.0163782051282063]
                          [                    ]
                          [-0.00113888888889210]

sol := eval( p,[seq( u[i]=C[i], i=1..7 )] );

                                                                              2
      sol := 15.0474358974359 + 1.63841491841497 x + 0.529261460761524 (x - 5) 

                                    3                             4
         - 0.306191724941740 (x - 5)  - 0.0192414529914619 (x - 5) 

                                     5                              6
         + 0.0163782051282063 (x - 5)  - 0.00113888888889210 (x - 5) 

# Let's test what the stock solver gets here (using SVD to get
# pseudo-inverse
Statistics:-LinearFit( p, X, Y, x );

                                                                          2
          15.0474358974359 + 1.63841491841492 x + 0.529261460761456 (x - 5) 

                                       3                             4
            - 0.306191724941725 (x - 5)  - 0.0192414529914522 (x - 5) 

                                        5                              6
            + 0.0163782051282051 (x - 5)  - 0.00113888888888892 (x - 5) 

And that candidate solution can be plotted, along with the original data points.

plots:-display( plots:-pointplot( <X|Y> ),
                plot( sol, x=1..9 ) );

It might be made more explicit, by computing the pseudo-inverse of Matrix M using singular value decomposition. Are you supposed to know how to use the singular values and left & right singular vectors to generate a pseudo-inverse (or do a least squares fit)?

[addendum, as per request for clarification]

Take the definitions for r and u from above.

r := Vector[row]([seq((x-5)^i,i=0..6)]);

       [                 2         3         4         5         6]
  r := [1, x - 5, (x - 5) , (x - 5) , (x - 5) , (x - 5) , (x - 5) ]


u := Vector([a,b,c,d,e,f,g]);

                                   [a]
                                   [ ]
                                   [b]
                                   [ ]
                                   [c]
                                   [ ]
                              u := [d]
                                   [ ]
                                   [e]
                                   [ ]
                                   [f]
                                   [ ]
                                   [g]

Observe that r . u produces the polynomial p in x, with unknowns being the names in u. You are going to try and fit this polynomial p to the X and Y data.

Let's just take a single numeric data pair, X[1] and Y[1]. Evaluate r at x=X[1].

eval( r, x=X[1] );

                 [1, -4, 16, -64, 256, -1024, 4096]

The dot product of that with u is the polynomial p evaluated at x=X[1], and we take that as equal to Y[1].

eval( r, x=X[1] ) . u = Y[1];

       a - 4 b + 16 c - 64 d + 256 e - 1024 f + 4096 g = 18.35

This can be done for all 9 data pairs taken from data Vectors X and Y.

That would produce 9 evaluations of p (at each of the X[i] for i=1..9. Each of those equals the corresponding Y[i], i=1..9. So you have 9 equations, with the LHSs being in terms of the 7 unknowns (ie. entries of u, or coefficients in unexpanded p).

Those 9 equations can be written more compactly, in Matrix-Vector notation, as M . u = Y. The ith row of Matrix M consists of the components of r evaluated as X[i].

In other words, the 9 data pairs have produced 9 instances of y=p(x). These are 9 equations that are linear in the 7 unknowns. That is an overdetermined system for which you can attempt to compute a least squares solution.

By multiplying Y by the pseudo-inverse of 9-by-7 Matrix M you obtain a Vector C whose 7 values are taken to be a least squares solution for u.

The least squares solution values C can be substituted in for the unknowns in p. Here are two ways to get that,

eval( p,[seq( u[i]=C[i], i=1..7 )] );

                                                                      2
     15.0474358974359 + 1.63841491841497 x + 0.529261460761524 (x - 5) 

                                  3                             4
        - 0.306191724941740 (x - 5)  - 0.0192414529914619 (x - 5) 

                                    5                              6
        + 0.0163782051282063 (x - 5)  - 0.00113888888889209 (x - 5) 

r.C;

                                                                      2
     15.0474358974359 + 1.63841491841497 x + 0.529261460761524 (x - 5) 

                                   3                             4
        - 0.306191724941740 (x - 5)  - 0.0192414529914619 (x - 5) 

                                    5                              6
        + 0.0163782051282063 (x - 5)  - 0.00113888888889209 (x - 5) 

You might now ask yourself why x-5 is used in the terms of u, and why a 6th degree polynomial was chosen.