decomposition used in parfor

In the following problem, where I need to solve many times same linear system "A" with different right sides "b" is very benefitial to use decomposition, which in this case works typically much faster, than simple A\b.

See this example:

clear all
N = 300;
M = 2000;
A = rand(N);      
b = rand(N,M);  
x = zeros(N,M); 

tic;
for ii=1:M
    x(:,ii) = A \ b(:,ii);
end
toc

tic;
dA = decomposition(A); 
dAc = parallel.pool.Constant(dA);
parfor ii=1:M
    x(:,ii) = dAc.Value \ b(:,ii);
end
toc

As I learned from here, TMW expert proposed to use parallel.pool.Constant, but this does not work at all!!!

Warning: Saving a decomposition is not supported. 

Is there any way how to use decomposition together with parfor or any other parallel evaluation construct???

[Bruno Luong]

Odd your test script works just fine on my PC, with R2022a Update4

Elapsed time is 1.290296 seconds.
Elapsed time is 0.148252 seconds.

[Matt J]

with R2022a Update4

Not on mine.

[Bruno Luong]

My bad, I have the Automatic creation parallel pool option turned off, so the parfor runs like like a for.

[Edric Ellis]

It should be possible to use decomposition with parallel.pool.Constant (even if it turns out you don't really need it for this case). Could you say more about what doesn't work there for you?

[Bruno Luong]

@Edric Ellis this is error we get

N = 10;
M = 4;
A = rand(N);      
b = rand(N,M);  
x = zeros(N,M); 

dA = decomposition(A); 
dAc = parallel.pool.Constant(dA);
parfor ii=1:M
    x(:,ii) = dAc.Value \ b(:,ii);
end

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 2).

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Warning: Saving a decomposition is not supported.

Error using \
Matrix dimensions must agree.

[Edric Ellis]

Ah OK, I should have been clearer. This case needs the function_handle constructor for parallel.pool.Constant.

N = 10;
M = 4;
A = rand(N);      
b = rand(N,M);  
x = zeros(N,M); 
dAc = parallel.pool.Constant(@() decomposition(A));
parfor ii=1:M
    x(:,ii) = dAc.Value \ b(:,ii);
end
disp('done.')

done.

This causes each worker to build the decomposition for itself from the value of A.

[Bruno Luong]

So this work around makes the same decomposition performed by all worker instead of once, right?

Answer 1

[Matt J]

Your example doesn't describe a situation where either decomposition() or parfor will be beneficial. You should instead just use x=A\b.

N = 300;
M = 2000;
A = rand(N);      
b = rand(N,M);  
x = zeros(N,M); 


tic;
x=A\b; 
toc%Elapsed time is 0.010876 seconds.


[L,U,P]=lu(A);
tic;
parfor ii=1:M
    x(:,ii) = U\(L\(P *b(:,ii)));
end
toc%Elapsed time is 0.284230 seconds.

If the problem is that the b(:,ii) are not simultaneously available, you can do the LU decomposition explicitly, as in my example above.

[Michal]

Thanks!!! Is there any scenario where use of decomposition is really benefitial?

[Bruno Luong]

When the rhs b is not available at once?

[Michal]

So, there is no reason to use parfor and decomposition at the same time? Is this the reason why decomposition does not support saving?

[Bruno Luong]

To exploit cpu of computer cluster?

[Bruno Luong]

" why decomposition does not support saving?"

I have no idea. I even don't know why there is a saving there.

[Walter Roberson]

If I recall correctly, at some point a Mathworks employee indicated the decompositions are created by calling into external libraries, and that the decompositions live in memory allocated in those libraries.

Those libraries do not provide serialization and deserialization functions, so in order to be able to save and load into a worker, Mathworks would have to write serialization and deserialization routines for the external libraries -- and would have to revise them any time the details of the external library changed.

[Michal]

This is my primary problem: I would like to use precomputed dA = decomposition(A) at parfor loop, where each b(:,ii) -> b_ii is computed independently, but this approach is impossible due to the fact, that parfor does not work in this case.

So in may case the problem looks like:

dA = decomposition(A); 
parfor ii=1:M
    b_ii = rand(N,1); % independent generation of the ii-th rhs b
    x(:,ii) = dA\b_ii;
end

[Walter Roberson]

Use parfevalOnAll to write the decomposition once per worker -- after which all iterations on that worker can access it.

Yes, some work is duplicated, being done as many times as you have workers. But it does not have to be done every iteration.

[Bruno Luong]

If you display dA you can get information about which kind of decomposition it is finally used in your matrix, so you can use just the standard linear algebra function too achieve the same thing, if you know a bit of algebra.

Personally I don't see a need to use decomposition in fact.

btw you should also consider using INV

iA = inv(A); 
parfor ii=1:M
    b_ii = rand(N,1); % independent generation of the ii-th rhs b
    x(:,ii) = iA*b_ii;
end

[Walter Roberson]

inv() is less robust than decomposition.

[Bruno Luong]

This is just a non-sense claim, blindly believed by too many people. The difference is negligible (in the order of effect of floating point truncation of the rhs or matrix).

[Michal]

Oops ... using INV is surprisingly fast in this case ... Thanks!

[Walter Roberson]

Cleve has written a book that recommends against matrix inversion. It is referenced in one of the lapack papers, http://www.netlib.org/lapack/lawnspdf/lawn27.pdf

[Bruno Luong]

I have another take, It is very simple

each column iA(:,j) where iA designate inv(A) is solution of

A*xj = ej

at the numerical precision error.

If you believe the algorithm behind inv is a "bad" one, the simply compute

iA = A \ eye(size(A)).

Personaly I use INV many many time and I nevr encount precision issue more than "\".

If you don't believe that the solution of

A*x = b

can be approximate by

x = sum(b(j)*xj)

from the most fundamental property required by the linear operator, then you don't believe the whole algebra computation with finite precision, or the problem is elsewhere (your system is ill-conditioning to start with).

[Matt J]

This is my primary problem: I would like to use precomputed dA = decomposition(A) at parfor loop, where each b(:,ii) -> b_ii is computed independently, but this approach is impossible due to the fact, that parfor does not work in this case.

If the b are computed independently, it would be better to postpone the inversion until after the parfor loop

parfor ii=1:M
    b_ii = rand(N,1); % independent generation of the ii-th rhs b
    B(:,ii) = b_ii;
end

x=A\B;