213 - A62 - Least Squares

29 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 =

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

b =

 -4
 -4
 4
 -4
 -2
 4
A =

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

b =

 -4
 -4
 4
 -4
 -2
 4
# 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 =

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

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 =

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

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
   1    3   -3    0 
  -1   -5   -4   -2 
  -1   -1    2   -1 
  -1   -5    3    1 
  -3    3   -2   -4 
   0   -4    0    1 
Now apply Gram-Schmidt
 1.000  2.615 -2.752  0.272 
-1.000 -4.615 -4.556 -0.347 
-1.000 -0.615  2.059 -0.684 
-1.000 -4.615  2.444  0.454 
-3.000  4.154 -0.900  0.283 
 0.000 -4.000 -0.615  0.452 
R =

   1.3000e+01   5.0000e+00   2.0000e+00   1.4000e+01
  -3.5527e-15   8.3077e+01  -1.2769e+01  -1.5385e+01
  -4.4409e-16   4.4409e-15   3.9730e+01   1.2481e+01
   2.2204e-16   5.1070e-15   4.3299e-15   1.1529e+00
Here are the original vectors
   1    3   -3    0 
  -1   -5   -4   -2 
  -1   -1    2   -1 
  -1   -5    3    1 
  -3    3   -2   -4 
   0   -4    0    1 
Now apply Gram-Schmidt
 1.000  2.615 -2.752  0.272 
-1.000 -4.615 -4.556 -0.347 
-1.000 -0.615  2.059 -0.684 
-1.000 -4.615  2.444  0.454 
-3.000  4.154 -0.900  0.283 
 0.000 -4.000 -0.615  0.452 
R =

   1.3000e+01   5.0000e+00   2.0000e+00   1.4000e+01
  -3.5527e-15   8.3077e+01  -1.2769e+01  -1.5385e+01
  -4.4409e-16   4.4409e-15   3.9730e+01   1.2481e+01
   2.2204e-16   5.1070e-15   4.3299e-15   1.1529e+00
# 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 =

    6.00000
   -0.30769
   27.02963
   -3.01016

x =

   3.14118
  -0.25657
   1.50061
  -2.61098

xother =

   3.14118
  -0.25657
   1.50061
  -2.61098
y =

    6.00000
   -0.30769
   27.02963
   -3.01016

x =

   3.14118
  -0.25657
   1.50061
  -2.61098

xother =

   3.14118
  -0.25657
   1.50061
  -2.61098
# Let's see how good we did by computing the residual r = norm(b - A*x) 
        
r =  7.4149
r =  7.4149
A\b 
        
ans =

   3.14118
  -0.25657
   1.50061
  -2.61098
ans =

   3.14118
  -0.25657
   1.50061
  -2.61098