Cholesky vs Eigs speed comparison?

I have many large (8000 x 8000) sparse PSD matrices for which I would like to verify that their largest eigenvalue is at most some constant C. If A denotes such a matrix, there are two ways to check this.

1. eigs(A,1) <= C
2. chol(C*speye(size(A,1)) - A) (and inspect the output flag)

Somehow, method 1 is much faster than 2 (this difference is not due to the computation time of C*speye(size(A,1)) - A) . Why is that? The general consensus is that cholesky decomposition is the fastest way of determining PSD-ness of matrices. Is chol not as fast for sparse matrices as eigs?

Answer 1

[Bruno Luong]

Your test of PSD is not right. You should do:

eigs(A,1, 'smallestreal') > smallpositivenumber % 'sr', 'sa' for older MATLAB

EIGS compute one eigen value by iterative method. It converges rapidly for

eigs(A,1)

[Lennart Sinjorgo]

I may have been unclear in my original question. Let me clarify:

I already know that my matrices are PSD, so I don't need to compute

eigs(A,1, 'smallestreal') > smallpositivenumber

I interested in the fastest way to verify that the largest eigenvalue of A is at most some constant C. To do so, my two proposed methods are valid: the largest eigenvalue of A is at most C, if and only if the matrix

B = C*speye(size(A,1)) - A

is positive(semi)definite. As such, B has a cholesky decomposition if and only if the largest eigenvalue of A is at most C. However, computing the cholesky decomposition of B is more costly than simply computing

eigs(A,1)

which suprises me. Do you have some explanation for that?

[Bruno Luong]

As I say eigs(A,1) is fast. It is no surprise to me; it is basically iterate

https://en.wikipedia.org/wiki/Power_iteration

A = rand(1000);

niter = 20;
tic
x = rand(size(A,2),1);
for i = 1:niter
   lambda = norm(x);
    x = x/lambda;
    x = A*x;
end
lambda = norm(x)

lambda = 499.7373

toc

Elapsed time is 0.088616 seconds.

few time

You coulld also use norm

tic, norm(A), toc

ans = 499.9048

Elapsed time is 0.079124 seconds.

tic, eigs(A,1), toc

ans = 499.7373

Elapsed time is 0.081136 seconds.

[Matt J]

You coulld also use norm

norm(A,2) doesn't work for sparse matrices, so you do have to use eigs(A,1). However, I just wanted to mention that norm(A,inf) is faster than eigs(A,1) and provides an upper bound on it. So, it might be worth checking if norm(A,inf)<C first.

A=sprand(8000,8000,1/8000*50);

tic;
 Ninf=norm(A,inf);
toc

Elapsed time is 0.008915 seconds.

tic;
N2=eigs(A,1);
toc

Elapsed time is 0.111935 seconds.

Ninf>N2

ans = logical
   1

[Bruno Luong]

norm(A,inf) can get quite larger. It depends on what is the goal it can be a problem or not

n = 2^13;
[A,~] = qr(randn(n));
norm(A,inf) / norm(A,2)

ans = 72.9242

We know that both norms are not equivament when the dimensions of the matrix change freely.

[Bruno Luong]

An example with sparse matrix based on Haar basis. The ratio of two matrix norms increases like sqrt(N)

N = 2^14

N = 16384

A = haar(N);
% eigs(A,1) must be 1
norm(A,inf) / eigs(A,1)

ans = 128.0000

function S=haar(N)
if (N<2 || (log2(N)-floor(log2(N)))~=0)
    error('The input argument should be of form 2^k');
end
p=[0 0];
q=[0 1];
n=nextpow2(N);
for i=1:n-1
    p=[p i*ones(1,2^i)];
    t=1:(2^i);
    q=[q t];
end

j = 1:N;
I = 1+0*j;
J = j;
V = 1+0*j;
for i=2:N
    P=p(1,i); Q=q(1,i);
    j = 1+N*(Q-1)/(2^P):N*(Q-0.5)/(2^P);
    I = [I, i+0*j];
    J = [J, j];
    V = [V, 2^(P/2)+0*j];
    j = 1+(N*((Q-0.5)/(2^P))):N*(Q/(2^P));
    I = [I, i+0*j];
    J = [J, j];
    V = [V, -(2^(P/2))+0*j];
end
S = sparse(I,J,V,N,N);
S=S/sqrt(N);
end
        

[Bruno Luong]

Here is an example where the ratio is huge (= N), discovered by codinng mistake .;)

N = 2^16

N = 65536

j = 1:N;
I = 1+0*j;
J = j;
V = 1+0*j;
A = sparse(I,J,V,N,N);
norm(A,inf) / eigs(A,1)       

ans = 6.5536e+04