A2 Solutions
P1
function [v,sigma] = StatVec1(theta)
% theta is a column n-vector with 0 <= theta(1) < ... < theta(n) < 2pi
% v is the normalized stationary vector for P(theta)
% sigma is the column n-vector of the singuar values of P(theta).
% Uses the SVD method;
n = length(theta);
P = SetUp(theta);
A = eye(n,n) - P;
[U,S,V] = svd(A);
v = V(:,n)/sum(V(:,n));
sigma = diag(S);
%-----------------------------------------------------------------------
% GRADING SCHEME: See under SetUp.m
%-----------------------------------------------------------------------
function [v,UCond] = StatVec2(theta)
% theta is a column n-vector with 0 <= theta(1) < ... < theta(n) < 2pi
% v is the normalized stationary vector for P(theta)
% Ucond is the condition of U(1:,n-1,1:n-1) where (I - P(theta)) = LU
% Uses the LU method
n = length(theta);
P = SetUp(theta);
A = eye(n,n) - P;
[L,U] = lu(A);
v = [-U(1:n-1,1:n-1)\U(1:n-1,n) ; 1];
v = v/sum(v);
UCond = cond(U(1:n-1,1:n-1));
%-----------------------------------------------------------------------
% GRADING SCHEME: See under SetUp.m
%-----------------------------------------------------------------------
function P = SetUp(theta)
% theta is a column n-vector with 0 <= theta(1) < ... < theta(n) < 2pi
% P is the n-by-n stochastic matrix P(theta)
n = length(theta);
v = [cos(theta) sin(theta)];
%Set up matrix of inverse pairwise distances. D(i,j) = 1/distance between ith and jth point
D = zeros(n,n);
for i=1:n
for j = 1:i-1
D(i,j) = 1/norm(v(i,:)-v(j,:));
D(j,i) = D(i,j);
end
end
% The Markov matrix..
P = zeros(n,n);
for j=1:n
P(:,j) = [D(1:j-1,j);0;D(j+1:n,j)]/sum(D(:,j));
end
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 1 pts for correctness of SetUp()
% 2 pts for correctness of StatVec1()
% 2 pts for correctness of StatVec2()
%-----------------------------------------------------------------------
P2
function d = PosDefD(A)
% A is an n-by-n symmetric positive definite matrix.
% d is a column n-vector with the property that A(1:k,1:k) - d(k)ek*ek'
% is singular where ek is the last column of eye(k,k).
G = chol(A)';
% Compare blocks in [A11 A12; A21 A22] = [G11 0; G21 G22]*[G11 0; G21 G22]'
% where A11 is k-by-k and conclude that
% G(1:k,1:k) is the Cholesky factor of A(1:k,1:k).
[n,n] = size(A);
d = zeros(n,1);
for k=1:n
% If we reduce A(k,k) by d(k) then when we compute G(k,k) in Cholesky we
% have this recipe:
% G(k,k) = sqrt( (A(k,k)-d(k) ) - (G(k,1)^2 + ... + G(k,k-1)^2) )
% Thus, by setting
d(k) = A(k,k) - G(k,1:k-1)*G(k,1:k-1)';
% We make G(1:k,1:k) singular.
end
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 2 pts for showing that there is a relationship between d(k) and Chol(A)
% 3 pts for correctness
%-----------------------------------------------------------------------
P3
function X = Solve(A)
% A is a n-by-(n+1)-by-m array such that A(1:n,1:n,k) is nonsingular, k=1:m
% X is a n-by-m with the property that A(1:n,1:n,k)X(:,k) = A(:,n+1,k), k=1:m
[n,np1,m] = size(A);
for k=1:n-1
% k-th step of Gaussian elimination
% Pivoting...
for p=1:m
% In the pth system, Swap row k with row q
v = reshape(A(k:n,k,p),n-k+1,1);
[mu,q0] = max(abs(v)); q = q0+k-1;
A([k,q],k:n+1,p) = A([q k],k:n+1,p);
end
% Eliminating the kth variable from equations k+1 through n...
for i=k+1:n
tau = A(i,k,:)./A(k,k,:);
for j=k:n+1
A(i,j,:) = A(i,j,:) - tau.*A(k,j,:);
end
end
end
% Vectorized back-substitution...
X = zeros(n,m);
for i = n:-1:1
X(i,:) = reshape(A(i,n+1,:),1,m);
for j=i+1:n
X(i,:) = X(i,:) - reshape(A(i,j,:),1,m).*X(j,:);
end
X(i,:) = X(i,:)./reshape(A(i,i,:),1,m);
end
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 3 pts for vectorizing backsubstitution correcly
% 2 pts for correct pivoting
%-----------------------------------------------------------------------
P4
function D = Tetra(theta,phi)
% theta and phi are column n-vectors with distinct compponents that satisfy
% with -90 <= theta(i) <= 90 and -180 <= phi(i) <= 180.
%
% Let C(i) be the point on the sphere of radius 3950 with latitude theta(i) and longitude phi(i),
% both in degrees. Assume that no C(i) is an antipode of another.
%
% Let T(p,q,r,s)) be the tetrahedron defined by the
% tangent planes to C(p), C(q), C(r), and C(s).
%
% D is an n-by-n-by-n-by-n array where for all p, q, r, and s that satisfy n >= p > q > r > s >= 1,
% D(p,q,r,s) is the diameter of T(p,q,r,s).
% Compute Cartesian coordinates of the surface points.
% (We work with the unit sphere until the very end when we scale the
% diameters by R.)
R = 3950;
theta = (pi/180)*theta; phi = (pi/180)*phi;
cTheta = cos(theta); sTheta = sin(theta);
cPhi = cos(phi); sPhi = sin(phi);
P = [-cTheta.*sPhi -cTheta.*cPhi sTheta]; % The i-th surface point is (P(i,1),P(i,2),P(i,3))
% Set up and solve the linear systems that define the m intersections of the tangent planes
% m = n(n-1)(n-2)/6 = n-choose-3.
n = length(theta);m = n*(n-1)*(n-2)/6;
X = cell(n,n,n);
for i1 = 1:n
for i2=1:i1-1
for i3 = 1:i2-1
xVec = [P(i1,:) ; P(i2,:) ;P(i3,:) ]\ones(3,1); % The intersection of T(i1), T(i2), and T(i3)
X{i1,i2,i3}= struct('sol',xVec,'norm',norm(xVec)); % Save solution and its norm in X(i1,i2,i3)
end
end
end
% Check all the tetrahedrons...
D = zeros(n,n,n,n);
for p = 1:n
for q = 1:p-1
for r = 1:q-1
for s = 1:r-1
% Is [0;0;0;0] a convex combination of the tetrahedron's four vertices
% X{p,q,r}.sol X{p,q,s}.sol X{p,r,s}.sol X{q,r,s}.sol?
% If so, then the tetrahedron encloses the unit sphere.
% To answer the convex combination question, we need to know if
% [X(p,q,r).sol X(p,q,s).sol X(p,r,s).sol X(q,r,s).sol]
% has a non-negative vector in its null space.
alfa = [X{p,q,r}.sol X{p,q,s}.sol X{p,r,s}.sol]\(-X{q,r,s}.sol);
if alfa >= 0
D(p,q,r,s) = R*max([X{p,q,r}.norm X{p,q,s}.norm X{p,r,s}.norm X{q,r,s}.norm]);
else
D(p,q,r,s) = -1;
end
end
end
end
end
%----------------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 3 pts for correctness of the algorithm.
% 2 pts for efficiency. See common mistake below.
%
% COMMON MISTAKE NOTE (1)
% -----------------------
% Most students implement this problem by doing 4 nested for-loop and
% then for each iteration of tuple (p,q,r,s), each vertex of the tetrahedron
% is computed by finding the intersection of 3 planes. The problem with this
% approach is that same equations are solved many times. For example,
% the intersection point of planes T(p), T(q), and T(r) are the same for
% all tuples (p,q,r,s) for every s. Therefore, the correct approach is
% to compute the intersection point of any three planes once and store
% the result. Then, when deciding whether (p,q,r,s) circumscribes the sphere,
% these saved values are retrieved for further computation. This is an order-
% of-magnitude difference
%-----------------------------------------------------------------------------