# decomposition used in parfor

27 views (last 30 days)
Michal on 18 Aug 2022
Commented: Matt J on 19 Aug 2022
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???
Edric Ellis on 19 Aug 2022
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 on 19 Aug 2022
So this work around makes the same decomposition performed by all worker instead of once, right?

Matt J on 18 Aug 2022
Edited: Matt J on 18 Aug 2022
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.
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

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!