Question

Write a MATLAB function [p, z] = proj(A,b) that computes the projection of b onto the Column Space of an m × n matrix A. Your program should allow for the possibility that the columns of A are not linearly independent. In order for the program to work, you will need to create a basis for Col A.

Answer 1

diary on

format compact

%Part IV. Orthogonal Projections & Least-Squares

Solutions. Exercise3_6

%shrink function

function B = shrink(A)

format compact,

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