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.

In particular, if A and B are both rank one, then

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

Get the MATLAB code

Published with MATLAB® R2020a