solve efficiently sparse large linear system of equations

Hello Community,

i have a large sparse block structured matrix.

with this quadratic sparse Matrix L:

is there any faster way to get a solution.

here is my code: (a tol of 1e-6 works the best for my problem)

opt.setup = milu;
[A,B]=ilu(L,opt.setup);
[sol,flag]=gmres(L,f,[],1e-6,size(L,1),A,B);

Best Regards,

Marko

Answer 1

[Fabio Freschi]

I would let the backslash operator work on this matrix

sol = L\f;

Can you share the original matrix to make other attempts?

[Fabio Freschi]

My understanding is that this is a small example of a larger problems. I deal with such systems when I use hybrid formulations coupling, for example, FEM and BEM. In these cases I try to exploit the structure of the block-matrix. In fact it seems that some blocks have a nice structure

• block 1,2 seems diagonal
• block 2,3 and 3,2 seem block diagonal
• the first ald last rows of block 1,1 seem a diagonal matrix (dirichlet boundary conditions?)

in these case it might help unfolding the calculation. If you can share more details about the matrix (and the matrix itself), I can give it a try.

[Marko]

Hello Fabio,

i use the chebyshev collocation method for solving the navier stokes equation. And because the equations are nonlinear, i have to solve them may times iteratively.

Some lines of the matrix are for the dirichlet b.c.

I have uploaded L and f as mat-file.

L\f is really expensive for solving (time and memory)

[John D'Errico]

So, no. Apparently it is NOT a problem where L in the real problem is much larger.

@Marko, you have not answered my questions. Does L vary? This matrix is relatively tiny. If it does not vary, then you shoild just factorize it ONCE.

Anyway, an important point is this matrix is singular.

>> cond(full(L))
ans =
      5.41307201046228e+18

[Fabio Freschi]

If the matrix comes from the numeric solution of navier stokes equations, I guess the OP has a different matrix any times he changes anything in the problem (geometry, mesh, material parameters), so my understanding is that this matrix may be also larger in more complex cases.

The singularity of the matrix is an important issue, especially if the singularity arises from a mistake in the formulation or a wrong assignment of the boundary conditions. I suggest @Marko to investigate this problem. However, there are cases of numerical methods where the matrices are singular if they are not gauged. And gauging may not be necessary when iterative solvers are used.

A short OT

I don't really like imposing dirichlet boundary conditions by zeroing the row and putting one in the diagonal. I prefer removing the fixed variables from the system, keeping the original properties of the matrix

clear variables, close all

% dummy random matrix
n = 10;
A = rand(n)+2*eye(n);
b = rand(n,1);

% fix 3 (random) unknowns to 1 (e.g. 1 5 7)
idFix = [1 5 7];
xFix = [1 1 1];

% move the fixed unknowns to the rhs
b = b-A(:,idFix)*xFix(:);

% create the vector of true unknowns
idUnk = setdiff((1:size(b,1)).',idFix(:));

% remove unused rows and cols from matrix and rhs
A = A(:,idUnk);
A = A(idUnk,:);
b = b(idUnk,:);

% preallocate unknown vector
x = zeros(size(b));

% solve your system
x(idUnk) = A\b;

% add the fixed unknowns
x(idFix) = xFix;

Working in this way could have the advantage of preserving some features of the matrix. For example block 1,2 is not diagonal because some entries are missing (I suppose due to the zeroing of the rows).

[Fabio Freschi]

I made some tests to with your matrix and rhs.

clear variables, close all

load('Lf.mat');
figure,spy(L)

% gmres w/o preconditioned
tic
[x0,flag0,relres0,iter0,resvec0] = gmres(L,f,[],1e-6,size(L,1));
toc

Elapsed time is 11.421778 seconds.

% gmres w/ preconditioning
tic
[M1,M2] = ilu(L);
[x1,flag1,relres1,iter1,resvec1] = gmres(L,f,[],1e-6,size(L,1),M1,M2);
toc

Elapsed time is 0.133606 seconds.

% convergence plot
figure
semilogy(0:length(resvec0)-1,resvec0/norm(f),'-o')
hold on
semilogy(0:length(resvec1)-1,resvec1/norm(M2\(M1\f)),'-o')

% backslash
tic
x2 = L\f;
toc

Elapsed time is 0.810359 seconds.

% error
norm(L*x2-f)

ans = 3.0881e-12

% LU factorization
tic
[L0,U0,ip,iq,D0] = lu(L,'vector');
x3(iq) = U0\(L0\(D0(ip,ip)\f(ip)));
toc

Elapsed time is 0.602693 seconds.

% error
norm(L*x3-f)

ans = 1.5603e-10

It seems that the sparse solver does not worry to much about the singularity (but, again, this issue must be investigated). The LU factorization looks fast enough to me, also because if the problem is nonlinear, you can reuse the factors L0 and U0 and change only the rhs f. Note that most of the computation time is for the LU factorization and not in the forward/backward solution.

% LU factorization
tic
[L0,U0,ip,iq,D0] = lu(L,'vector');
for i = 1:100
   x3(iq) = U0\(L0\(D0(ip,ip)\f(ip)));
end
toc

Elapsed time is 1.346799 seconds.

As you can see you can solve the same system 100 times with a little more effort.

In case you can produce a block matrix without zeroing the rows, it could be possible to try to exploit the block-matrix nature of the matrix. See for example section 2.5 of the preprint of the paper attached. A part from the formulation details, that are not of interest here, in that case I have a 2x2 block matrix with a diagonal block (like yours) and I exploit the fact by calculating the Schur complement, and providing the result to GMRES.

Your case is 3x3 block matrix but you can give a try to check if you have benefits.

By the way: sometimes the Schur complement improves the conditioning of the matrix, with a resulting speedup in the iteration of GMRES.

Answer 2

[John D'Errico]

The fact is, a 3kx3k system is not that slow to solve. So unless you are doing this solve many thousands of times, keeping it sparse may well slow you down.

It is not even really that sparse, since for a 3100x3100 matrix with 250000 nonzeros, there are on average

249155/3100

ans = 80.3726

So roughly 80 non-zeros per row. That is not very sparse.

Have you tried a direct solve in sparse form, using a column ordering to reduce fill-in?

Have you compared the solve time if the matrix were a full matrix?

For example...

A = randn(3100);
b = randn(3100,1);
timeit(@() A\b)
ans =
           0.1878352492155

That test was run on my computer. I was watching the actuivity monitor on my computer as it ran, and the test itself was so fast that MATLAB never even fired up multiple cores in the solve, at least that the activity monitor showed in the time slice it was watching.

Even for for a full matrix of that size, it only took a small fraction of a second to solve as a full matrix.

Are you truly needing to solve this probem many thousands of times? If not, then this entire question seems moot.

If so, then is your matrix fixed, so just different right hand sides? If so, then a matrix factorization would be right, or if you have many right hands sides, known in advance, then a direct call to backslash is absolutely correct.

Or is your real problem orders of magnitude larger?

The point of my answer is this really is a small probem, and sparse form may not even be a good choice. The amount of time tken by the incomplete LU there is often a significant fraction of the time it would take to perform the matrix factorization itself.

So if you want help, then it does help if you explain all of these things, because your question as it is does not completely make sense. Why do you think you need to solve this problem as you wish to solve it?

Edit:

I see that you have provided the matrix L. Must I point out that L is singular? You claim that \ is slow here. The probem is, that backslash, when applied to a FULL matrix of that size, is actually roughly 6 times faster for the full matrix than it is for the sparse matrix.

To remove the issue of your matrix also being singular, I'll perturb the diagonal by a small amount.

Lhat = L + speye(3139)*1e-7;
rank(full(Lhat))
ans =
        3139
        
timeit(@() Lhat\f)
ans =
           1.3129801282155
           
Lhatfull = full(Lhat);
timeit(@() Lhatfull\f)
ans =
           0.2029861762155

So if you just store the matrix in a full form, \ is significantly faster for the full matrix! More than 6x faster.