# 213 - A16 - Least Squares

## 836 days ago by Professor213

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