ACM11 Homework 1 Solutions

Problem 1: simple substitution cypher
Problem 2 (original): vectorization
Problem 2 (version 2): vectorization
Problem 3
Problem 4: Solving linear systems (accuracy)
Problem 5: Solving linear systems (time)

Problem 1: simple substitution cypher

Problem: decode the message. Below is how the message was encoded:

plainText = 'All work and no play makes Jack a dull boy';
k = 37;
cipherText = mod( plainText + k, 128 );

To recover it, using the fact that the original plain text contains the word "work", we can do the following (first, import the data into the variable "cipherText" and make sure "cipherText" is a row vector):

for kk = 1:128
    test = mod( cipherText - kk, 128 );
    if strfind(test,'work')
        recoveredText = test;
        disp(['At k = ',num2str(kk),' found the word ''work''']);
        break;
    end
end
disp('The recovered message is:');
disp( char(recoveredText) )

At k = 37 found the word 'work'
The recovered message is:
All work and no play makes Jack a dull boy

so we find that k = 37. A value of k = 91 is also acceptable, since perhaps it wasn't clear if we shifted UP by k or DOWN by k (and 91 = -37 modulo 128 ).

Problem 2 (original): vectorization

There were two versions of this, after we noticed different results on different computers and versions of MATLAB. Here's what the code should look like:

fprintf('On my computer, a 3 GHz Pentium 4 running Redhat linux,\n')
fprintf('  and with MATLAB version %s, the results are:\n',version);
randn('state',2008);
N  = 1e7;
a = randn(N,1); b = randn(N,1);
disp('... via my ''for'' loop:')
tic; s = 0;
for i = 1:N
    s = s + a(i)*b(i);
end
toc

disp('... via the command a''*b:')
tic
s2 = a'*b;
toc
disp(['Difference in answers is ',num2str(abs(s-s2)) ]);
disp('... and via command dot(a,b):');
tic
s3 = dot(a,b);
toc

On my computer, a 3 GHz Pentium 4 running Redhat linux,
  and with MATLAB version 7.2.0.294 (R2006a), the results are:
... via my 'for' loop:
Elapsed time is 16.660016 seconds.
... via the command a'*b:
Elapsed time is 0.036960 seconds.
Difference in answers is 0
... and via command dot(a,b):
Elapsed time is 0.165563 seconds.

Problem 2 (version 2): vectorization

Matrix multiplication in different ways. Using the "for" loop, there are several ways to do this. On the innermost loop, if you use a "for" loop for the third time, it probably will take minutes or hours to run, so I won't demonstrate that for a 500 x 500 matrix.

We will discuss in class the benefits of pre-allocating memory for a variable. Methods 1 and 2 below differ only in that the second method has already allocated space in RAM for the variable.

clear
randn('state',2008);
N  = 500;
B = randn(N); C = randn(N);
fprintf('On my computer, a 3 GHz Pentium 4 running Redhat linux,\n')
fprintf('  and with MATLAB version %s, the results are:\n',version);
disp('... method 1: via ''for'' loop, without pre-allocating memory:')
tic
for i = 1:N
    for j = 1:N
        A(i,j) = B(i,:) * C(:,j);
    end
end
toc

disp('... method 2: via ''for'' loop, with pre-allocating memory:')
tic
A2 = zeros(N);  % pre-allocate space for the variable
for i = 1:N
    for j = 1:N
        A2(i,j) = B(i,:) * C(:,j);
    end
end
toc
% Method 3: really slow!
% tic
% A3 = zeros(N);
% for i = 1:N
%     for j = 1:N
%         for k = 1:N
%             A3(i,j) = A3(i,j) + B(i,k) * C(k,j);
%         end
%     end
% end
% toc
disp('... method 4: using MATLAB''s builtin command:');
% Method 4: matlab's version
tic
A4 = B*C;
toc
fprintf('And the value of A(1,1) is %.6f\n',A(1,1) );

On my computer, a 3 GHz Pentium 4 running Redhat linux,
  and with MATLAB version 7.2.0.294 (R2006a), the results are:
... method 1: via 'for' loop, without pre-allocating memory:
Elapsed time is 6.156748 seconds.
... method 2: via 'for' loop, with pre-allocating memory:
Elapsed time is 5.010969 seconds.
... method 4: using MATLAB's builtin command:
Elapsed time is 0.317947 seconds.
And the value of A(1,1) is 0.734259

Problem 3

This shows how functions behave differently depending on the input. The function "diag", with with a matrix input, returns a vector, and if given a vector input, returns a matrix.

The simplest answer is B = diag(diag(A)); for example:

A = randn(5)
B = diag(diag(A))

A =

   -0.9762    0.6231    0.7027    1.9495   -0.8712
    0.3872    0.8008   -0.5929   -0.7131   -0.7103
   -1.8757   -0.3361   -0.6677    0.9442    0.0630
   -0.6952   -0.0010    1.8529    0.5051    0.0511
   -2.0533    0.0425    1.0980    0.3940   -0.2299


B =

   -0.9762         0         0         0         0
         0    0.8008         0         0         0
         0         0   -0.6677         0         0
         0         0         0    0.5051         0
         0         0         0         0   -0.2299

Problem 4: Solving linear systems (accuracy)

We find that none of the methods do a great job finding "x", and only the builtin linear solver does well at finding a vector xhat that gives a small residual.

N = 11;
A = hilb(N);
x = ones(N,1);
b = A*x;
% Solve, method 1:
xhat1 = A\b;
% xhat1 = linsolve(A,b)  % is very similar
fprintf('Error, via method 1, is %.2e\n', norm(xhat1-x) );
fprintf('Residual, via method 1, is %.2e\n\n', norm(A*xhat1-b) );

% Solve, method 2:
xhat2 = inv(A) * b;
fprintf('Error, via method 2, is %.2e\n', norm(xhat2-x) );
fprintf('Residual, via method 2, is %.2e\n\n', norm(A*xhat2-b) );

% Solve, method 3:
xhat3 = invhilb(N) * b;
fprintf('Error, via method 3, is %.2e\n', norm(xhat3-x) );
fprintf('Residual, via method 3, is %.2e\n\n', norm(A*xhat3-b) );

fprintf('Actual x \t xhat1 \t xhat2 \t xhat3\n')
[x, xhat1, xhat2, xhat3]

Error, via method 1, is 1.05e-02
Residual, via method 1, is 5.98e-16

Error, via method 2, is 2.46e-02
Residual, via method 2, is 8.18e-04

Error, via method 3, is 3.33e-02
Residual, via method 3, is 4.79e-03

Actual x 	 xhat1 	 xhat2 	 xhat3

ans =

    1.0000    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000    1.0000
    1.0000    1.0001    1.0003    1.0001
    1.0000    0.9994    1.0021    0.9976
    1.0000    1.0022    0.9993    1.0088
    1.0000    0.9951    0.9814    0.9844
    1.0000    1.0066    1.0127    1.0273
    1.0000    0.9945    1.0000    0.9961
    1.0000    1.0025    1.0098    1.0044
    1.0000    0.9995    0.9999    0.9990

Problem 5: Solving linear systems (time)

Here's what a small Hadamard matrix looks like:

hadamard(8)

ans =

     1     1     1     1     1     1     1     1
     1    -1     1    -1     1    -1     1    -1
     1     1    -1    -1     1     1    -1    -1
     1    -1    -1     1     1    -1    -1     1
     1     1     1     1    -1    -1    -1    -1
     1    -1     1    -1    -1     1    -1     1
     1     1    -1    -1    -1    -1     1     1
     1    -1    -1     1    -1     1     1    -1

This matrix is well-conditioned ( you can check with "cond( )" ). For both methods, we get no error (to within machine precision).

N = 2048;
A = hadamard(N);
x = ones(N,1);
b = A*x;
clc
% Solve, method 1:
tic;xhat = A\b;toc
fprintf('Error, via method 1, is %f\n', norm(xhat-x) );
fprintf('Residual, via method 1, is %f\n\n', norm(A*xhat-b) );

% Solve, method 2:
tic;xhat = inv(A) * b;toc
fprintf('Error, via method 2, is %f\n', norm(xhat-x) );
fprintf('Residual, via method 2, is %f\n\n', norm(A*xhat-b) );

Elapsed time is 0.900489 seconds.
Error, via method 1, is 0.000000
Residual, via method 1, is 0.000000

Elapsed time is 2.195640 seconds.
Error, via method 2, is 0.000000
Residual, via method 2, is 0.000000