Unexpected results using SVD to separate components that make up a function
75 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
John D'Errico
2025 年 8 月 14 日 12:44
編集済み: John D'Errico
2025 年 8 月 14 日 14:29
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]
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)
Indeed, R1 is a rank 1 matrix. Another rank 1 matrix is just the matrix of all ones.
rank(ones(5))
And the sum of those two matrices would have rank no more than 2.
rank(R1 + ones(5))
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)
rank(R5)
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)
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)
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.
その他の回答 (2 件)
Chuguang Pan
2025 年 8 月 14 日 11:59
As far as I konw, the SVD of a matrix P is equal to the summation of rank-1 components.

1 件のコメント
John D'Errico
2025 年 8 月 14 日 14:34
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.
Torsten
2025 年 8 月 14 日 14:12
編集済み: Torsten
2025 年 8 月 14 日 14:37
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")
3 件のコメント
Torsten
2025 年 8 月 14 日 17:00
編集済み: Torsten
2025 年 8 月 14 日 17:22
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.
参考
カテゴリ
Help Center および File Exchange で Annotations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!