UW AMath

AMath 584, Autumn Quarter, 2011

Applied Linear Algebra and

Introductory Numerical Analysis

Table Of Contents

Matlab

Demo 1

See also the Python version of Demo 1 for comparison.

[demo1.m]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
% Illustrate a few things for a real or complex matrix...

caseno = 1   % set to 1 or 2

switch caseno
    case 1,
       A = [1,2,1; 2,4,1; 3,6,1; 4,8,1];
    case 2,
       A = [1+1j, 1; 2+2j, 2; 3+3j, 3];
    end

% print A:
A

disp('The adjoint of A is:')
A'

disp('The matrix (adjoint of A) * A is:')
A'*A

rankA = rank(A);
disp(sprintf('The rank of A is %g', rankA))

disp('Basis for the null space of A:')
null(A, 'r')     % 'r' means make rational

x = ans          % results of previous command
disp('b = A*x:')
b = A*x          % * is matrix multiplication!

disp(sprintf('2-norm of A is %g', norm(A)))
disp(sprintf('1-norm of A is %g', norm(A,1)))
disp(sprintf('max-norm of A is %g', norm(A,inf)))

Image processing with the SVD

These examples require this function for displaying matrices as images:

[showimage.m]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
function showimage(A)
    % Show the matrix A as a grayscale image.
    % Elements should lie between 0 and 1, or will be truncated to
    % this range.

    k = rank(A);

    % restrict elements to interval [0,1]:
    A = min(A, 1);
    A = max(A, 0);

    % Put in as components of RGB map:
    % All three components equal ==> gray scale
    [m,n] = size(A);
    C = zeros(m,n,3);
    C(:,:,1) = A;
    C(:,:,2) = A;
    C(:,:,3) = A;

    image(C)
    title(sprintf('Rank %i matrix', k))

First example:

[image1.m]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
A = [0 0 0 0 0 0 0 0 0;
     0 1 1 0 0 0 0 0 0;
     0 1 1 0 0 0 0 0 0;
     0 0 0 1 1 1 1 0 0;
     0 0 0 1 1 1 1 0 0;
     0 0 0 1 0 0 1 0 0];


figure(1)
showimage(A)


rankA = rank(A)
junk = input('Hit return for low rank approximations');

[U,S,V] = svd(A)

figure(2)
for k=1:rankA
    Ak = U(:,1:k) * S(1:k,1:k) * V(:,1:k)';
    showimage(Ak)
    junk = input('Hit return to continue...');
end

Hello example from Exercise 9.3:

[helloloop.m]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
%
% Sets up the hello matrix A that is
% show in Trefethen and Bau, page 68.
%
A=zeros(15,40);
for i=2:9
     A(i,2)=1;
     A(i,3)=1;
     A(i,6)=1;
     A(i,7)=1;
end;
A(5:6,4:5)=[1 1;1 1];

for i=3:10
     A(i,10)=1;
     A(i,11)=1;
end;
A(3:4,12:15)=ones(2,4);
A(6:7,12:15)=ones(2,4);
A(9:10,12:15)=ones(2,4);

for i=4:11
     A(i,18)=1;
     A(i,19)=1;
end;
A(10:11,20:23)=ones(2,4);

for i=5:12
     A(i,26)=1;
     A(i,27)=1;
end;
A(11:12,28:31)=ones(2,4);
%
for i=6:13
     A(i,34)=1;
     A(i,35)=1;
     A(i,38)=1;
     A(i,39)=1;
end;
 A(6:7,36:37)=ones(2,2);
 A(12:13,36:37)=ones(2,2);
%
%

% Show the nonzero structure of A:
figure(1)
spy(A)

B = 0.2 + 0.6*A;  % shift values away from 0, 1
                  
[U,S,V] = svd(B);

figure(2)
showimage(B)

junk = input('Hit return for rank k approximations')

for k=1:11
    Bk = U(:,1:k) * S(1:k,1:k) * V(:,1:k)';
    showimage(Bk)
    disp(sprintf('k = %i', k))
    junk = input('Hit return to continue...');
end

Recommender systems

Example shown in class

[recommend1.m]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
% Matrix with rankings:
A = [2 2 -2 -2; 2 2 -2 -2; -2 -2 2 2; -2 -2 2 2];

% add another row:
B = [A; 0 0 2 0];

disp('Matrix of ratings:')
B

[U,S,V] = svd(B);

B1 = U(:,1)*S(1,1)*V(:,1)';

disp('Rank 1 approximation:')
B1