Krylov Subspace Methods for Dense Systems

Question

I am currently researching on the viability of using KS methods for solving large dense systems. What I wish to prove (or disprove) is that methods like CG, BiCG and QMR are as good (if not better) than the generic LU or QR decomposition methods in place today.

Till now I have tried non-preconditioned versions of CG on symmetric positive definite matrix (duh!) and compared it against LU(DGESV). (I have written so in Fortran and MATLAB as well) The results are not pleasing. CG loses by orders of magnitude.

I have written (or found and optimized) codes for QMR, BiCG, CGNE, CGNR and others as well but their results are worse.

Now, I have two choices

1. Try improve BLAS Level 2 routines to yield a better convergence. ( I use Intel MKL in Fortran). But, it turns out, there is nothing I can do to make this better. Or can I? Profiling shows bad use of threads and matrix-vector performance. Which is not surprising for BLAS 1/2 operations.

2. Work on preconditioners and hope that I come across a good one which solves a linear system faster than LU/QR etc.

3. (Bonus Option) Leave this idea and try something else for my thesis. (I have time left to do this). Maybe something else in the domain of Krylov Subspace Methods.

Any help is much appreciated. P.S. I had no clue this group existed. This is awesome!

Answer 1

I am extremely surprised that there is no mention of conditioning or the shape of the spectrum in your discussion, as it will be the decisive property in whether or not iterative methods can beat dense methods.

As an extreme example, suppose that your dense matrix is some small perturbation of the identity matrix. Then most iterative methods will converge in $O(1)$ iterations, and dense matrix-vector multiplication is $O(n^2)$, so the solve cost is $O(n^2)$ rather than $O(n^3)$. There is clearly a large class of well-conditioned dense matrices that similar arguments could be made for. Whether or not they show up in practice is another question...

In my experience, when I generate random dense matrices, the performance of GMRES (and CG for HPD matrices) is rather terrible relative to an LU or QR solver. However, iterative methods would certainly be faster for very large matrices when only a handful of iterations are required.

EDIT: Since I was asked about experimental results, I thought I would show some MATLAB code and some timings which make the point a little more clear. The following generates a random $2000 \times 2000$ real dense SPD matrix with eigenvalues uniformly sampled from [1,2] and compares the timings of unpreconditioned CG (to 6 digits of accuracy) and a Cholesky solver. On my two year old iMac, CG is already four times faster.

    %
    % Compares CG and a dense Cholesky solve for a well-conditioned SPD matrix
    %
    n=2000;
    S=1;
    T=2;
    tol=1e-6;

    % Build an n x n dense SPD matrix with eigenvalues uniformly sampled from 
    % [S,T]
    D=diag(random('uniform',S,T,n,1));
    U=randn(n,n);
    [U,R]=qr(U);
    clear R;
    A=U*D*U';

    % Set up a random right-hand side
    b=randn(n,1);

    % Solve with CG to the specified tolerance, ensuring that the maximum
    % number of iterations is not reached
    fprintf(1,'Running CG...\n');
    tic; [x,flag,relres,iter]=pcg(A,b,tol,10*n); toc;
    iter

    % Solve using a Cholesky decomposition
    fprintf(1,'Running a Cholesky solver...\n');
    opts.LT=false; opts.UT=false; opts.UHESS=false;
    opts.SYM=true; opts.POSDEF=true;
    opts.RECT=false; opts.TRANSA=false;
    tic; x=linsolve(A,b,opts); toc;

I saw the following results:

Running CG... Elapsed time is 0.074933 seconds.

iter =

 8

Running a Cholesky solver... Elapsed time is 0.279619 seconds.

Now, if $A$ is instead generated as the symmetric product a matrix with normally distributed entries, e.g.,

    A=randn(n,n);
    A=A'*A;

then the situation is not so good for CG. I saw the following output:

Running CG... Elapsed time is 28.568865 seconds.

iter =

    3892

and the condition number was 1.0587e+09.

Answer 2

In most cases I know of where Krylov methods are used for dense problems, the operator is a low-rank perturbation of the identity (obtained by discretizing a continuum operator which is a compact perturbation of the identity). Such operators appear frequently in boundary integral equations as discretized Fredholm integral operators of the second kind. These operators have a few extreme eigenvalues, but the spectrum decays rapidly to the identity, so Krylov methods can be very effective even though the system is ill-conditioned due to the extreme eigenvalue.

For these problems, the operator might be applied using some fast method such as FMM, but are sometimes assembled using dense summation (especially for 2D problems). If the boundary and/or coefficient structure is very complicated, the decay of the spectrum might not be particularly fast, in which case you might consider using "multigrid of the second kind".

In other examples, the "dense" operator might be obtained as a Schur complement, in which case there is likely a fast way to apply it without assembling, and approximations that are useful for preconditioning might be obtained through approximate commutator arguments.

As a general rule, Krylov methods are not good general-purpose methods for dense systems, but if your problem has exploitable structure such as good conditioning, decay of the spectrum, fast application, or effective preconditioners, then Krylov methods can certainly be useful.

Answer 3

CG was originally abandoned because of this loss of accuracy that you observed. It was revived when its use as an approximate solver was appreciated. So I don't think you want to revisit this problem. And preconditioning is very useful, diagonal preconditioning at the very least.

Iterative methods are very problem dependent so I would recommend finding an application that needs to solve dense systems, or even better move to sparse systems, to focus your research. N.B., Krylov methods have been extensively investigated in the past 30-40 years, it might be hard to find a good topic here unless you have a lead of some sort or an expert adviser. There are often properties of your application that can be exploited with preconditioning, which can provide the kernel of a dissertation.