213 - A16 - Least Squares

147 days ago by Professor213

# Do this page in Octave mode 
       
A = floor(rand(6,4)*10)-5 b = floor(rand(6,1)*10)-5 
       
A =

 2 -4 -2 1
 3 -3 -3 -1
 -5 4 0 3
 1 4 -4 3
 -2 -1 -1 -4
 0 1 3 -2


b =

 0
 -5
 -1
 3
 0
 -2
A =

 2 -4 -2 1
 3 -3 -3 -1
 -5 4 0 3
 1 4 -4 3
 -2 -1 -1 -4
 0 1 3 -2


b =

 0
 -5
 -1
 3
 0
 -2
# The goal is to solve Ax = b but likely this system is over-determined and therefore inconsistent. Let's see... format short mat = [A b] rref(mat) 
       
mat =

   2  -4  -2   1   0
   3  -3  -3  -1  -5
  -5   4   0   3  -1
   1   4  -4   3   3
  -2  -1  -1  -4   0
   0   1   3  -2  -2

ans =

   1   0   0   0   0
   0   1   0   0   0
   0   0   1   0   0
   0   0   0   1   0
   0   0   0   0   1
   0   0   0   0   0
mat =

   2  -4  -2   1   0
   3  -3  -3  -1  -5
  -5   4   0   3  -1
   1   4  -4   3   3
  -2  -1  -1  -4   0
   0   1   3  -2  -2

ans =

   1   0   0   0   0
   0   1   0   0   0
   0   0   1   0   0
   0   0   0   1   0
   0   0   0   0   1
   0   0   0   0   0
# So, let's do Gram-Schmidt and generate a QR factorization v1 = A(:,1); v2 = A(:,2); v3 = A(:,3); v4 = A(:,4); disp("Here are the original vectors") V = [v1 v2 v3 v4]; printf ("% 4.0f % 4.0f % 4.0f % 4.0f \n", V'); disp("Now apply Gram-Schmidt") u1 = v1; u2 = v2 - (v2'*u1)/(u1'*u1)*u1; u3 = v3 - (v3'*u1)/(u1'*u1)*u1 - (v3'*u2)/(u2'*u2)*u2; u4 = v4 - (v4'*u1)/(u1'*u1)*u1 - (v4'*u2)/(u2'*u2)*u2 - (v4'*u3)/(u3'*u3)*u3; Q = [u1 u2 u3 u4]; printf ("% 5.3f % 5.3f % 5.3f % 5.3f \n", Q'); R = Q'*A 
       
Here are the original vectors
   2   -4   -2    1 
   3   -3   -3   -1 
  -5    4    0    3 
   1    4   -4    3 
  -2   -1   -1   -4 
   0    1    3   -2 
Now apply Gram-Schmidt
 2.000 -2.558 -1.708  2.136 
 3.000 -0.837 -2.086 -0.883 
-5.000  0.395 -1.681  1.607 
 1.000  4.721 -2.902 -0.643 
-2.000 -2.442 -2.085 -3.528 
 0.000  1.000  3.159 -1.492 
R =

   43.00000  -31.00000  -15.00000   -5.00000
   -0.00000   36.65116   -5.81395   21.39535
    0.00000   -0.00000   32.84518  -11.35025
    0.00000   -0.00000   -0.00000   23.00665
Here are the original vectors
   2   -4   -2    1 
   3   -3   -3   -1 
  -5    4    0    3 
   1    4   -4    3 
  -2   -1   -1   -4 
   0    1    3   -2 
Now apply Gram-Schmidt
 2.000 -2.558 -1.708  2.136 
 3.000 -0.837 -2.086 -0.883 
-5.000  0.395 -1.681  1.607 
 1.000  4.721 -2.902 -0.643 
-2.000 -2.442 -2.085 -3.528 
 0.000  1.000  3.159 -1.492 
R =

   43.00000  -31.00000  -15.00000   -5.00000
   -0.00000   36.65116   -5.81395   21.39535
    0.00000   -0.00000   32.84518  -11.35025
    0.00000   -0.00000   -0.00000   23.00665
# And now, let's try to solve the original Ax=b problem by using QRx = b y = Q'*b x = inv(R)*y # and also by using A'Ax = A'b xother = inv(A'*A)*A'*b 
       
y =

   -7.0000
   15.9535
   -2.9112
    3.8671

x =

   0.085672
   0.332312
  -0.030548
   0.168086

xother =

   0.085672
   0.332312
  -0.030548
   0.168086
y =

   -7.0000
   15.9535
   -2.9112
    3.8671

x =

   0.085672
   0.332312
  -0.030548
   0.168086

xother =

   0.085672
   0.332312
  -0.030548
   0.168086
# Let's see how good we did by computing the residual r = norm(b - A*x) 
       
r =  5.4780
r =  5.4780
A\b 
       
ans =

   0.085672
   0.332312
  -0.030548
   0.168086
ans =

   0.085672
   0.332312
  -0.030548
   0.168086