decomposition used in parfor
Show older comments
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???
7 Comments
Bruno Luong
on 18 Aug 2022
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
on 18 Aug 2022
with R2022a Update4
Not on mine.
Bruno Luong
on 18 Aug 2022
My bad, I have the Automatic creation parallel pool option turned off, so the parfor runs like like a for.
Edric Ellis
on 19 Aug 2022
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?
@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
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.')
This causes each worker to build the decomposition for itself from the value of A.
Bruno Luong
on 19 Aug 2022
So this work around makes the same decomposition performed by all worker instead of once, right?
Answers (1)
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.
15 Comments
Michal
on 18 Aug 2022
Bruno Luong
on 18 Aug 2022
When the rhs b is not available at once?
Bruno Luong
on 18 Aug 2022
To exploit cpu of computer cluster?
Bruno Luong
on 18 Aug 2022
Edited: Bruno Luong
on 18 Aug 2022
" why decomposition does not support saving?"
I have no idea. I even don't know why there is a saving there.
Walter Roberson
on 18 Aug 2022
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.
Walter Roberson
on 18 Aug 2022
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
on 18 Aug 2022
Edited: Bruno Luong
on 18 Aug 2022
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
on 18 Aug 2022
inv() is less robust than decomposition.
Bruno Luong
on 18 Aug 2022
Edited: Bruno Luong
on 18 Aug 2022
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
on 18 Aug 2022
Walter Roberson
on 18 Aug 2022
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
on 18 Aug 2022
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
on 19 Aug 2022
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;
Categories
Find more on Linear Algebra in Help Center and File Exchange
on 19 Aug 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
