-3

I wanna solve #n linear equations AX=bi(for #n b's) in Matlab which b changes in a loop and A is constant.

One way which is fast, is to compute the inverse of A before the loop and in the loop body just get X from inv(A)*b, but because the matrix A is singular, I get an awful answer! Of course, the numerical solution A/b gives a good answer, but the point is that it takes a long time to compute #n different X's in #n loops.

What I want is a solution which can be both accurate and fast.

Chris Barlow
  • 3,274
  • 4
  • 31
  • 52
Mary
  • 19
  • 6
  • `A/b` doesn't make sense for what you've written. Do you mean `A\b`? Also, if A is actually singular, there IS no good solution, and `A\b` is better than `inv(A)`, but the problem is with your formulation, not the code. – Peter Oct 01 '13 at 21:27
  • sorry! I did mean A\b and actullay i've used that and got a good answer, but my problem is the slowness of my code. coz i wanna solve the equation for 10000 different b's which A is the same. – Mary Oct 02 '13 at 18:56
  • In that case, I think you will find my answer fits the bill. You can factorize `A` once and solve very quickly (back substitution via multiple possible methods) when you get a new `b`. – chappjc Oct 02 '13 at 20:39

1 Answers1

3

I actually think this is a good question, typos and issues of matrix singularity aside. There are a few good ways to handle this, and Tim Davis' factorize submission on MATLAB Central covers all the angles.

However, just for reference, let's do it on our own in native MATLAB, starting with the case where A is square. First, there are the two methods you suggested (inv and \,mldivide):

% inv, slow and inacurate
xinvsol = inv(A)*b;
norm(A*xinvsol - b ,'fro')

% mldivide, faster and accurate
xref = A\b;
norm(A*xref - b ,'fro')

But if like you said A does not change, just factorize A and solve for new b! Say A is symmetric positive definite:

L = chol(A,'lower'); % Cholesky factorization

% mldivide, much faster (not counting the chol factorization) and most accurate
xcholbs= L'\(L\b); %'
norm(A*xcholbs - b ,'fro')

% linsolve, fastest (omits checks for matrix configuration) and most accurate
sol1 = linsolve(L, b, struct('LT',true));
xcholsolv = linsolve(L, sol1, struct('LT',true,'TRANSA',true));
norm(A*xcholsolv - b ,'fro')

If A is not symmetric positive definite, then you'd use LU decomposition for a square matrix or QR otherwise. Again, you can do it all yourself, or you can just use Tim Davis' awesome factorize functions.

chappjc
  • 30,359
  • 6
  • 75
  • 132
  • Good to hear. What solution did you adopt? – chappjc Oct 03 '13 at 23:24
  • just I used the method factorize existing in the link u said. – Mary Nov 22 '13 at 16:50
  • @user2374418 Nice. I use it too. Could you accept the answer since that was the primary suggestion? – chappjc Nov 22 '13 at 17:50
  • what do u mean? used that coz it worked good! but still serach for a better solution to solve massive size matrices. – Mary Nov 23 '13 at 18:19
  • @user2374418 Well, I thought you were just looking for a way to solve when b_i changes but not A, and it sounds like this worked. I didn't realize there was an open issue of massive matrixes (not mentioned in the question). I'd suggest creating a new question emphasizing massive matrixes, explaining specifically how your current solution fails. But why not accept this answer since it addresses the problem as described in the question as you have said twice now? That's usually how it works - specific questions, specific answers. – chappjc Nov 23 '13 at 18:58
  • yeah, that time my question was what I said, but now have another question.I'll ask it when completely know what i want. tnx for ur attention. – Mary Nov 24 '13 at 20:14
  • 1
    @Mary - Don't be a tease. Re-accept his answer. If you want to ask another question, **make another post**. – rayryeng Jul 03 '14 at 00:05
  • I don't get your meaning? :| – Mary Jul 04 '14 at 09:28
  • @Mary - Look at the conversation that you had with chappjc. You accepted his answer, then unaccepted his answer. You unaccepted his answer because you said you have another question to ask. Focus on **one question** at a time. Accept his answer as he has adequately answered your inquiries **on this post**. If you have another question, **make another post**. – rayryeng Jul 05 '14 at 07:16
  • 1
    There is but [one appropriate response](http://i.imgur.com/uF8vKog.jpg) to that. Step 5: Quit Stack Overflow. :( – chappjc Jul 06 '14 at 17:26
  • @chappjc - I'm sorry. I tried my best to get this person to reverse their decision and re-accept your answer. This person obviously can't tell her head from her @$$ so I'll just leave it as it is. FWIW, I do thoroughly enjoy reading your answers, and I upvote whenever I can. One of the best things I learned from you was `accumarray`. Mary - do us a favour and do Step #5 of chappjc's instructions. Please and thank you. – rayryeng Jul 09 '14 at 22:57