Unexpected results using SVD to separate components that make up a function

As preliminary to a larger project, I was trying out the SVD function with a test matrix made up by adding three 2D gaussians p1, p2, p3. I was hoping the SVD would be able to separate the function into these three components, and they would have the largest singular values in the S matrix.

I then tried to remake these components by setting everything except one singuar value in the S matrix to zero and calculating U*S_new*V'.

This does not make anything resembling my original components. Have I misunderstood how SVD should work, or have I misunderstood how separating the components should work?

I have attached my code:

% SVD test

clear all; 

[X,Y] = meshgrid(-2:.1:4);

mu1 = [1 2];
Sigma1 = [5 3;3 5];
p1 = mvnpdf([X(:) Y(:)],mu1,Sigma1);
p1 = reshape(p1,size(X));

mu2 = [0 1];
Sigma2 = [3 0;0 1];
p2 = mvnpdf([X(:) Y(:)],mu2,Sigma2);
p2 = reshape(p2,size(X));

mu3 = [3 0];
Sigma3 = [2 1;1 2];
p3 = mvnpdf([X(:) Y(:)],mu3,Sigma3);
p3 = reshape(p3,size(X));

p = p1+p2+p3;

[U, S, V] = svd(p, "econ");

p_remade = U*S*V';

figure; 
subplot(3,1,1); imagesc(p1)
subplot(3,1,2); imagesc(p2)
subplot(3,1,3); imagesc(p3)
sgtitle("original components making p")

figure; 
subplot(2,1,1); imagesc(p)
title("original p")
subplot(2,1,2); imagesc(p_remade)
title("SVD p")

num_comp = 3;
Si = zeros(size(S));
figure;
comp = zeros(length(S),length(S),num_comp);
for i = 1:num_comp
    Si(i,i) = S(i,i);
    comp(:,:,i) = U*Si*V';
    subplot(num_comp,1,i); imagesc(squeeze(comp(:,:,i)));
end

Answer 1

[John D'Errico]

Yes. You do misunderstand how svd works, and what it produces. As well, I think you misunderstand what a rank 1 matrix means.

A rank 1 matrix means that every row of a matrix is a simple multiple of a given vector. The same applies to every column. The following is a rank 1 matrix:

R1 = [1:5]'.*[1:5]

R1 = 5×5

1 2 3 4 5 2 4 6 8 10 3 6 9 12 15 4 8 12 16 20 5 10 15 20 25

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Do you see that every row or column is a scalar multiple of the other rows? Now, test it out. Is this a rank 1 matrix?

rank(R1)

ans = 1

Indeed, R1 is a rank 1 matrix. Another rank 1 matrix is just the matrix of all ones.

rank(ones(5))

ans = 1

And the sum of those two matrices would have rank no more than 2.

rank(R1 + ones(5))

ans = 2

But even something as simple as a diagonal matrix is no longer rank 1. It will be full rand, so here, rank 5, since you CANNOT write any column or row as a linear combination of the other rows or columns.

R5 = eye(5)

R5 = 5×5

1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rank(R5)

ans = 5

Now consider mvnpdf. Interestingly, a multivariate normal, with a purely diagonal covariance matrix also has the property that it would create a rank 1 matrix, because in that case, the multivariate pdf is separable. You can split it into the product of two univariate PDFs for the bivariate case. This claim is simple enough to test (the algebra is pretty simple too, but I'm feeling too lazy to do that for now. Ask, if you don't trust me. Hey you can buy a used car from me, really. Seriously, ask if it seems important and my numerical example is not sufficient, and I can show how that works in a comment.)

[x,y] = meshgrid(linspace(-3,3,25));
z = reshape(mvnpdf([x(:),y(:)]),size(x));
rank(z)

ans = 1

And indeed the bivariate pdf with no off-diagonal compoents in the covariance matrix results in a rank 1 matrix.

In your case however, you made one big initial mistake. The covariance matrix you chose will mathematically result in a full rank matrix. We can show the problem here:

mu1 = [1 2];
Sigma1 = [5 3;3 5];
z1 = reshape(mvnpdf([x(:),y(:)],mu1,Sigma1),size(x));
rank(z1)

ans = 13

As it turns out, numerically, it had only rank 13. But that was due to floating point issues in double precision. Mathematically, it is actually full rank. You CANNOT write any column of that matrix as a linear combination of the rest of the columns. However, in floating point arithmetic, z1 is not truly numerically full rank.It is still rank 13 though, so it has an intermediate computed rank.

Next, the sum of two such PDFs, even if they were each rank 1, will not allow you to recover the original PDFs in a way that will look like a pair of PDFs. This is because the SVD produces singular vectors that are orthogonal. And that means the two vectors you would get out, even if the component matrices in the sum were rank 1, will not look like a normal PDF! This was your second major mistake.

In the end, for your project, I'm not certain what it is that you wanted to do in your project. But an SVD may not be the approriate tool to solve your problem. Since we are not told what problem you really wanted to solve and what is your project, it is difficult to know. I would strongly recommend you sit down with your advisor/teacher and discuss this with them.

[Nasim]

Thank you, this makes a lot of sense.

Answer 2

[Chuguang Pan]

As far as I konw, the SVD of a matrix P is equal to the summation of rank-1 components.

Since the matrix is not rank-1, I think svd function can not seperate them.

[John D'Errico]

Yes, you are correct. The SVD will essentially split any matrix into linear combinations of rank 1 terms.

But even at that, the components it will generate will be orthogonal. As such, it would not reproduce the original PDFs, since they were not orthogonal themselves. The result would be arbitrarily random looking, and would not resemble a PDF at all.

Answer 3

[Torsten]

Others have already explained why you cannot recover p1, p2 and p3 from an SVD of p = p1 + p2 + p3.

All you can do is write

p = sum_{i=1}^{i=61} S(i,i)*U(:,i).*(V(:,i).')

and it will turn out that the quality of approximation of p is already remarkable if you don't sum up to 61, but only up to 3. But note that the number 3 has nothing to do with the fact that p is the sum of 3 different PDFs - it's only based on the magnitude of the first 3 singular values.

% SVD test
clear all; 
[X,Y] = meshgrid(-2:.1:4);
mu1 = [1 2];
Sigma1 = [5 3;3 5];
p1 = mvnpdf([X(:) Y(:)],mu1,Sigma1);
p1 = reshape(p1,size(X));
mu2 = [0 1];
Sigma2 = [3 0;0 1];
p2 = mvnpdf([X(:) Y(:)],mu2,Sigma2);
p2 = reshape(p2,size(X));
mu3 = [3 0];
Sigma3 = [2 1;1 2];
p3 = mvnpdf([X(:) Y(:)],mu3,Sigma3);
p3 = reshape(p3,size(X));
p = p1+p2+p3;
[U, S, V] = svd(p, "econ");
p_remade = U*S*V';
% Use only first three singular values to approximate p
num_comp = 3;
p_approx = zeros(size(p));
for i = 1:num_comp
    p_approx = p_approx + S(i,i)*U(:,i).*(V(:,i).');
end
figure; 
subplot(3,1,1); imagesc(p1)
subplot(3,1,2); imagesc(p2)
subplot(3,1,3); imagesc(p3)
sgtitle("original components making p")

figure; 
subplot(3,1,1); imagesc(p); colorbar
title("original p")
subplot(3,1,2); imagesc(p_remade) ; colorbar
title("SVD p")
subplot(3,1,3); imagesc(p_approx);colorbar
title("approximate p")

[John D'Errico]

@Torsten

You seem to be missing the point though. Suppose you start with a single PDF. Can you then reconstruct that PDF from only the first component using an SVD? And the answer is: Only very poorly, if that multivariate PDF has non-zero covariances.

% SVD test
clear all; 
[X,Y] = meshgrid(-2:.1:4);
mu1 = [1 2];
Sigma1 = [5 3;3 5];
p1 = mvnpdf([X(:) Y(:)],mu1,Sigma1);
p1 = reshape(p1,size(X));
p = p1;
[U, S, V] = svd(p, "econ");
p_remade = U*S*V';
% Use only first three singular values to approximate p
num_comp = 1;
p_approx = zeros(size(p));
for i = 1:num_comp
    p_approx = p_approx + S(i,i)*U(:,i).*(V(:,i).');
end
figure; 
imagesc(p1);colorbar
sgtitle("original p1 component")

figure; 
imagesc(p_approx);colorbar
title("approximate p")

With only one component available, the result will be only a rank 1 matrix. And therefore, we will see now a circular PDF, instead of the elliptical one of the original. Again, this ends up as a poor approximation. That p1 array is actually a mathematically rank 61 matrix (though it will be far lower when computed in double precision. Probably only numerically rank 15 or so.)

Using 3 components when you had 3 terms in the sum is not the same thing as extracting them independently from the sum.

[Torsten]

@John D'Errico

You showed what cannot be done using SVD, I showed what can be done. Of course you are right: it doesn't help for the OP's task. But I think it's also worth to mention what SVD is good for in a different context. Especially I refer to the last part of the OP's code

num_comp = 3;
Si = zeros(size(S));
figure;
comp = zeros(length(S),length(S),num_comp);
for i = 1:num_comp
    Si(i,i) = S(i,i);
    comp(:,:,i) = U*Si*V';
    subplot(num_comp,1,i); imagesc(squeeze(comp(:,:,i)));
end

corrected it and placed it in the correct context.

Do you know of another mathematical method that can do what the OP aims at ? My guess is that it's not possible to decompose the signal p = p1 + p2 + p3 in an adequate way.

[Nasim]

Thank you for the discussion. Luckily, the problem is only that I chose a poor example to test my functions, not that I need an alternate solution for my actual project. We (and many others in the field) have been using SVDs to remove noise in our data and reduce it to the main components for ages. But my data is of the shape (8,600,400,100) and I have to reshape it to do SVD (8x600x400, 100) - and eventhough I can make it work in python, in Matlab I was not getting the correct first few components. I thought making up simpler data would make debugging easier, and at least I got to clear up a few misunderstandings :)