Linear System Solve in MATLAB
显示 更早的评论
Dear all,
Thanks in advance for your help!
I am dealing with a di-banded lower triangular system where the full right hand side vector is not available to begin with. For instance, I would solve the first equation and the right hand side of the second equation will depend on the solution from the first equation.
I understand that we could use for-loop to loop over this but it seems like the computational expense would be much worse, when compares with the vectorized version...
Do you think the performance will suffer greatly if we do for-loop in this case?
Thanks for your help once again!
Sincerely,
Tae
4 个评论
John D'Errico
2022-3-21
Not totally clear.
You have a banded matrix A. Is A VERY large? Is it stored in sparse form? How large is the matrix A? How much time was required to perform the initial solve?
You solve the problem A*x = b1, where b1 is known. Given the solutino to this problem, you create a new right hand side vector, call it b2, and then solve the problem, with the same matrix A, where A*y = b2.
Is that correct?
Taehun Kim
2022-3-21
编辑：Taehun Kim
2022-3-21
Thanks John for your reply! I am sorry that my explanayion was not clear. So, I have n (let's say around 200) equations where for a given row n, we have an equation where 
. The claim is that at 
, we have the equation 
and we know 
(given as a boundary condition), and 
. Then, the idea is to simply solve for 
. Again, 
's and 's are parameters and 
's are the unknowns in this case. For 
now, we have 
. The problem is that 
is a nonlinear function of and can only be obtained after solving the equation for case. So, I am forced to solve each equation sequentially, rather than being able to vectorize and solve all at once, if the vector 
is available to begin with... I am worried if the sequential solve would be very slow when comapred against the backslash?
. The claim is that at , we have the equation and we know (given as a boundary condition), and . Then, the idea is to simply solve for . Again, 's and 's are parameters and 's are the unknowns in this case. For now, we have . The problem is that is a nonlinear function of and can only be obtained after solving the equation for case. So, I am forced to solve each equation sequentially, rather than being able to vectorize and solve all at once, if the vector is available to begin with... I am worried if the sequential solve would be very slow when comapred against the backslash?
John D'Errico
2022-3-21
But b_n is a nonlinear function of x_(n-1). So regardless of whether it might be slow or fast, I don't see you having any choice. The looping time is almost irrelevant. If you are worried about that, use the profiling tool to look at where the time is spent. I would conjecture it will be mostly in the nonlinear solve when you need to compute bn at each iteration.
Unless you are doing this entire operation millions of times, I would just say to get a cup of coffee and enjoy the few seconds of relaxation.
Taehun Kim
2022-3-21
采纳的回答
So, I have n (let's say around 200) equations...I am worried if the sequential solve would be very slow when comapred against the backslash?
The two are not comparable because backslash cannot handle non-linear equations, which is what you have. In any case, the performance impact due to a 200-iteration loop seems insiginificant. If the b_i() are very simple the loop will run very fast and 200 steps is nothing. If the b_i() are complex, then the overhead of the loop will be insignificant next to the nonlinear operations anyway.
3 个评论
Taehun Kim
2022-3-21
Does this sound like a "light computation"?
Hard to say, because your example doesn't illustrate the nonlinear part of your computation. But as I said also, it doesn't matter whether the computation is light or heavy. In either case, it won't be the for-loop that is slowing you down.
Taehun Kim
2022-3-21
更多回答（1 个）
If you're trying to some the same linear system repeatedly for different right hand side vectors consider using decomposition on the coefficient matrix (so MATLAB performs the analysis and preparation work for your system once) and solve it repeatedly.
b = rand(200, 1e4);
A = randi([-10 10], 200);
x1 = zeros(size(b));
x2 = zeros(size(b));
tic
for k = 1:width(b)
x1(:, k) = A\b(:, k);
end
toc
tic
D = decomposition(A);
for k = 1:width(b)
x2(:, k) = D\b(:, k);
end
toc
11 个评论
Taehun Kim
2022-3-21
Humm I must miss the context but reading Steve's code I though "why not just invert directly?"
b = rand(200, 1e4);
A = randi([-10 10], 200);
x1 = zeros(size(b));
x2 = zeros(size(b));
tic
for k = 1:width(b)
x1(:, k) = A\b(:, k);
end
toc
tic
x3 = A \ b;
toc
norm(x3-x1,'fro')/norm(x1,'fro')
Taehun Kim
2022-3-21
Thanks Bruno! Yes. I think the problem that I motivated in this thread was as the following: We have a set of n equations where in the matrix form, we can write the coefficient matrix as a lower triangular matrix with two diagonal entries (the main and the one below). For the system of 
, we have a situation where the right hand side vector 
actually depends on the solution from the previous equation. For instance, from the first equation, we get 
but for the second equation, 
depends on 
and so forth! So, I had to go through the for-loop and I was worried that I was going to get a hit in the performance. From the answers and responses from many above, the performance hit from the for loop is minimal when compared to the performance hit from evaluating the non-linear right hand side! I hope this was a good recapture of the information discussed above! :)
, we have a situation where the right hand side vector actually depends on the solution from the previous equation. For instance, from the first equation, we get but for the second equation, depends on and so forth! So, I had to go through the for-loop and I was worried that I was going to get a hit in the performance. From the answers and responses from many above, the performance hit from the for loop is minimal when compared to the performance hit from evaluating the non-linear right hand side! I hope this was a good recapture of the information discussed above! :)But then why not just use inv() v(despite all the warning I think it's just a bullshit, sorry TMW)
b = rand(200, 1e4);
A = randi([-10 10], 200);
x1 = zeros(size(b));
x2 = zeros(size(b));
tic
for k = 1:width(b)
x1(:, k) = A\b(:, k);
end
toc
tic
D = decomposition(A);
for k = 1:width(b)
x2(:, k) = D\b(:, k);
end
toc
x4 = zeros(size(b));
tic
E = inv(A);
for k = 1:width(b)
x4(:, k) = E*b(:, k);
end
toc
Taehun Kim
2022-3-21
Matt J
2022-3-21
I do not see how inv() helps, since the equation supposedly has the nonlinear form A*x=b(x).
Taehun Kim
2022-3-21
Taehun Kim
2022-3-21
Bruno Luong
2022-3-21
I simply comment on Steve for-loop answer with repeat A \ b for the same A, so I replace that statement by inv(A)*b which is strictly equivalent to his code and fatest so far. I don't know the context of non-linear and how the equation is discretized.
Steven Lord
2022-3-22
In my example I knew all the b vectors from the start, but if you were solving for one vector and using that solution to create the next right hand side vector you could use decomposition.
Bruno Luong
2022-3-22
编辑：Bruno Luong
2022-3-22
But inv is twice faster than decomposition. At the end when using "\" on decomposition it numerically close to applied multipy to inv. If someone disagree they ough to explain the math to me.
If for some reason inv using worse algorithm then just call
E = A \ eye(size(A); % = inv(A)
The warning of MATLAB on using inv is jist non-sense to me.
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
