A1 Solutions
P1
function B = MatProd(A)
% A is a length t cell array with the property that A{1},...,A{t} are
% matrices and A{1}*A{2}*...*A{t} is defined
% B = A{1}*A{2}*...*A{t}
global OpCount
t = length(A);
for k=1:t
[m(k),n(k)]= size(A{k});
end
if t==1
B = A{1};
elseif t==2
B = A{1}*A{2};
OpCount = OpCount + 2*m(1)*n(1)*n(2);
else
[nmin,p] = min(n(1:t-1));
Aleft = cell(p,1);
for k=1:p
Aleft{k} = A{k};
end
Aright = cell(t-p,1);
for k=1:t-p
Aright{k} = A{p+k};
end
B = MatProd(Aleft)*MatProd(Aright);
OpCount = OpCount + 2*m(1)*n(p)*n(t);
end
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 2 pts for correct implementation of D&C method given by the question.
% 1 pt for suggesting and implementing an improved version. Changing
% a few statements do not count.
% 1 pt if the suggested improve version gives better answer than
% the D&C method in every test case.
% 1 pt for submitting all required plots.
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
% Below is one of the interesting non-dynamic programming answer for
% the improved version. This is due to Bernard Wong.
function B = MatProdGreedy(A)
% As described in the assignment:
% A is a length t cell array with the property that A{1},...,A{t} are
% matrices and A{1}*A{2}*...*A{t} is defined
% B = A{1}*A{2}*...*A{t}
global OpCount;
% First determine the size of the block matrix
[NumRowsInBlock NumColsInBlock] = size(A);
% Error case, return nothing
% We expect A to be a column cell array
if NumColsInBlock ~= 1 || NumRowsInBlock == 0
B=[];
return
end
% If there are only 1 or 2 blocks left, there we are at the end of the
% recursion and we should just calculate either A{1,1}*A{2,1} for 2
% blocks or A{1,1} for 1 block
if NumRowsInBlock == 1
B = A{1,1};
return
elseif NumRowsInBlock == 2
B = A{1,1} * A{2,1};
% Count the ops
[s q] = size(A{1,1});
[q r] = size(A{2,1});
OpCount = OpCount + 2*s*q*r;
return
end
% Iterative approach. Iterate through the arrays to see which consecutive
% array pair will create the smallest resulting matrix size after
% multiplication. If there is a tie, choose the one with the largest inner
% indices (i.e. m1 and n2), as this can reduce the size of subsequent
% matrices that are created.
% First find the cheapest pair(s)
SmallestSizeProduct = size(A{1,1}, 1) * size(A{2,1}, 2);
ArrayOfSmallest = [1];
for i=2:NumRowsInBlock-1
CurrentSizeProduct = size(A{i,1},1)*size(A{i+1,1},2);
if CurrentSizeProduct == SmallestSizeProduct
ArrayOfSmallest = [ArrayOfSmallest i];
elseif CurrentSizeProduct < SmallestSizeProduct
ArrayOfSmallest = [i];
SmallestSizeProduct = CurrentSizeProduct;
end
end
% Out of the cheapest pairs, find the cheapest base it it(sic) having
% the largest inner indices
BiggestInnerOfArrayOfSmallest = size(A{ArrayOfSmallest(1,1)},2);
IndexOfSmallestArray = 1;
for i=2:size(ArrayOfSmallest, 2)
if size(A{ArrayOfSmallest(1,i)}, 2) < BiggestInnerOfArrayOfSmallest
IndexOfSmallestArray = i;
end
end
p = ArrayOfSmallest(1,IndexOfSmallestArray);
% Now perform the multiplication
Aproduct = A{p,1}*A{p+1,1};
% Count the ops
[s q] = size(A{p,1});
[q r] = size(A{p+1,1});
OpCount = OpCount + 2*s*q*r;
% And recursively call itself
if p == 1
B = MatProdGreedy({Aproduct A{p+2:NumRowsInBlock,1}}');
elseif p == NumRowsInBlock - 1
B = MatProdGreedy({A{1:p-1,1} Aproduct}');
else
B = MatProdGreedy({A{1:p-1,1} Aproduct A{p+2:NumRowsInBlock,1}}');
end
%-----------------------------------------------------------------------
% Below is one of the clear and well-commented dynamic programming
% solution. This is due to Michael Jennings.
%
function B = DynMatProd(A)
% A is a length t cell array with the property that A{1},...,A{t} are
% matrices and A{1}*A{2}*...*A{t} is defined.
% B = A{1}*A{2}*...*A{t}
% Let Opt(i,j) denote the minimum cost of multiplying matrices A{i},
% ...,{j}. Let r(i) and c(i) be the # of rows and # of columns
% respectively of A{i}. The dynamic programming recurrence relation is:
% Opt(i,j) = min{ Opt(i,k) + Opt(k+1,j) + 2*r(i)*c(k)*c(j)}, where the
% minimum is taken over k=i:j-1. The base cases are:
% Opt(i,i) = 0, and Opt(i,i+1) = 2*r(i)*c(i)*c(i+1).
% This is an O(t^3) algorithm, and so it does not involve much more
% work than the O(t^2) DC algorithm.
global OpCount;
global Opt;
len = size(A,1);
for i=1:len
Opt(i,i)=0;
r(i) = size(A{i},1);
c(i) = size(A{i},2);
partition(i,i) = i;
end
for i=1:len-1
Opt(i,i+1) = 2*r(i)*c(i)*c(i+1);
partition(i,i+1) = i+1;
end
for i=len-2:-1:1
for j=i+2:len % i+1 is already stored
min = inf; % reset min
for k=i:j-1
temp = Opt(i,k) + Opt(k+1,j)+2*r(i)*c(k)*c(j);
if min > temp
min = temp;
% keep track of which k minimizes Opt
partition(i,j) = k;
end
end
Opt(i,j) = min;
end
end
OpCount = Opt(1,len);
% send a complete description of the problem
% to a simple recursive procedure
B = Mult(A, partition, r, c, 1, len, len);
% end of Michael's quote.
%
% Mult.m performs the multiplication according to the division
% found by DynMatProd() method. The code for Mult.m is not cited here.
P2
% P2: Floating Point Numbers
% Since 10^j = (5^j)*(2^j), we see that the mantissa must be long enough so that
% it can hold 5^j.
% There is no point considering values of j that are bigger than jmax where
jmax = ceil(64*log(2)/log(5)) % Smallest j so 5^j > 2^64
% How many bits are required to represent 5^j where j=1:jmax?
alpha = floor(log2(5.^(1:jmax)'))+1; % alpha(j) = # bits in binary expansion of 5^j.
% If 2^i-1 <= 5^j < 2^i, then 10^j = (5^j / 2^i) * 2^(i+j) and so the exponent must be long
% enough so that it can hold the base-2 representation of (i+j).
beta = floor(log2(alpha + (1:jmax)'))+1; % beta(j) = # bits required for the exponent of 10^j
% To store 10^j we need a mantissa that is at least alpha(j) bits long and an
% exponent that is at least beta(j) bits long.
clc
disp(' t k')
disp('----------')
for t=1:63
% t = number of bits allocated for the mantissa
% Determine largest k so that the floating point representation for 10^k
% can fit in a floating point word that has a t-bit mantissa and a 64-t bit exponent.
% That is, find the largest k so that alpha(k) <= t AND beta(k) <= 64 - t.
k = sum(alpha<=t & beta<=64-t);
disp(sprintf(' %2d %2d',t,k))
end
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 3 pts for correct calculation of # of mantissa bits needed.
% 2 pts for correct calculation of # of exponent bits needed.
%
% Common mistakes:
% 1) Thinking that e bits are needed to store 2^e. In fact, you need
% e+1 bits to store 2^e. This mistake is mostly in the form of
% coding that ceil(log2(<# of mantissa bits> + k)) exponent bits are
% needed for 10^k. The correct one is floor(log2(<# of mantissa bits> + k))+1.
% 1 pt is deducted if this is the only mistake.
%
% 2) Do not take into account that in order to make mantissa into the
% correct form (0.xxxxx) in base 2, extra exponents are needed to shift
% the decimal point to correct place. For example,
% 10 = 5 * 2 = (101)_2 * 2 = (0.101) * 2^4. Solutions containing
% this mistake usually code that floor(log2(k))+1 bits are needed
% for the exponent instead of floor(log2(<# of mantissa bits> + k))+1.
% 2 pts are deducted for this mistake.
%-----------------------------------------------------------------------
P3
function K = HamSq(H) % H is a 2n-by-2n Hamiltonian Matrix % K = H^2 % Extract A, F, and G... [m,m] = size(H); n = m/2; A = H(1:n,1:n); F = H(1:n,n+1:2*n); G = H(n+1:2*n,1:n); % [A F ; G -A']*[A F ; G -A'] = [K11 K12; K21 K22] where % K11 = A*A + F*F, K12 = A*F - F*A', K21 = G*A -A'*G , K22 = G*G + A'*A' K11 = A*A + F*G; AF = A*F; GA = G*A; K = [K11 AF-(AF)' ; GA-(GA)' K11']; %----------------------------------------------------------------------- % GRADING SCHEME % -------------- % 2 pts for calculating K efficiently. % 2 pts for vectorizing code completely. 1 pt is deducted for each % nested for-loop. % 1 pt for submitting required output, the resulting K's for each % method and the plot. %-----------------------------------------------------------------------
P4
function theta = MakeSing(A,p,q)
% A is an n-by-n matrix with unit column sums and entries that satisfy 0 <= A(i,j) < 1.
% p and q are integers that satisfy 1 <= p,q <= n.
% The matrix A^(pq)(theta) is singular.
[n,n] = size(A);
eq = zeros(n,1); eq(q) = 1;
xq = A'\eq;
theta = (1 - A(p,q)*xq(p))/(A(p,q) - A(p,q)*xq(p));
%-----------------------------------------------------------------------
% GRADING SCHEME
% --------------
% 3 pts for computing theta correctly.
% 2 pts for efficiency. Mainly, solve linear system when needed
% without forming inverse of A explicitly.
%-----------------------------------------------------------------------