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

[John D'Errico]

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]

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?

[John D'Errico]

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]

Thanks John for your insightful comment! I was worried about this but your comment really gave me some sense of relief and assurance! I did run some small test as below where I applied to one of the answers! It seems like the performance hit is there but not too much significant as I was worrying over about it! Many thanks for your help! :)

Answer 1

[Matt J]

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.

[Taehun Kim]

Thanks Matt for your answer! I am greatly appreciative to hear from all of your perspective and get some sense or ideas on the problem! In this case, I do have some components of the right hand side pre-computed before the for-loop and I will just have to add some number to the righ hand side to update it. Does this sound like a "light computation"? Thanks alot! I also did some simple randomized case study as below... I hope when I solve the actual problem, it behaves nice too... Thanks for your help in getting me to a working direction!

%Initialize the random matrix A

a1 = rand(200,1);

a2 = rand(199,1);

A = zeros(200,200);

A= A + diag(a1,0);

A= A + diag(a2,-1);

%Initialize the random right hand side vector b

b = rand(200,1);

%Solve using backslash

x1 = zeros(200,1);

tic;

x1 = A\b;

toc;

%Solve using a loop

x2 = zeros(200,1);

tic;

x(1) = b(1)/A(1,1);

for i = 2 : 200

x(i) = (b(i) - A(i,i-1) * x(i-1)) / A(i,i);

end

toc;

Elapsed time is 0.001877 seconds.

Elapsed time is 0.006078 seconds.

[Matt J]

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]

Ah, I see... Thanks Matt for the insightful answer and sharing your experience! I really appreciate it!

Answer 2

[Steven Lord]

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

Elapsed time is 3.731868 seconds.

tic
D = decomposition(A);
for k = 1:width(b)
    x2(:, k) = D\b(:, k);
end
toc

Elapsed time is 0.216534 seconds.

[Taehun Kim]

Thanks Steven! The decomposition.m seems to be a beautifully written function!

[Bruno Luong]

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

Elapsed time is 5.297187 seconds.

tic
x3 = A \ b;
toc

Elapsed time is 0.025828 seconds.

norm(x3-x1,'fro')/norm(x1,'fro')

ans = 6.2704e-15

[Taehun Kim]

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! :)

[Bruno Luong]

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

Elapsed time is 4.255249 seconds.

tic
D = decomposition(A);
for k = 1:width(b)
    x2(:, k) = D\b(:, k);
end
toc

Elapsed time is 0.233741 seconds.

x4 = zeros(size(b));
tic
E = inv(A);
for k = 1:width(b)
    x4(:, k) = E*b(:, k);
end
toc

Elapsed time is 0.091995 seconds.

[Taehun Kim]

Inv.m clearly is the fastest in the example that you generated above! Thanks and I will take a note on this!

[Matt J]

I do not see how inv() helps, since the equation supposedly has the nonlinear form A*x=b(x).

[Taehun Kim]

Yes, Matt you are right! For this particular problem, it won't help. However, I will definietely take note of this and may apply it in any other types or problem!

[Taehun Kim]

Thanks everyone for your help!!!

[Bruno Luong]

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]

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]

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.