# Krylov Subspace Methods for Dense Systems

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!

• What do you consider "large"? What motivated you to tackle this research problem? Jan 8 '12 at 11:04
• Large is 100000 or possibly more if we get a RAM upgrade. Our real problem is dense, non symmetric and non positive definite. We are using CG as guinea pig since it is the best iterative algorithm around. If CG can't solve faster for a Symmetric PD matrix over the generic LU, nothing can (We believe that). The motivation for the project was simply "Why Not". We were discussing the applications of KS methods and realized that apart from the inability to solve Matrix Vector swiftly, there was no real reason to hold back. We decided to try. In hindsight, I'm not sure we had any other motivation. Jan 8 '12 at 11:05
• Only some random thoughts: (i) Have you tried out and compared other CG-codes? I believe there are pretty much of them out there. (ii) Still, you can profile how CG compares to LU/QR for varying size and density ( and sparsity patterns ). Depending on how well CG performs for moderate density, you could judge whether there is hope. Jan 8 '12 at 13:31
• (i)I wouldn't be exaggerating if I said I have tried every known implementation of CG in literature (including places like Mathworks Exchange) but it doesn't help. (ii) I have tried LU/QR vs CG on a for loop running from 3 to 80,000 with an increment of 5! Near the higher values, the disparity between CG and LU/QR is humongous! LU is a 100 times faster. I haven't tried it for sparsity patterns because our application itself would never be sparse. Not the least. What exactly do you mean by density? Thanks a lot for your comments Jan 8 '12 at 13:44
• I mean, let us start with a diagonal matrix. Then a three band matrix, a five band matrix, then broader bandwidth. CG will beat LU on a diagonal matrix by orders of magnitude, but how does this develop for more and more density, where roughly is the turning point? Now you have a two-parameter space to profile your solvers. PS: Ok, as Paulson added, the spectrum of the matrix is generally more important than density. One more parameter for your profiling. Jan 8 '12 at 14:17

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.

• Just out of curiosity, Jack, how bad are the condition numbers on the random matrices where LU outperforms CG (or the like)? Also, what's the stopping criteria for comparing the two (i.e., does CG have to take the answer down to 10 times machine precision, or do you let it have more slack?)? It seems like the convergence criteria would have a huge impact on the comparison. Jan 8 '12 at 15:24
• Jack, thanks a lot for your reply. I have Deeply investigated conditioning. In fact, all my results pass through 3 filters: one of random matrices with modest conditioning, one with vandermonde type conditioning and third with hilbert type conditioning. All three have more or less similar results: CG loses to LU/QR Jan 8 '12 at 16:19
• Bill, I added some more results that show where it does and does not work. Essentially, generating a matrix with normally distributed entries and then squaring it, in the symmetric outer product sense, leads to very poor conditioning. More carefully crafting the spectrum to lie in a nice interval, like [1,2], results in very fast convergence. Jan 8 '12 at 17:34
• Thanks, Jack. Interesting results. Presumably the difference between the iterative and direct solves is order 1e-6? I wonder if MATLAB is calling the level 2 BLAS for the MATVECs. I'm pretty sure it calls the level 3 BLAS for the operations in Cholesky, but I don't know what it does for CG. Jan 8 '12 at 18:47
• Bill: Roughly, the direct solve had a residual $\|b-Ax\|_2 \approx$ 1e-13, whereas it was just over 1e-5 for the iterative solve, so the difference was almost eight digits. Jan 8 '12 at 18:53

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.

• I suspect that he meant unstructured dense matrices. I think that the standard terminology needs to change, and that the typical usage of "dense linear algebra" should be replaced with "unstructured dense linear algebra". Jan 10 '12 at 19:22

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.