write a MATLAB function p projAb which computes the project

write a MATLAB function p = proj(A,b)

which computes the projection of b onto the Column Space of an m x n matrix A. Your program should allow a possibility that the columns of A are not linearly independent. In order for the algorithm to work, you will need to create a basis for Col A.

Your function proj(A,b) should begin with

format compact,

A=shrink(A);

b=transpose(b);

Then, the program has to check if the input vector b has exactly m entries, where m is the number of columns in A-transpose. If it doesn’t, the program breaks with a message: “No solution: dimensions of A and b disagree” and returns “p=empty_vector”.

If b has exactly m components, proceed to the next step:

Determine whether it is the case that the vector b is an element of Col A. If b is an element of Col A, the program should have outputs (without any computations!): the vector p, the corresponding vector z, and a message “b is in the Col A”. After that, the program terminates.

If b is not an element of Col A, proceed to the next step.

Determine whether it is a case when b is orthogonal to Col A. If it is the case, the program should have outputs (without any computations!): the vector p, the corresponding vector z, and a message “b is orthogonal to Col A”. After that, the program terminates. If it is not the case, proceed to the next part:

Find the solution to the normal equations

Your output for this part should include:

vector p;

vector z = b – p, which is the component of b orthogonal to Col A;

One more command should be written to check whether the vector p is, indeed, orthogonal to the vector z = b – p - check whether the inner product of these vectors is zero. (For purpose of this program, we assume that a number is zero if its absolute value is less than 10^-7).

If your code confirms that p bottom z, please display a message: “Yes, p and z are orthogonal! Great Job!” Otherwise, the message should be something like: “Oops! Is there a bug in my code?”

Solution

diary on format compact %Part IV. Orthogonal Projections & Least-Squares Solutions. Exercise3_6 %shrink function function B = shrink(A) format compact, [~, pivot] = rref(A); B = A(: , pivot); %proj function function p = proj(A,b) format compact, A = shrink(A); b = transpose(b); if size(A, 1) == size(transpose(b), 2) if rank(A) == rank([A b]) p = transpose(b); disp(p) z = transpose(b) - p; disp(z) disp(\'b is in the Col A\') quit else R = colspace(sym(A)); T = 0; for i = 1:size(R,2) T = T + dot(R(:,i),b); end T = closetozeroroundoff(T); if T == 0 z = transpose(b); p = z - transpose(b); disp(p) disp(z) disp(\'b is orthogonal to Col A\') quit else p = (A * ((transpose(A) * A) \\ transpose(A)))*b; disp(p) z = transpose(b) - transpose(p); D = dot(p, z); if abs(D) < 1E-7 disp(\'Yes, p and z are orthogonal! Great Job!\') else disp(\'Oops! Is there a bug in my code?\') end end end else disp(\'No solution: dimensions of A and b disagree\') p = sym(\'empty_vector\'); %return end %Exercise3_6 (a) A=magic(6);A=A(:,1:4),b=(1:6) A = 35 1 6 26 3 32 7 21 31 9 2 22 8 28 33 17 30 5 34 12 4 36 29 13 b = 1 2 3 4 5 6 proj(A,b) 0.9492 2.1599 2.9492 3.9180 5.1287 5.9180 Yes, p and z are orthogonal! Great Job! ans = 0.9492 2.1599 2.9492 3.9180 5.1287 5.9180 %Exercise3_6 (b) A=magic(6),E=eye(6);b=E(6,:) A = 35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11 b = 0 0 0 0 0 1 proj(A,b) -0.2500 0.0000 0.2500 0.2500 0.0000 0.7500 Yes, p and z are orthogonal! Great Job! ans = -0.2500 0.0000 0.2500 0.2500 0.0000 0.7500 %Exercise3_6 (c) A=magic(4),b=(1:5) A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 b = 1 2 3 4 5 proj(A,b) No solution: dimensions of A and b disagree ans = empty_vector %Exercise3_6 (d) A=magic(5),b=rand(1,5) A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 b = 0.8147 0.9058 0.1270 0.9134 0.6324 proj(A,b) 0.8147 0.9058 0.1270 0.9134 0.6324 0 0 0 0 0 b is in the Col A diary on format compact %Exercise3_6 (e) A=ones(6);A(:)=1:36,b=[1,0,1,0,1,0] A = 1 7 13 19 25 31 2 8 14 20 26 32 3 9 15 21 27 33 4 10 16 22 28 34 5 11 17 23 29 35 6 12 18 24 30 36 b = 1 0 1 0 1 0 proj(A,b) 0.7143 0.6286 0.5429 0.4571 0.3714 0.2857 Yes, p and z are orthogonal! Great Job! ans = 0.7143 0.6286 0.5429 0.4571 0.3714 0.2857 %Exercise3_6 (f) A=ones(6);A(:)=1:36;A=null(A,\'r\'),b=ones(1,6) A = 1 2 3 4 -2 -3 -4 -5 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 b = 1 1 1 1 1 1 proj(A,b) 0 0 0 0 0 0 1 1 1 1 1 1 b is orthogonal to Col A %Exercise3_7 %shrink function function B = shrink(A) format compact, [~, pivot] = rref(A); B = A(: , pivot); %solvemore function function X = solvemore(A,b) format long, A = shrink(A); [~,n] = size(A); if rank(A) == rank([A b]) disp(\'The equation is consistent--look for the exact solution\') B = closetozeroroundoff(A\' * A - eye(n)); if B == zeros(size(B)) if size(A,1) == size(A,2) x1 = A \\ b; x2 = transpose(A)*b; N = norm(x1 - x2); disp(\'A is orthogonal\') X = [ x1, x2 ]; disp(X) disp(\'The norm of the difference between two solutions is N =\') disp(N) else x1 = A \\ b; disp(\'A has orthonormal columns but is not orthogonal\') X = x1; disp(X) end else x1 = A \\ b; disp(\'A does not have orthogonal columns\') X = x1; disp(X) end else disp(\'The system is inconsistent--look for the least-squares solution\') x3 = inv(transpose(A)*A)*transpose(A)*b; n1 = norm(b - A * x3); disp(\'The solution of the normal equations is\') disp(x3) disp(\'The least-squares error of the approximation is n1 = \') disp(n1) B = closetozeroroundoff(A\' * A - eye(n)); if B == zeros(size(B)) disp(\'A has orthonormal columns: an orthonormal basis for Col A is U = A\') U = A; else U = orth(A); disp(\'An orthonormal basis for Col A, U = \') disp(U) b1 = U*transpose(U)*b; disp(\'The projection of b onto Col A is\') disp(b1) x3 = inv(transpose(A)*A)*transpose(A)*b; x4 = A \\ (b1); disp(\'The least-squares solution by using the projection onto Col A is x4 =\') disp(x4) n2 = norm(b - A * x4); disp(\'The least-squares error os this approximation is n2 = \') disp(n2) n3 = norm(x3 - x4); x = rand(n,1); n4 = norm(b - A*x); disp(\'An error of approximation of b by Ax for a random vector x in the real coordinate space is\') disp(n4) X = [x3,x4]; end end diary on format compact %Exercise3_7 (a) A=magic(4);b=A(:,1),A=orth(A) b = 16 5 9 4 A = -0.5000 0.6708 0.5000 -0.5000 -0.2236 -0.5000 -0.5000 0.2236 -0.5000 -0.5000 -0.6708 0.5000 solvemore(A,b) The equation is consistent--look for the exact solution A has orthonormal columns but is not orthogonal -16.999999999999993 8.944271909999161 3.000000000000000 ans = -16.999999999999993 8.944271909999161 3.000000000000000 %Exercise3_7 (b) A=magic(5);A=orth(A),b=rand(5,1) A = Columns 1 through 3 -0.447213595499958 -0.545634873129949 0.511667273601712 -0.447213595499958 -0.449758363151204 -0.195439507584839 -0.447213595499958 -0.000000000000023 -0.632455532033676 -0.447213595499958 0.449758363151190 -0.195439507584871 -0.447213595499958 0.545634873129986 0.511667273601673 Columns 4 through 5 0.195439507584855 -0.449758363151197 -0.511667273601693 0.545634873129968 0.632455532033676 0.000000000000000 -0.511667273601693 -0.545634873129968 0.195439507584854 0.449758363151197 b = 0.814723686393179 0.905791937075619 0.126986816293506 0.913375856139019 0.632359246225410 solvemore(A,b) The equation is consistent--look for the exact solution A is orthogonal -1.517501961599936 -1.517501961599937 -0.096093467150122 -0.096093467150122 0.304574206628231 0.304574206628231 -0.567677934732546 -0.567677934732546 -0.086157982822827 -0.086157982822827 The norm of the difference between two solutions is N = 1.035640085178848e-15 ans = -1.517501961599936 -1.517501961599937 -0.096093467150122 -0.096093467150122 0.304574206628231 0.304574206628231 -0.567677934732546 -0.567677934732546 -0.086157982822827 -0.086157982822827 %Exercise3_7 (c) A=magic(4),b=ones(4,1) A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 b = 1 1 1 1 solvemore(A,b) The equation is consistent--look for the exact solution A does not have orthogonal columns 0.058823529411765 0.117647058823529 -0.058823529411765 ans = 0.058823529411765 0.117647058823529 -0.058823529411765 %Exercise3_7 (d) A=magic(4),b=rand(4,1) A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 b = 0.097540404999410 0.278498218867048 0.546881519204984 0.957506835434298 solvemore(A,b) The system is inconsistent--look for the least-squares solution The solution of the normal equations is -0.001240611967882 -0.031004149639098 0.087551437445384 The least-squares error of the approximation is n1 = 0.372331330756435 An orthonormal basis for Col A, U = -0.363225569906992 -0.839773278980323 0.335928601456289 -0.511952614823082 0.201051551215513 -0.497476425501395 -0.413098103234259 -0.228706948321817 -0.571876812690966 -0.659789104673461 0.449502219631666 0.559129763025005 The projection of b onto Col A is 0.180796221571844 0.528265668584352 0.297114069487679 0.874251018861863 The least-squares solution by using the projection onto Col A is x4 = -0.001240611967882 -0.031004149639099 0.087551437445385 The least-squares error os this approximation is n2 = 0.372331330756435 An error of approximation of b by Ax for a random vector x in the real coordinate space is 34.842978007481335 ans = -0.001240611967882 -0.001240611967882 -0.031004149639098 -0.031004149639099 0.087551437445384 0.087551437445385 %Exercise3_7 (e) A=magic(4);A=orth(A),b=rand(4,1) A = -0.500000000000000 0.670820393249937 0.500000000000000 -0.500000000000000 -0.223606797749979 -0.500000000000000 -0.500000000000000 0.223606797749979 -0.500000000000000 -0.500000000000000 -0.670820393249937 0.500000000000000 b = 0.957166948242946 0.485375648722841 0.800280468888800 0.141886338627215 solvemore(A,b) The system is inconsistent--look for the least-squares solution The solution of the normal equations is -1.192354702240901 0.617321717584815 -0.093301415370740 The least-squares error of the approximation is n1 = 0.028942288916205 A has orthonormal columns: an orthonormal basis for Col A is U = A diary off