Mount St. Helens and Matrix Rank
What do Mount St. Helens and the rank of a matrix have in common? The answer is the MATLAB function peaks. Let me explain. Please bear with me -- it's a long story.
Contents
Mount St. Helens
Mount St. Helens is an active volcano in the Cascade Range near Portland, Oregon and Seattle, Washington. She had a dramatic eruption in May 1980. Four years after the eruption, I moved to Beaverton, Oregon to join a spinoff from Intel that was developing one of the world's first commercial parallel computers, the iPSC. In Beaverton, I met two guys from Tektronix, a company famous for oscilloscopes that was then trying to get into computer graphics. These guys, whose names I don't remember, were plotting a surface that resembled Mount St Helens. That surface eventually became the MATLAB function peaks, which is a function of two variables that is useful for demonstrating surf, mesh, pcolor and contour.
small_fig
clear
P = peaks;
surf(P)
axis tight
Rank Revealing Factorizations
Two months ago, I wrote a blog post about rank revealing matrix factorizations, It is hard to think of the array P that we obtain from peaks as the representation of a linear transformation and a test matrix, but it is. In fact, for any positive integer n, peaks(n) is an n -by- n matrix. The default has n = 49.
I was very surprised when I discovered that for any n > 2, the rank of peaks(n) is only 3.
rankP = rank(P)
rankP = 3
C*inv(W)*R
MATLAB's rank counts nonzero singular values, but we can also try the factorization P = C*inv(W)*R. We intend to use car only on small matrices that have integer entries. This P is neither small nor integer, but that doesn't stop us. car also finds the rank is 3.
[C,W,R] = car(P); whos
Name Size Bytes Class Attributes C 49x3 1176 double P 49x49 19208 double R 3x49 1176 double W 3x3 72 double rankP 1x1 8 double
The 49-by-3 matrix C is the first three columns of P and the 3-by-49 matrix R is the first three rows of P. The 3-by-3 matrix W is their intersection, P(1:3,1:3). These three matrices lie safely out of the reach of volcanic eruptions. Here they are, in color.
d = .1; n = 49; r = 3; s = [patch([1-d n+d n+d 1-d 1-d],[1-d 1-d r+d r+d 1-d],[d d d d d],'c') patch([1-d r+d r+d 1-d 1-d],[1-d 1-d n+d n+d 1-d],[d d d d d],'y') patch([1-d r+d r+d 1-d 1-d],[1-d 1-d r+d r+d 1-d],[d d d d d],'m')]; legend(s,{'C','R',' W'}, 'position',[.75 .62 .12 .12]);
Plots of C and R do not easily reveal the fact that they contain enough information to reconstruct the entire surface P.
subplot(2,1,1) plot(C,'.-','linewidth',1,'markersize',12) title('C') subplot(2,1,2) plot(R','.-','linewidth',1,'markersize',12) title('R')
peaks
The function peaks begins by generating matrices x and y that sample a continuous real variable on the interval [-3,3]. Setting the meshsize dx = 1/8 leads to n = 49. Smaller mesh sizes lead to larger matrices, but they all have rank 3.
dx = 1/8; [x,y] = meshgrid(-3:dx:3);
The core of peaks is this formula, which must be something like the one that the Tektronix guys were using for their Mt. St. Helens It is not obvious that the formula always generates a matrix of low rank.
z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2);
Rank
The rest of this post is devoted to proving that rank(z) = 3.
Here are two important facts about the rank of sums and pointwise products.
- rank(A + B) <= rank(A) + rank(B)
- rank(A.*B) <= min(rank(A),rank(B))
In particular, if A and B are both rank one, then
- rank(A + B) <= 2
- rank(A.*B) <= 1
Reformulate
Let's factor out the exponentials using this function.
E = @(x,a) exp(-(x+a).^2);
So, here is z again.
z = 3*(1-x).^2.*E(x,0).*E(y,1) ... - 10*(x/5 - x.^3 - y.^5).*E(x,0).*E(y,0) ... - 1/3.*E(x,1).*E(y,0);
Let z1 be the first term.
z1 = (3*(1-x).^2).*E(x,0).*E(y,1);
rank(z1) = 1 because z1 is the pointwise product of three matrices. There is no sum.
Let
zx = -10*(x/5 - x.^3).*E(x,0) - 1/3.*E(x,1); zy = (10*y.^5).*E(x,0);
rank(zx) = 1 because zx does not involve y.
rank(zy) = 1 because zy is a pointwise product.
z2 = (zx + zy).*E(y,0);
rank(z2) = (1 + 1)*1 = 2
z = z1 + z2;
rank(z) = 1 + 2 = 3
Check
My original P is still in the workspace.
check = norm(P - z,inf)
check = 1.3517e-14
评论
要发表评论,请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。