MATLAB Answers

0

Solving linear equations with large matrices

Asked by Srijeet Tripathy on 12 Apr 2019
Latest activity Answered by John D'Errico
on 12 Apr 2019
I have a simple equation of the form Ax = B; where A is a large dense matrix(nXn) and B is a column matrix (nX1).
A is hermitian but not symmetric positive definite.
I am proceeding through simple \ + 'qr' decomposition. However this method is quite inefficient as it is quite time consuming.
Moreover, the command is nested within two for loops each running ~300 times. This I effectively have to solve the linear equation 300x300 times.
Can anyone suggest a faster and better method?

  6 Comments

John asked several good questions; I want to focus on one. Is the matrix constant, with only the right hand side changing in the loops? If so, consider creating a decomposition of the matrix before the loop and solving using that decomposition inside.
Alternately if the matrix is changing inside the loops but only slightly and you expect the solution at one iteration is "close to" the solution at the next, maybe an iterative solver like gmres would help, using the solution from one iteration as the initial guess for the next. You may be able to take adantage of the structure of the coefficient matrix as well. If you can compute the product of the matrix with a vector without explicitly creating the matrix, you may be able to save some memory by passing a function handle into gmres or the other iterative solvers.
MATLAB already does the linear algebra on big problems in parallel, so it will automatically farm out the computation to as many cores as you have, or at least as many as you allow it to do so.
So check this at the command line:
maxNumCompThreads
It will tell you how many cores your computer has available.
And do you have a good CPU fan? This because running all of your CPU cores flat out will cause it to kick on full, and excessive heat will throttle down your CPU.
Do you have access to a more powerful computer? Are you willing to spend money to solve the problem, perhaps renting time on somethign seriously bigger? You need to recognize that sometimes a necessary tradeoff exists between $ and your own time.
@ John •The order of the matrix varies from 800x800 up-to 4000x4000 •for a 800x800 matrix each solution takes up-to 20 to 30 seconds in comparison a cholesky factorization for a matrix of same order takes about 5 to 10 seconds. •Ram is on the plus side, I have around 40gb, no issues with memory either. •I am initially decomposing the matrix through least squares then solving it. •Matrix is not constant but value of every element in the matrix changes slightly after every iteration. @Steven •Thank you for pointing out about gmres. I will check the iterative solvers! But my matrix is dense so I hope there will be no issues there.

Sign in to comment.

1 Answer

Answer by John D'Errico
on 12 Apr 2019
 Accepted Answer

Enough information now that I'll hazard an answer.
By the way, the logic behind suggesting GMRES is NOT that your matrix is sparse, but that if it changes only by a tiny amount each time, then you might be able to gain using GMRES, because you supply a starting value. It might be worth a shot, but I'm not sure you would gain. Depending on the condition number of your matrices, the solution could change a fair amount based on only a small perturbation.
But something in your response makes no sense. A single solve for an 800x800 matrix takes 20-30 seconds? Or 90000 of them? I assume you mean the latter. A quick test on my computer shows:
A = rand(800,800);
tic,[L,U,P] = lu(A);toc
Elapsed time is 0.024238 seconds.
And a factorization will be the lion's share of the time. But that also makes no sense.
90000*.024
ans =
2160
So not 20-30 seconds, but 2000 seconds. Your computer must be either REALLY old and slow, or REALLY fast. A really slow computer makes no sense, as you say you have 40GB of RAM, and a computer old enough to be that slow would be unlikely have that much RAM.
Anyway, one thing I would point out is that QR is not necesssaily a good choice to do a square matrix solve. If these matrices are square, then why are you using QR? A pivoted LU will be nicely stable, and will be considerably faster than QR.
A = rand(4000);
>> tic,[L,U,P] = lu(A);toc
Elapsed time is 0.300080 seconds.
>> tic,[Q,R,E] = qr(A);toc
Elapsed time is 8.192923 seconds.
So perhaps when you said 20-30 seconds, you meant for a 4kx4k matrix took that long for a solve, if your computer is 3x as slow as mine. That is possible, since mine should be putting up some decent bench scores. (All 8 cores on my CPU were running on that QR factorization.)
But the point is, why are you using QR here at all? The matrix is square! Just use a column pivoted LU. The difference is huge, certainly on a 4kx4k matrix.
In fact, since your matrices are changing on every iteration anyway, then why are you using QR or even LU anyway? Just use backslash, and let it worry about then doing the final solve too. Its probably faster than doing an LU, then applying those factors to a column vector for the solve.
tic,x = A\b;toc
Elapsed time is 0.349814 seconds.
That was pretty fast, in comparison to the time the LU itself took.
My point in this is that using QR to solve the problem may well be a major factor in why your code is so slow. Just use backslash, unless there is some reason why you REALLY, REALLY need a QR. But if you would be willing to consider an iterative solver, that also makes no sense. So use backslash.

  0 Comments

Sign in to comment.