Solve large linear systems with Parallel computing toolbox

Dear all,

I need to solve something of this form , where are just 7 distinct doubles and I is the identity matrix. Of course, those linear systems can be solved in parallel, and I want to do that in Matlab with the PCT. The matrix A is A = gallery('poisson',n).

In my cluster, I have a node with 16 CPUs, and I want to use this fact to boost the performance. I wrote the following code to see if the parfor gives an improvement w.r.t the classical for. I started a parallel pool of 7 workers, and when I run it on the cluster I specified to use 7 CPU cores, according to the phylosophy "1 worker per CPU core", but my performance does not get better.

Here's the code with the following output:

clear all
close all

m = 70^2;
A = gallery('poisson',70);
I = speye(m);
v = ones(m,1);
x = zeros(m,7);
theta = [1.1,0.2,5.6,0.2,6,8,9.9];
tic
for i=1:7
	x(:,i) = (A - theta(i)*I)\v;
end
toc

parpool(7)
tic
parfor i=1:7
    x(:,i) = (A - theta(i)*I)\v;
end
toc

The results are:

Elapsed time is 0.184104 seconds. (with for loop)

Elapsed time is 0.451166 seconds. (with parfor loop)

So my questions are:

• is there something wrong in how I wrote my code to run in parallel? How can I improve my performance? (iterative solvers, or differen methods)
• why the parfor considerably slower than the classical for? I've seen that linear algebra operations are already multithreaded and hence there could be no gain with a parfor

Answer 1

[Dana]

First of all, when I run the code you posted, Matlab gives warnings that A-theta(i)*I is singular to working precision for i=2,4. Not sure if that's expected.

Second, I actually do find parfor faster when I run your code (though only barely). Not sure why you're finding otherwise, but it could be something specific to your processor.

[GiuliaC]

Thanks Dana for the answer!

That is strange. Where are you running your code? And how many physical cores does your machine has?

You're right about the theta. I mistyped the coefficients for i=2,4. Indeed, they shift the eigenvalues towards zero.

[Dana]

The computer I ran it on has 4 cores, but it's an Intel processor with hyperthreading, so effectively it has 8 "virtual" cores. I have MATLAB set up to use 7 workers when running parallel loops.

The thing with running parfor loops is that it requires a bunch of extra "overhead", i.e., extra stuff MATLAB must do just to set itself up to run the calculations you want. It's only worthwhile to run a parfor loop instead of a regular for loop if the computational savings while running the loops more than offsets the extra overhead. That cost-benefit analysis will depend on a number of things that are particular to your machine: what's worth it on one machine may not be worth it on another.

[GiuliaC]

Thanks Dana for your comment. Yes, then of course for small dimensions (and just 7 iterations), the for loop is the choiche, as it has no extra overhead.

I tried to increase m from 70 to 500, i.e. the matrix now is 500x500, and indeed I find:

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Elapsed time is 5.922038 seconds. (FOR LOOP)

Starting parallel pool (parpool) using the 'local' profile ... connected to 7 w$

Pool with properties:

Connected: true

NumWorkers: 7

Cluster: local

AttachedFiles: {}

IdleTimeout: 30 minute(s) (30 minutes remaining)

SpmdEnabled: true

Elapsed time is 2.192346 seconds. (PARFOR FOR LOOP)

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

So as the matrix size increases, it seems worthwile to use a parfor loop.

Just another question: do you have any adivise on how to solve those linaer systems? I just used MatLab backslash \ and since the matrix are saved in sparse format it works fine, but I don't know if I should consider using an iterative method.

Thanks

[Dana]

Backslash is almost always the best way to go as long as the matrix you're "inverting" (A - theta(i)*I in this case) is non-singular, so that a unique solution is guaranteed to exist. In this case, all the different methods you could use will give the same answer, but the backslash method is almost always the fastest/most accurate.

Where the different methods differ is when the matrix you're inverting is singular. In that case, depending on the precise structure of the problem, there are either no solutions, or an infinite number of solutions. When there are no solutions, each method returns an answer that solves the (unsolvable) system as closely as possible, but the different methods differ in how they define "closeness" and therefore may provide different answers. On the other hand, when there are an infinite number of solutions, each method essentially "selects" one of those solutions according to some criterion, and the different methods again use different criteria, and therefore return different answers. In either of these scenarios, you'd need to think carefully about which method would best suit your needs.

[Bruno Luong]

Just another question: do you have any adivise on how to solve those linaer systems? I just used MatLab backslash \ and since the matrix are saved in sparse format it works fine, but I don't know if I should consider using an iterative method.

It depends on A. If A is sparse or symmetric or symmetric definite positive, or hermitian etc... you might gain with iterative solver (relative cheap to evaluate A*x, and fast converging).

Note that inv(I - M) is (I +M + M^2 + M^3 ....) under some condiition. You might also able to use this if you know the eigen decomposition of M (proportinal to A).

Someone might give you better hints if you are willing tell us more about your matrix.

[GiuliaC]

Thanks Bruno for your comment.

Yes, A is sparse and symmetric, and negative definite. As you can see from the code, it's the finite difference discretization matrix of the laplacian in 2D.

The fact is that I don't know if iterative methods are suitable to run in parallel, and if this gives a speed-up w.r.t sparse direct solvers. Do you hav e any advice?

[Bruno Luong]

You should then definitively try look into using one of the iterative solvers such as pcg, cgs, and friends.

In such method you can provide "A" through a function, in your case it's come down to computing

y = (-A*x +theta_i*x) 

for any arbitrary given vector x, where A is sparse.

This must speed up, and if furthermore you could provide and cheap approximation of inv(A) for preconditioning, it will speed up even more.

I have no idea about efficiency of par-for since I do not own the parallel computing toolbox.

EDIT: Some reading for you on iterative solver and multigrid solver.