Module for PA = LU Factorization with Pivoting

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

Background

Definition (LU-Factorization).  The nonsingular matrix A has an LU-factorization if it can be expressed as the product of a lower-triangular matrix L and an upper triangular matrix U:  

		  A=LU.  

When this is possible we say that A has an LU-decomposition.  It turns out that this factorization (when it exists) is not unique.  If L has 1's on it's diagonal, then it is called a Doolittle factorization.  If U has 1's on its diagonal, then it is called a Crout factorization.  When  L=(U^T), it is called a Cholesky decomposition.  In this module we will develop an algorithm that produces a Doolittle factorization.

Theorem (LU Factorization with NO pivoting).  If row interchanges are not needed to solve the linear system AX = B, then A has the LU factorization (illustrated with 4×4 matrices).  

		(({{(Subscript[a, 1, 1]), (Subscript[a, 1, 2]), (Subscript[a, 1, 3]), (Subscript[a, 1, 
                      4])}, {(Subscript[a, 2, 1]), (Subscript[a, 2, 2]), (Subscript[a, 2, 3]), (Subscript[a, 2, 
                      4])}, {(Subscript[a, 3, 1]), (Subscript[a, 3, 2]), (Subscript[a, 3, 3]), (Subscript[a, 3, 
                      4])}, {(Subscript[a, 4, 1]), (Subscript[a, 4, 2]), (Subscript[a, 4, 3]), (Subscript[a, 4, 
                      4])}})) = (({{1, 0, 0, 0}, {(Subscript[l, 2, 1]), 1, 0, 0}, {(Subscript[l, 3, 1]), (Subscript[l, 3, 2]), 1, 0}, {(Subscript[l, 4, 1]), (Subscript[l, 4, 2]), (Subscript[l, 4, 3]), 1}}))(({{(Subscript[u, 1, 1]), (Subscript[u, 1, 2]), (Subscript[u, 1, 3]), (Subscript[u, 1, 
                    4])}, {0, (Subscript[u, 2, 2]), (Subscript[u, 2, 3]), (Subscript[u, 2, 4])}, {0, 0, (Subscript[u, 3, 3]), (Subscript[u, 3, 4])}, {0, 0, 0, (Subscript[u, 4, 4])}})).
		
Remark 1.  This is not a linear system. 
Remark 2.  The easy solution uses row vectors and is a modification of limited Gauss Jordan elimination. 
Remark 3.  A sufficient condition for the factorization to exist is that all principal minors of A are nonsingular.

Proof  LU Factorization  LU Factorization

Mathematica Subroutine (Limited Gauss-Jordan Elimination).

GaussJordan[A0_, 
        n_]:=
Module[ {A = A0, i, 
          p},  
Print[ 
          MatrixForm[A] ];  
For[ p = 1,  
          p ≤ 
            n,  p++,  
(Subscript[A, [[p]]])  =  ((Subscript[A, [[p]]])/(Subscript[A, [[p, p]]]));  
For[  i = 1,  
            i ≤ 
              n,  i++,  
If[  i ≠ 
                  p,  
(Subscript[A, [[i]]])  =  
                  (Subscript[A, [[i]]])  -  (Subscript[A, [[i, p]]]) 
                      (Subscript[A, [[p]]])  ]; ];  
Print[ 
            MatrixForm[A] ];  ];  ]

Computer Programs  LU Factorization  LU Factorization

Mathematica Subroutine (LandU).

\(\(\(LandU[A0_, 
        n_]\)\(:=\)\(
\)\(Module[\ {A = A0, i, 
          p}, \ 
U = A0; \ 
L = 
          IdentityMatrix[n]; \ 
Print[MatrixForm[L], 
          MatrixForm[U], " = ", 
          MatrixForm[A0]]; \ 
For[\ p = 1, \ 
          p ≤ 
            n - 1, \ \(p++\), \ 
\(For[\ i = p + 1, \ 
              i ≤ n, \ \(i++\), \ 
m = 
                A\_\(\([[\)\(i, p\)\(]]\)\)\/A\_\(\([[\)\(p, p\)\(]]\)\); 
L\_\(\([[\)\(i, p\)\(]]\)\) = 
                m; \ 
A\_\(\([[\)\(i\)\(]]\)\)\  = \ 
                A\_\(\([[\)\(i\)\(]]\)\)\  - \ m\ A\_\(\([[\)\(p\)\(]]\)\); 
U\  = \ A; \ 
Print[MatrixForm[L], 
                MatrixForm[U], " = ", 
                MatrixForm[A0]];\ \ ];\)\ \ ];\ \ ]\)\(\ \)\)\)

Example 1. Given  A=(({{4, 2, 3}, {-3, 1, 4}, {2, 4, 5}})).  Find matrices L and U so that LU = A.

Solution 1.

A = (({{4, 2, 3}, {-3, 1, 4}, {2, 4, 5}}));

Use the LandU subroutine and construct matrices L and U so that LU = A.

LandU[A, 3];

Does  L.U = A ?

Print["A = ", MatrixForm[A]];
Print["L = ", MatrixForm[L]];
Print["U = ", MatrixForm[U]];
Print["Does LU = A ?"]
Print[MatrixForm[L], MatrixForm[U], " = ", 
        MatrixForm[A]];
Print[L . U == A];

Example 2. Given  A=(({{4, 2, -1, 3}, {3, -4, 2, 5}, {-2, 6, -5, -2}, {5, 1, 6, -3}})).  Find matrices L and U so that LU = A.

Solution 2.

A = (({{4, 2, -1, 3}, {3, -4, 2, 5}, {-2, 6, -5, -2}, {5, 1, 6, -3}}));

Use the LandU subroutine and construct matrices L and U so that LU = A.

LandU[A, 4];

Does  L.U = A ?

Print["A = ", MatrixForm[A]];
Print["L = ", MatrixForm[L]];
Print["U = ", MatrixForm[U]];
Print["Does LU = A ?"]
Print[MatrixForm[L], MatrixForm[U], " = ", 
        MatrixForm[A]];
Print[L . U == A];

Example 3. Given  A=(({{1, 2, 6}, {4, 8, -1}, {-2, 3, 5}})).  Can A be factored  A = LU?

Solution 3.

A = (({{1, 2, 6}, {4, 8, -1}, {-2, 3, 5}}));

Use the LandU subroutine and see what happens.

LandU[A, 3];

We get a division by zero error, hence there is no LU-decomposition of A. 

Remark.  The matrix A is nonsingular.

We are done.  This shows that nonsingular is not sufficient to guarantee an  LU-decomposition. 

Aside.  We can interchange rows and the new matrix B will have an   LU-decomposition.

A = (({{1, 2, 6}, {4, 8, -1}, {-2, 3, 5}})); 
B = 
        (Subscript[A, [[{1, 3, 2}]]]);
A = 
        B;
Print["                  B = ", 
        MatrixForm[B]];
LandU[A, 3]; 
Print[MatrixForm[L], 
        MatrixForm[U], " = ", MatrixForm[B]];
Print[L . U == B];

Print["B = ", MatrixForm[B]];
Print["L = ", MatrixForm[L]];
Print["U = ", MatrixForm[U]];
Print["Does LU = B ?"]
Print[MatrixForm[L], MatrixForm[U], " = ", 
        MatrixForm[B]];
Print[L . U == B];

Theorem (A = LU;  Factorization with NO Pivoting).  Assume that A has a Doolittle, Crout or Cholesky factorization.  The solution X to the linear system  AX = B,  is found in three steps:  

	1.  Construct the matrices  L and  U, if possible.  
	2.  Solve  LY = B  for  Y  using forward substitution.
	3.  Solve  UX = Y  for  X  using back substitution.

Proof  LU Factorization  LU Factorization

The above theorem assumes that there are no row interchanges.  We have seen in Example 3 an example of a nonsingular matrix A could not be directly factored as A = LU.  If row interchanges are permitted then a factorization of a "permuted matrix" will be obtained.  A permutation of the first n positive integers  1,2,...,n.  is an arrangement  (Subscript[k, 1]),(Subscript[k, 2]),...,(Subscript[k, n])  of these integers in a definite order.  For example,  1,4,2,3,5  is a permutation of the five integers  1,2,3,4,5.  The standard base vectors  (Subscript[E, i]) = [0,...,(Subscript[1, i]),0,...,0]  for  i = 1,2,...,n  are used in the next definition.

Definition (Permutation Matrix).  An n×n permutation matrix P is a matrix with precisely one entry whose value is "1" in each column and row, and all of whose other entries are "0."  The rows of  P  are a permutation of the rows of the identity matrix and  P  can be written as  

		P =(([(Subscript[E, (Subscript[k, 1])]^T),(Subscript[E, (Subscript[k, 2])]^T),...,(Subscript[E, (Subscript[k, n])]^T)])^T).

The elements of  P =[(Subscript[p, i, j])]  have the form  

		(Subscript[p, i, j]) = {({{1,   for  j =(Subscript[k, i]) }, {0, otherwise}})}

Theorem.  Suppose that  P =(([(Subscript[E, (Subscript[k, 1])]^T),(Subscript[E, (Subscript[k, 2])]^T),...,(Subscript[E, (Subscript[k, n])]^T)])^T)  is a permutation matrix.  The product  PA  is a new matrix whose rows consists of the rows of  A  rearranged in the new order  (Subscript[row, (Subscript[k, 1])])A,(Subscript[row, (Subscript[k, 2])])A,...,(Subscript[row, (Subscript[k, n])])A.

For example, the permutation matrix  P = (({{0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}))  will interchange rows 1 and 2 and also interchange rows 3 and 4.

Clear[a];
MatrixForm[(({{0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}})).(({{(Subscript[a, 1, 1]), (Subscript[a, 1, 2]), (Subscript[a, 1, 3]), (Subscript[a, 1, 
                      4])}, {(Subscript[a, 2, 1]), (Subscript[a, 2, 2]), (Subscript[a, 2, 3]), (Subscript[a, 2, 
                      4])}, {(Subscript[a, 3, 1]), (Subscript[a, 3, 2]), (Subscript[a, 3, 3]), (Subscript[a, 3, 
                      4])}, {(Subscript[a, 4, 1]), (Subscript[a, 4, 2]), (Subscript[a, 4, 3]), (Subscript[a, 4, 
                      4])}})) ]

Theorem.  If  A  is a nonsingular matrix, then there exists a permutation matrix  P  so that  PA  has an LU-factorization 

		PA = LU.

Theorem (PA = LU;  Factorization with Pivoting).  Given that A is nonsingular.  The solution X to the linear system  AX = B, is found in four steps:  

	1.  Construct the matrices  L,U  and  P.  
	2.  Compute the column vector   PB.
	3.  Solve  LY = PB  for  Y  using forward substitution.
	4.  Solve  UX = Y  for  X  using back substitution.

Proof  LU Factorization  LU Factorization

Computer Programs  LU Factorization  LU Factorization

Mathematica Subroutine (LUfactor).
Remark.  The matrix  A  is a global variable and elements are changed when the  LUfactor  subroutine is executed.  Hence, it is important to save a copy of  A  in the variable  A0.
Use one of the following two versions of the  LUfactor  subroutine and the  SolveLU  subroutine.  The first version of  LUfactor  uses parallel programming and "row operations" instead of loops.

\(\(\(LUfactor[
        n_]\)\(:=\)\(
\)\(Module[\ {i, j, k, p, 
          Z}, 
Print["A=", 
          MatrixForm[A]]; 
r\  = \ Table[j, {j, n}]; 
        Print["r=", r]; 
For[\ p = 1, \ 
          p ≤ n - 1, \ \(p++\), 
For[\ 
            j = p + 1, \ 
            j ≤ 
              n, \ \(j++\), 
\(If[\ 
                Abs[A\_\(\([[\)\(r\_\(\([[\)\(j\)\(]]\)\), p\)\(]]\)\)] > 
                  Abs[A\_\(\([[\)\(r\_\(\([[\)\(p\)\(]]\)\), p\)\(]]\)\)], 
r\_\(\([[\)\({j, p}\)\(]]\)\)\  = \ 
                  r\_\(\([[\)\({p, j}\)\(]]\)\)\ \ ];\)\ ]; 
          Print["r=", r]; 
For[\ k = p + 1, \ 
            k ≤ 
              n, \ \(k++\), 
A\_\(\([[\)\(r\_\(\([[\)\(k\)\(]]\)\), p\)\(]]\)\)\  = \ 
              A\_\(\([[\)\(r\_\(\([[\)\(k\)\(]]\)\), p\)\(]]\)\)\/A\_\(\([[\)\(r\_\(\([[\)\(p\)\(]]\)\), p\)\(]]\)\); \ 
A\_\(\([[\)\(\ \)\(r\_\(\([[\)\(k\)\(]]\)\), Range[p + 1, n]\)\(\ \)\(]]\)\)\  = \ 
              A\_\(\([[\)\(\ \)\(r\_\(\([[\)\(k\)\(]]\)\), Range[p + 1, n]\)\(\ \)\(]]\)\)\  - \ 
                A\_\(\([[\)\(r\_\(\([[\)\(k\)\(]]\)\), p\)\(]]\)\)\ A\_\(\([[\)\(\ \)\(r\_\(\([[\)\(p\)\(]]\)\), Range[p + 1, n]\)\(\ \)\(]]\)\);\ \ ];\ ]; 
L\  = \ \(P\  = \ 
            IdentityMatrix[n]\); 
Z\  = \ 
          Table[0, {j, n}]; 
P\  = \ 
          P\_\(\([[\)\(r\)\(]]\)\); 
U\  = \ 
          A\_\(\([[\)\(r\)\(]]\)\); 
For[\ i = 1, \ 
          i ≤ 
            n, \ \(i++\), 
L\_\(\([[\)\(\ \)\(i, Range[1, i - 1]\)\(\ \)\(]]\)\)\  = \ 
            A\_\(\([[\)\(\ \)\(r\_\(\([[\)\(i\)\(]]\)\), Range[1, i - 1]\)\(\ \)\(]]\)\); \ 
U\_\(\([[\)\(\ \)\(i, Range[1, i - 1]\)\(\ \)\(]]\)\)\  = \ 
            Z\_\(\([[\)\(\ \)\(Range[1, i - 1]\)\(\ \)\(]]\)\);\ \ ]; \ Print["L=", MatrixForm[L]]; 
        Print["U=", MatrixForm[U]];]\)\(\ \ \)\)\)

Use the subroutine SolveLU which is similar to the forward substitution and back substitution subroutines.

Mathematica Subroutine (SolveLU).

SolveLU[
          n_] := 
Module[{k}, 
X = Y = 
              Table[0, {n}]; 
For[ k = 1,  
            k ≤ 
              n,  k++, 
(Subscript[Y, [[k]]])  =  ((Subscript[B, [[(Subscript[r, [[k]]])]]]) - (∑_(j =  1)^(k - 1))(Subscript[A, [[(Subscript[r, [[k]]]), j]]]) 
                        (Subscript[Y, [[j]]])); ]; 
For[ k = n,  
            1 ≤ 
              k,  k--, 
(Subscript[X, [[k]]])  =  ((Subscript[Y, [[k]]]) - (∑_(j =  k + 1)^n) 
                            (Subscript[A, [[(Subscript[r, [[k]]]), 
                            j]]]) 
                        (Subscript[X, [[j]]]))/
                (Subscript[A, [[(Subscript[r, [[k]]]), k]]]) ]; ];

Remark.  Everything has been carefully set up so that  L,  U, and  P  can all be studied.

Example 4.  Use  PA = LU  factorization with pivoting to solve the linear system  (({{2.0, 4.0, 1.0}, {4.0, -10, 2.0}, {1.0, 2.0, 4.0}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}})) = (({{5.0}, {-8.0}, {13.0}})).

Solution 4.

Remark.  The  LU  factorization part does not depend on  B.  To emphasize this, we only enter the matrix  A  at this time.

A =(({{2.0, 4.0, 1.0}, {4.0, -10.0, 2.0}, {1.0, 2.0, 4.0}}));
A0  =  A;

Execute the  LUfactor  subroutine.  Note that the dimension of the system must be used, in this case n = 3.

LUfactor[3];

Now print out the matrices  L,  U,  L.U,  L.U and  P.A.  (Remember we stored a copy of  A  in A0).

Print[" L = ", MatrixForm[L] ];
Print[" U = ", MatrixForm[U] ];
Print[" P = ", MatrixForm[P] ];
Print[""];
Print[" A = ", MatrixForm[A0] ];
Print[""];
Print["LU = ", MatrixForm[L . U] ];
Print["PA = ", MatrixForm[P . A0] ];

Now enter the matrix  B.

B = (({{5.0}, {-8.0}, {13.0}}));

Solve the system  AX = B for  X.

SolveLU[3];
\(\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])""};\)\)
X = vars;
Print["Solve AX = B"];
Print[MatrixForm[A0] , MatrixForm[vars], " = ", 
        MatrixForm[A0 . X]];
Print[" X = ", MatrixForm[vars], " = ", 
        MatrixForm[X]  ];
Print[""];
Print["AX = ", MatrixForm[A0] , MatrixForm[X], " = ", 
        MatrixForm[A0 . X]];
Print[" B = ", MatrixForm[B] ];

We are done.

Aside.  For curiosity, we can check to see the "in between" steps that were made in the solution.  First, we solved  LY = PB  for  Y  using forward substitution.  This is just for fun

Print["L = ", MatrixForm[L] ];
Print["Y = ", MatrixForm[Y] ];
Print["P = ", MatrixForm[P] ];
Print["LY = ", MatrixForm[L], MatrixForm[Y], " = ", 
        MatrixForm[L . Y] ];
Print["PB = ", MatrixForm[P], MatrixForm[B], " = ", 
        MatrixForm[P . B] ];

Second, we solved  UX = Y  for  X  using back substitution.

Print[" U = ", MatrixForm[U] ];
Print[" X = ", MatrixForm[X] ];
Print["UX = ", MatrixForm[U] , MatrixForm[X], " = ", 
        MatrixForm[U . X]];
Print[" Y = ", MatrixForm[Y] ];

Finally, we can verify that the residual is  zero, i. e. use matrix multiplication to see if    AX - B  =  0.

Print["AX - B" ];
Print["= ", MatrixForm[A0] , MatrixForm[X], " - ", 
        MatrixForm[B]];
Print["= ", MatrixForm[A0 . X], " - ", 
        MatrixForm[B], " = ", MatrixForm[A0 . X - B]];

Aside.  Suppose a new linear system  AX = B  is to be solved, where  A  is the same we used above, but  B  is a new vector.
Then we enter  B  and only invoke the  SolveLU  subroutine.

B = (({{1.0}, {5.0}, {4.0}}));

Solve the new system  AX = B for  X.

SolveLU[3];
\(\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])""};\)\)
Print["Solve AX = B"];
Print[MatrixForm[A0] , MatrixForm[vars], " = ", 
        MatrixForm[A0 . X]];
Print[" X = ", MatrixForm[vars], " = ", 
        MatrixForm[X]  ];
Print[""];
Print["AX = ", MatrixForm[A0] , MatrixForm[X], " = ", 
        MatrixForm[A0 . X]];
Print[" B = ", MatrixForm[B] ];

Example 5.  Use  PA = LU  factorization with pivoting to solve the linear system  (({{4.0, 2.0, -1.0, 3.0}, {3.0, -4.0, 2.0, 5.0}, {-2.0, 6.0, -5.0, -2.0}, {5.0, 1.0, 6.0, -3.0}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{16.9}, {-14.0}, {25.0}, {9.4}})).

Solution 5.

Remark.  The  LU  factorization part does not depend on  B.  To emphasize this, we only enter the matrix  A  at this time.

A = (({{4.0, 2.0, -1.0, 3.0}, {3.0, -4.0, 2.0, 5.0}, {-2.0, 6.0, -5.0, -2.0}, {5.0, 1.0, 6.0, -3.0}})); 
A0  =  A;

Execute the  LUfactor  subroutine.  Note that the dimension of the system must be used, in this case n = 3.

LUfactor[4];

Now print out the matrices  L,  U,  L.U,  L.U and  P.A.  (Remember we stored a copy of  A  in A0).

Print[" L = ", MatrixForm[L] ];
Print[" U = ", MatrixForm[U] ];
Print[" P = ", MatrixForm[P] ];
Print[""];
Print[" A = ", MatrixForm[A0] ];
Print[""];
Print["LU = ", MatrixForm[L . U] ];
Print["PA = ", MatrixForm[P . A0] ];

Now enter the matrix  B.

B = (({{16.9}, {-14.0}, {25.0}, {9.4}}));

Solve the system  AX = B for  X.

SolveLU[4];
\(\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)\)
Print["Solve AX = B"];
Print[MatrixForm[A0] , MatrixForm[vars], " = ", 
        MatrixForm[A0 . X]];
Print[" X = ", MatrixForm[vars], " = ", 
        MatrixForm[X]  ];
Print[""];
Print["AX = ", MatrixForm[A0] , MatrixForm[X], " = ", 
        MatrixForm[A0 . X]];
Print[" B = ", MatrixForm[B] ];

We are done.

Aside.  For curiosity, we can check to see the "in between" steps that were made in the solution.  First, we solved  LY = PB  for  Y  using forward substitution.  This is just for fun

Print["L = ", MatrixForm[L] ];
Print["Y = ", MatrixForm[Y] ];
Print["P = ", MatrixForm[P] ];
Print["LY = ", MatrixForm[L], MatrixForm[Y], " = ", 
        MatrixForm[L . Y] ];
Print["PB = ", MatrixForm[P], MatrixForm[B], " = ", 
        MatrixForm[P . B] ];

Second, we solved  UX = Y  for  X  using back substitution.

Print[" U = ", MatrixForm[U] ];
Print[" X = ", MatrixForm[X] ];
Print["UX = ", MatrixForm[U] , MatrixForm[X], " = ", 
        MatrixForm[U . X]];
Print[" Y = ", MatrixForm[Y] ];

Finally, we can verify that the residual is  zero, i. e. use matrix multiplication to see if    AX - B  =  0.

Print["AX - B" ];
Print["= ", MatrixForm[A0] , MatrixForm[X], " - ", 
        MatrixForm[B]];
Print["= ", MatrixForm[A0 . X], " - ", 
        MatrixForm[B], " = ", MatrixForm[A0 . X - B]];

Aside.  Suppose a new linear system  AX = B  is to be solved, where  A  is the same we used above, but  B  is a new vector.
Then we enter  B  and only invoke the  SolveLU  subroutine.

B = (({{-12.0}, {34.0}, {27.0}, {-19.0}}));

Solve the new system  AX = B for  X.

SolveLU[4];
\(\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)\)
Print["Solve AX = B"];
Print[MatrixForm[A0] , MatrixForm[vars], " = ", 
        MatrixForm[A0 . X]];
Print[" X = ", MatrixForm[vars], " = ", 
        MatrixForm[X]  ];
Print[""];
Print["AX = ", MatrixForm[A0] , MatrixForm[X], " = ", 
        MatrixForm[A0 . X]];
Print[" B = ", MatrixForm[B] ];

Application to Polynomial Curve Fitting

Theorem (Least-Squares Polynomial Curve Fitting).  Given the  n  data points  ((Subscript[x, 1]),(Subscript[y, 1])),((Subscript[x, 2]),(Subscript[y, 2])),...,((Subscript[x, n]),(Subscript[y, n])),  the least squares polynomial of degree  m  of the form  

		(Subscript[P, m]) (x)  =  (Subscript[c, 1])  +  
              (Subscript[c, 2]) x +  (Subscript[c, 3]) (x^2)  +   ...  + (Subscript[c, m]) 
            (x^(m - 1))  +  (Subscript[c, m + 1]) (x^m) 

that fits the n data points is obtained by solving the following linear system
 

		(({{n, (∑_(i = 1)^n) (Subscript[x, i]), \(∑\+\(i = 1\)\%n x\_i\^2\), ⋯, (∑_(i = 1)^n) ((Subscript[x, i])^m)}, {(∑_(i = 1)^n) (Subscript[x, i]), (∑_(i = 1)^n) ((Subscript[x, i])^2), (∑_(i = 1)^n) ((Subscript[x, i])^3), ⋯, (∑_(i = 1)^n) ((Subscript[x, i])^(m + 1))}, {(∑_(i = 1)^n) ((Subscript[x, i])^2), \(∑\+\(i = 1\)\%n x\_i\^3\), (∑_(i = 1)^n) ((Subscript[x, i])^4), ⋯, (∑_(i = 1)^n) ((Subscript[x, i])^(m + 2))}, {⋮, ⋮, ⋮, ⋱, ⋮}, {(∑_(i = 1)^n) ((Subscript[x, i])^m), \(∑\+\(i = 1\)\%n x\_i\^\(m + 1\)\), (∑_(i = 1)^n) ((Subscript[x, i])^(m + 2)), ⋯, (∑_(i = 1)^n) ((Subscript[x, i])^(2  m))}}))(({{ (Subscript[c, 1])}, { (Subscript[c, 2])}, { (Subscript[c, 3])}, {⋮}, { (Subscript[c, m + 1])}})) = (({{ (∑_(i = 1)^n) (Subscript[y, i])}, { (∑_(i = 1)^n) (Subscript[x, i]) (Subscript[y, i])}, { (∑_(i = 1)^n) ((Subscript[x, i])^2) (Subscript[y, i])}, {⋮}, {(∑_(i = 1)^n) ((Subscript[x, i])^m) (Subscript[y, i])}}))  

for the m+1 coefficients {(Subscript[c, 1]), (Subscript[c, 2]),  ... , (Subscript[c, m]), (Subscript[c, m + 1])}.  These equations are referred to as the "normal equations".

Proof  Least Squares Polynomials  Least Squares Polynomials

Example 6.  Find the "least squares parabola" that for the four data points  (2,8),(3,27),(4,64), and (5,125).

Solution 6.

(a). Write down the linear system AC = B to be solved.

XY  =  {{2, 8}, {3, 27}, {4, 64}, {5, 125}}; 
n = 4; 
x = (Subscript[(Transpose[XY]), [[1]]]); 
y = (Subscript[(Transpose[XY]), [[2]]]); 
A=(({{n, (∑_(k = 1)^n) (Subscript[x, [[k]]]), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^2)}, {(∑_(k = 1)^n) (Subscript[x, [[k]]]), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^2), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^3)}, {(∑_(k = 1)^n) ((Subscript[x, [[k]]])^2), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^3), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^4)}}));
B=(({{(∑_(k = 1)^n) (Subscript[y, [[k]]])}, {(∑_(k = 1)^n) 
                          (Subscript[x, [[k]]]) 
                      (Subscript[y, [[k]]])}, {(∑_(k = 1)^n)((((Subscript[x, [[k]]])))^2) (Subscript[y, [[k]]])}})); 
\(vars = {\*""(Subscript[c, 1])"", \*""(Subscript[c, 2])"", \*""(Subscript[c, 3])""};\)
Print["Solve AC = B"]; 
Print[MatrixForm[A], 
        MatrixForm[vars], " = ", MatrixForm[B]];

(b). Solve the linear system for the coefficients {(Subscript[c, i])} using our  LUfactor[n]  and  SolveLU[n]  subroutines.

LUfactor[3];
SolveLU[3];
c = Flatten[X]; 
Print[MatrixForm[vars], " = ", MatrixForm[c]];

(c). Construct the polynomial  p[x].  The coefficients are stored in the array  c  and the elements are (Subscript[c, [[i]]]).

Clear[x];
p[
          x_] =  (∑_(i = 1)^3) 
                (Subscript[c, [[i]]]) (x^(i - 1)); 
Print["p[x] = ", p[x]];

Of course we could do all this work in two lines by using Mathematica's built in  Fit[XY, {1, x, (x^2)}, x]  procedure.

q[x_] = Fit[XY, {1, x, (x^2)}, x];
Print["q[x] = ", q[x]];

We are done.

We can graph the polynomial, this is just for fun !

Needs["Graphics`Colors`"]; 
dots = 
        ListPlot[XY, PlotStyle -> {Red, PointSize[0.02]}, 
          DisplayFunction -> Identity]; 
gr1 = 
        Plot[p[x], {x, 0, 6}, PlotStyle -> Magenta, 
          DisplayFunction -> Identity];
Show[gr1, dots, PlotRange -> All, 
        DisplayFunction -> $DisplayFunction]; 
Print["The `least-squares` polynomial"]; 
Print["p[x] = ", p[x]]; 
Print["Points = ", XY];

Example 7.  Find the "least squares cubic" that for the four data points  (2,8),(3,27),(4,64), and (5,125).

Solution 7.

(a). Write down the linear system AC = B to be solved.

XY  =  {{2, 8}, {3, 27}, {4, 64}, {5, 125}}; 
n = 4; 
x = (Subscript[(Transpose[XY]), [[1]]]); 
y = (Subscript[(Transpose[XY]), [[2]]]); 
A=(({{n, (∑_(k = 1)^n) (Subscript[x, [[k]]]), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^2), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^3)}, {(∑_(k = 1)^n) (Subscript[x, [[k]]]), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^2), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^3), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^4)}, {(∑_(k = 1)^n) ((Subscript[x, [[k]]])^2), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^3), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^4), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^5)}, {(∑_(k = 1)^n) ((Subscript[x, [[k]]])^3), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^4), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^5), (∑_(k = 1)^n) ((Subscript[x, [[k]]])^6)}}));
B=(({{(∑_(k = 1)^n) (Subscript[y, [[k]]])}, {(∑_(k = 1)^n) 
                          (Subscript[x, [[k]]]) 
                      (Subscript[y, [[k]]])}, {(∑_(k = 1)^n)((((Subscript[x, [[k]]])))^2) (Subscript[y, [[k]]])}, {(∑_(k = 1)^n)((((Subscript[x, [[k]]])))^3) (Subscript[y, [[k]]])}})); 
A1 = 
        A;
\(vars = {\*""(Subscript[c, 1])"", \*""(Subscript[c, 2])"", \*""(Subscript[c, 3])"", \*""(Subscript[c, 4])""};\)
Print["Solve AC = B"]; 
Print[MatrixForm[A], 
        MatrixForm[vars], " = ", MatrixForm[B]];

(b). Solve the linear system for the coefficients {(Subscript[c, i])} using our  LUfactor[n]  and  SolveLU[n]  subroutines.

LUfactor[4];
SolveLU[4];
c = Flatten[X]; 
Print[MatrixForm[vars], " = ", MatrixForm[c]];

(c). Construct the polynomial  p[x].  The coefficients are stored in the array  c  and the elements are (Subscript[c, [[i]]]).

Clear[x];
p[
          x_] =  (∑_(i = 1)^4) 
                (Subscript[c, [[i]]]) (x^(i - 1)); 
Print["p[x] = ", p[x]];

Of course we could do all this work in two lines by using Mathematica's built in  Fit[XY, {1, x, (x^2), (x^3)}, x]  procedure.

q[x_] = Fit[XY, {1, x, (x^2), (x^3)}, x];
Print["q[x] = ", q[x]];

Notice that there appears to be a little "round-off" error creeping into the computation.  Lets "chop" it off.

q[x_] = Chop[q[x]]; 
Print["q[x] = ", q[x]];

We are done.

We can graph the polynomial, this is just for fun !

Needs["Graphics`Colors`"]; 
dots = 
        ListPlot[XY, PlotStyle -> {Red, PointSize[0.02]}, 
          DisplayFunction -> Identity]; 
gr1 = 
        Plot[p[x], {x, 0, 6}, PlotStyle -> Magenta, 
          DisplayFunction -> Identity];
Show[gr1, dots, PlotRange -> All, 
        DisplayFunction -> $DisplayFunction]; 
Print["The `least-squares` polynomial"]; 
Print["p[x] = ", p[x]]; 
Print["Points = ", XY];

Why is the "least-squares" polynomial a "perfect" power   p[x]  =  (x^3)?  Because the four data points are "perfect powers."
Of course we could do all this work in two lines by using Mathematica's built in  InterpolatingPolynomial[XY,x]  procedure.

q[x_] = 
        Expand[InterpolatingPolynomial[XY, x]];
Print["q[x] = ", q[x]];

The  InterpolatingPolynomial[XY,x]  procedure produces more accurate results if the answer should be an interpolating polynomial instead of a "least-squares" fit polynomial.

Caveat.  

	In numerical analysis we are concerned with round off error.  In the next example we will explore the "what if" scenario.  What will happen if the element  (Subscript[A, [[4, 
          4]]]) = 20514  is replaced with (Subscript[A, [[4, 
          4]]]) = 20515 in the matrix ?  This change is merely (1/20514) =  0.0000487  =  0.00487 %.  We will see that this will result is tremendous changes in the coefficients {(Subscript[c, i])}.  The numerical algorithms are correct, it is a problem with the matrix, it has a large condition number and is ill conditioned.  We need to be aware of these pitfalls.

Needs["LinearAlgebra`MatrixManipulation`"];
C1 = MatrixConditionNumber[N[A1]];
Print["A = ", MatrixForm[A1]];
Print["Condition number for A is ", C1];

Example 8.  Consider the linear system AC=B in Example 7, i.e.  (({{4, 14, 54, 224}, {14, 54, 224, 978}, {54, 224, 978, 4424}, {224, 978, 4424, 20514}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}, {(Subscript[c, 3])}, {(Subscript[c, 4])}})) = (({{224}, {978}, {4424}, {20514}})).

Replace the element  (Subscript[A, [[4, 
          4]]]) = 20514  with (Subscript[A, [[4, 
          4]]]) = 20515 in matrix A and solve the new system AC=B.

Solution 8.

(a). Enter the new matrix A and vector B.

XY  =  {{2, 8}, {3, 27}, {4, 64}, {5, 125}}; 
A=(({{4, 14, 54, 224}, {14, 54, 224, 978}, {54, 224, 978, 4424}, {224, 978, 4424, 20515}}));
B=(({{224}, {978}, {4424}, {20514}})); 
A2 = 
        A;
\(vars = {\*""(Subscript[c, 1])"", \*""(Subscript[c, 2])"", \*""(Subscript[c, 3])"", \*""(Subscript[c, 4])""};\)
Print["Solve AC = B"]; 
Print[MatrixForm[A], 
        MatrixForm[vars], " = ", MatrixForm[B]];

(b). Solve the linear system for the coefficients {(Subscript[c, i])} using our  LUfactor[n]  and  SolveLU[n]  subroutines.

LUfactor[4];
SolveLU[4];
c = Flatten[X]; 
Print[MatrixForm[vars], " = ", MatrixForm[c], " = ", 
        MatrixForm[N[c]]];

Remark.  The solution is quite different from that of Example 7.  This is not a problem in our subroutines, Mathematica will also get this answer.

A=(({{4, 14, 54, 224}, {14, 54, 224, 978}, {54, 224, 978, 4424}, {224, 978, 4424, 20515}}));
B=(({{224}, {978}, {4424}, {20514}})); 
c1 = 
        LinearSolve[A, B];
Print[MatrixForm[vars], " = ", MatrixForm[c1], " = ", 
          MatrixForm[N[c1]]];

These solution values are so far from the previous ones in Example 7 that they appear to be a worthless answer, and this type of ill-conditioning is not easy to detect.  One would normally expect that a small change in a coefficient in the matrix would result in a small change in the computed solution.  We changed the element  (Subscript[A, [[4, 4]]]) a small fraction,  (1/20514) =  0.0000487,  and the changes in the coefficients went 

from   (({{"\!\(c\_1\)"}, {"\!\(c\_2\)"}, {"\!\(c\_3\)"}, {"\!\(c\_4\)"}}))=(({{"0"}, {"0"}, {"0"}, {"1"}}))    to    (({{"\!\(c\_1\)"}, {"\!\(c\_2\)"}, {"\!\(c\_3\)"}, {"\!\(c\_4\)"}}))=(({{12.75`}, {-12.392857142857142`}, {3.75`}, {0.6428571428571429`}})).

The element  (Subscript[A, [[4, 
          4]]]) = 20514  was replaced with (Subscript[A, [[4, 
          4]]]) = 20515 in the matrix ?  This change is merely (1/20514) =  0.0000487  =  0.00487 %.  The percentage change in the solution (Subscript[c, 4]) is 36%.

Needs["LinearAlgebra`MatrixManipulation`"];
C2 = MatrixConditionNumber[N[A2]];
Print["A = ", MatrixForm[A2]];
Print["Condition number for A is ", C2];

Print["A = ", MatrixForm[A1]];
Print["For the `original matrix`"];
Print["Condition number for A is ", C1];

Continue with the explorations.

(c). Construct the polynomial  p[x].  The coefficients are stored in the array  c  and the elements are (Subscript[c, [[i]]]).

Clear[x];
q[
          x_] =  (∑_(i = 1)^4) 
                (Subscript[c, [[i]]]) (x^(i - 1)); 
Print["q[x] = ", q[x]];

We are done.

We can graph the polynomial, this is just for fun !

Needs["Graphics`Colors`"]; 
dots = 
        ListPlot[XY, PlotStyle -> {Red, PointSize[0.02]}, 
          DisplayFunction -> Identity]; 
gr2 = 
        Plot[q[x], {x, 0, 6}, PlotStyle -> Green, 
          DisplayFunction -> Identity];
Show[gr2, dots, PlotRange -> All, 
        DisplayFunction -> $DisplayFunction]; 
Print["The `least-squares` polynomial"]; 
Print["q[x] = ", q[x]]; 
Print["Points = ", XY];

Behold, do our eyes see anything significantly different from the previous graph in Example 7 ?

Show[gr1, gr2, dots, 
        DisplayFunction -> $DisplayFunction]; 
Print["p[x] = ", p[x]]; 
Print["q[x] = ", q[x]];

Show[gr1, gr2, dots, PlotRange -> {{2, 5}, {8, 125}}, 
        DisplayFunction -> $DisplayFunction]; 
Print["p[x] = ", p[x]]; 
Print["q[x] = ", q[x]];

If you were looking for the polynomial  p[x] = (x^3) and found the polynomial  q[x] = 
        (51/4) - ((347 x)/28) + ((15 (x^2))/4) + ((9 (x^3))/14)  would you be happy ?

Warning.  Remember to check the condition number of your matrix, beware if it is large and seek alternate ways to find compute a solution.  In approximation theory there are other methods and algorithms which can be used.

Old Lab Project (LU Factorization  LU Factorization).  Internet hyperlinks to an old lab project.

Research Experience for Undergraduates

LU Factorization  LU Factorization  Internet hyperlinks to web sites and a bibliography of articles.