Matlab code to compute the smallest nonzero singular value of the matrix without using SVD

Question

I want to compute the smallest nonzero singular value of the matrix A, which is defined as follows.

Let $B = rand(500, 250)$, $A = B*B^t$, where $t$ denotes the transpose of the matrix.

I found the following matlab code to compute singular values of the matrix A which is based on the Singular value decomposition of the matrix.

svds = svd(A);                             
s = min(svds);  % smallest singular value

I want to know is there any other efficient way to smallest singular value?

Thank you in advance

Answer 1

Since $A=BB^T$, it is symmetric and positive semidefinite, and its singular values are identical to its eigenvalues. Furthermore, its nonzero eigenvalues are identical to the nonzero eigenvalues of $B^TB$. To get the smallest eigenvalue of $B^TB$, use eigs(B'*B,1,'sm').

Answer 2

Summary

The answer is...use svds.

What are the singular values?

There may be some confusion over how you get the singular values. The command svd computes the singular values and the components that you don't want. The command svds only computes the singular values.

As explained here Computing pinv, the first step in computing the full singular value decomposition of a matrix $\mathbf{A}$ is to compute the eigenvalue spectrum of the product matrix $\mathbf{A}^{*}\mathbf{A}$. $$ \sigma = \sqrt{\lambda\left( \mathbf{A}^{*}\mathbf{A} \right)} $$ These are precisely the numbers you want (after they are ordered and $0$ values are culled). These values are returned by svds.

If you want to continue and compute the eigenvectors, and resolve the domain matrices, then execute svd.

Background

For background, the SVD is very powerful, and very expensive. You just want part of the information, and you hope to avoid the performance penalty. But the heart of the complexity of the SVD is the eigenvalue problem. This demands finding the roots of the characteristic polynomial. In your case, a polynomial of order 500. The task of finding roots to general polynomials is a demanding numeric problem. So by asking for the singular values, you have incurred most of the computational cost.

Caution

As an inside, make sure you understand how to handle small singular values. There is a tough issue of deciding if a small eigenvalue is a valid element of the spectrum, or a zero eigenvalue disguised as machine noise. Some discussion is Number of Singular Values and Kernel mapping from feature space to input space.

It may be reasonable to change your requirement from finding the smallest eigenvalue to setting a threshold for smallest eigenvalue.

Keep posting

As your problem and understanding evolve, keep posting and keep the discussion going.

@Rahul's answer

User @Rahul has a better solution because he skips the unneeded step of forming the product matrix. Almost certainly, eigs, svds, and svd call the same routine to find the roots of the characteristic polynomial, and in this instance the time savings may be imperceptible. Failure to recognize that we can bypass the product matrix is a critical oversight.

Answer 3

Thanks to all answers/comments, I suggest a method based on @Rahul's method (which has converted the problem from SVD to EVD).

Summary

The solution in short: s = 1/eigs(inv(B'*B),1,'lm');

Description

There exist a Power Method for calculation of the eigenvector corresponding to the largest eigenvalue of matrix $M$. It is to compute a sequence of unit vectors:

$\frac{Mv}{||Mv||}$, $\frac{M^2v}{||M^2v||}$, ..., $\frac{M^kv}{||M^kv||}$

starting from an initial random vector $v$ which almost always converges to the eigenvector with largest eigenvalue of $M$.

If Power Method requires few iterations to converge and also calculating the inverse matrix $(B^tB)^{-1}$ requires relatively lower computations, then we can use it to compute eigenvector $v$ with largest eigenvalue of $(B^tB)^{-1}$ which is the eigenvector with smallest eigenvalue of $B^tB$. finally use the inverse of largest eignvalue of $(B^tB)^{-1}$.

Having the eigenvector $v$ its corresponding eigenvalue can be obtained using $Mv=\lambda v$ (i.e: with $M = B^tB$ eigenvalue of v is equal to all elements of vector $M*v ./ v$ hence with less calculations $M(1,:)*v / v(1)$).

I have conducted a comparison of results and running times of all above methods in MATLAB.

n = 100; 
t1 = 0; t2 = 0; t3 = 0; t4 = 0; t5 = 0; % for measuring running times

B = rand(500,250); 
A = B*B'; 
for i = 1:n,
 tic; svds = svd(A); s1 = min(svds(svds > 1e-6)); t1 = t1 + toc;   % @srijan
 tic; s3 = eigs(B'*B,1,'sm'); t2 = t2 + toc;                       % @Rahul
 tic; svds = svd(B); s2 = min(svds(svds > 1e-6))^2; t3 = t3 + toc; % @Hans Engler
 tic; s4 = 1/eigs(inv(B'*B),1,'lm'); t4 = t4 + toc;                % @Alireza Kazemi(*)
end; 
% Computed smallest non-zero signular value of A
[s1,s2, s3, s4]
%-> RESULT: 3.403670094444867  3.403670094444888  3.403670094444812  3.403670094444802
% average times
[t1/n , t2/n, t3/n, t4/n]
%-> RESULT: 0.1193    0.0410    0.0317    0.0202 (seconds)
% ration of average times to the initial question of @srijan
[t1/t1, t2/t1, t3/t1, t4/t1]
%-> RESULT: 1.0000    0.3436    0.2655    0.1693

It is clear that the proposed method is fastest among compared methods, even though it needs a matrix inversion.

I must also mention that I did not implemented the Power Method here but just guessed that Matlab $eigs(...,'lm')$ must have implemented it.