Solving linear equations with large matrices

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?

[Alex Mcaulley]

Have you tried this?

https://es.mathworks.com/help/matlab/ref/mldivide.html

[Srijeet Tripathy]

Just as I mentioned in the question, I am using mldivide or \.

I aqm actually looking for a factorization as a preconditioner to the mldivide as th qr factorization is quite slow.

[John D'Errico]

How large of a matrix are you talking about? Note that a significant number of people think that a 100x100 matrix is large. So just calling it large tells us nothing.

What does time consuming mean to you? I.e., how much time?

Do you have sufficient memory? Is your problem that you have insuifficient memory, and are thus incurring a problem because you just need more RAM? RAM is CHEAP these days, and getting cheaper every minute.

Is your matrix fairly well-conditioned? If it is, AND you are running low on memory, then use of single precision MIGHT offer some relief.

You could use a pivoted LU as an alternative to QR.

If your matrix solve is happening inside a loop, is the matrix constant? If so, then you really don't want to be computing a factorization multiple times. Instead, you just have one problem with 90000 right hand sides. Or, if you compute the factorization in advance, then appy it to each right hand side.

Sorry if these questions seem silly, or perhaps don't apply to you, but how can we know what the real issue is unless you provide sufficient information? And there are a lot of people who might ask exactly the questions you have asked, yet have no clue what they are doing.

Remember that an nxn matrix solve is an O(n^3) operation on a dense matrix, that will require a lot of memory. If you really do have a big problem with a matrix that changes at every point in those loops, you will just need to accept the big time it will require.

[Steven Lord]

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.

[John D'Errico]

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.

[Srijeet Tripathy]

@ 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.

Answer 1

[John D'Errico]

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.