Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

Hi I need to get to the solution of Ax = b ;

Since my A is non-square, I need to use QR which is built-in in both \ and linsolve commands. I need to specify the accuracy (tolerance for error) could someone help plz. Thanks!

Answer 1

[John D'Errico]

What do you mean by "tolerance"?

A QR solution has no need for a tolerance, nor does linsolve, as they are not iterative methods.

Anyway, you cannot specify the desired accuracy of a solution in a linear system like this, solved using QR. Just wanting a more accurate solution is not relevant. The solution exists, and is unique, or there are infinitely many solutions, all of which are equally good (or poor), or there is no exact solution, and that which qr/linsolve returns is entirely valid.

Perhaps you have a singular system, or one that is nearly singular. In that case, a tolerance is still not really relevant. Yes, you could decide which singular values to discard, in case of a pinv style solution. But that is not a "tolerance" in the way tolerances are traditionally used. And if your system is singular, then again, tolerances have no real meaning.

As far as your matrix being non-square, that just means (assuming you have an over-determined system) that you will never be able to achieve an exact solution. (Well, not entirely true. It IS possible that an over-determined system has an exact solution. But that is terribly rare.) So if no exact solution exists, then specifying a tolerance is again a meaningless goal.

If your problem is non-square in the sense that it is under-determined, then again, a tolerance is irrelevant. A solution may exist, or not. If one solution does exist, then there will be infinitely many solutions, all equally good. A tolerance is meaningless here.

I think that perhaps you are being confused, thinking that the solution to a linear system is like that of a nonlinear system, where you will employ convergence tolerances, etc. Or, perhaps you think that simply by cranking down on a tolerance, you can always gain an arbitrarily good solution. (Many people seem to think that. They are flat out wrong.) Or perhaps you are thinking that a solution like that provided by linsolve/qr is like an iterative solver, like that from lsqr.

If this has not cleared things up, then I will suggest you need to explain CLEARLY what you want to do, and why you think that a tolerance is a meaningful thing to ask for here. As well, you need to explain why it is that you think what I have said above does not apply.

[Sam Sun]

Thank you very much for your time! The moment I saw QR I assumed that it was iterative so your first sentence cleared it up! I asked this question since I need to get a speed up somehow and I assumed that turning down the precision would do the trick -- now just in case you are interested-- I need to speed up this 'inversion' somehow since it gets called iteratively in my code.

All in all I really need to speed up the '\' operation here. Either I need to use 'A\B' where they are both large, sparse matrices or I need to do $A\x$ iteratively to get to my goal. Any ideas?

Once again, thanks for your time!

[John D'Errico]

If A is sparse, then is it stored in sparse form? Far too often, I see someone claim their matrix is sparse, yet they have not chosen to use sparse storage for the matrix. MATLAB cannot know that you want your matrix to be considered as sparse, unless you tell it that fact.

Next, many people think that a matrix that has perhaps 50% zeros, or even 80 or 90% zeros, is sparse, and that sparse linear algebra will then gain them speed. NOPE. Not gonna happen.

Regular (non-sparse) linear algebraic solutions are blindingly fast compared to a sparse solve. This is because a sparse algorithm needs to worry about where the non-zeros fall. A full code just pipelines everything, so that it flies like a bat out of Texas. So unless your system is truly sparse, you don't gain much. Why not? Fill-in costs you. The simple need to know where the non-zeros are, and when you can skip a multiply means you cannot easily pipeline these sparse calls. A sparse solve factorizes your problem, but in the end, you will create a much less sparse matrix. For example...

A = sprand(1000,1000,.1);
[L,U] = lu(A);
nnz(A)
ans =
       95160

nnz(L)
ans =
      484900

nnz(U)
ans =
      487464

So A is roughly 90% zero. But the factored form of A is essentially full. That means any solve using A as a sparse matrix would almost certainly be much slower than a full matrix solve.

b = rand(1000,1);

timeit(@() A\b)
ans =
          0.100780254952

Af = full(A);
timeit(@() Af\b)
ans =
          0.033592990952

So even though A is 90% zero, a sparse solve is considerably slower than the same thing using the full form of A. This is because A was not really very sparse.

A typical sparse matrix, where sparsity can benefit you, may be on the order of 99% zeros and very often even that is not very sparse. So lets try it again, for a truly sparse matrix.

A = sprand(1000,1000,.002) + speye(1000);
nnz(A)
ans =
        2997

[L,U] = lu(A);
nnz(L)
ans =
       27804

So every non-zero in A now turns into 10 non-zeros in L. Fill-in!

timeit(@() A\b)
ans =
            0.002649831952

Af = full(A);
timeit(@() Af\b)
ans =
            0.034667779952

Here the sparse solution is now better than 10x faster. Too much depends on the matrices involved to really know how to help you more here.

So, IF your matrix truly is sparse, then make sure it is stored in sparse form. (Sorry if that is already the case.) Then you need to decide how to solve it. Is A fixed? If so, then factorize it once, in advance. Is A not fixed, but not varying a lot? Now you might consider a carefully chosen pre-conditioner matrix, and an iterative method. (Not sure how much this would gain you. It has been many years since I used those tools, and I am sure they have changed.)

Next, consider tools like symamd, amd, symrcm, etc. These can decrease the amount of fill-in during a solve, and that will make a sparse solve go faster. Again, it has been many years since I was regularly working with huge sparse problems, and those tools have probably seen theoretical improvement since.

[Sam Sun]

Wow thanks again for the elaborated answer!

Indeed along with the size increase my matrix is more than 99% sparse (I did not pay attention to how sparse it was until you pointed out but it is indeed very sparse). And yes I told Matlab to do it-- otherwise I run out of RAM even before I start the 'inversion' step.

It is nice to hear advice from someone who has far more years of experience than me. A is indeed fixed-- might I ask how pre-factor is excuted and dare I ask why it is faster?

Thank you so much ! you should probably publish a book on these things (if you already did refer me there lol)

Answer 2

[Christine Tobler]

Hi Sam,

If you compute A\B for two large, sparse matrices (as you mention in your comment to John), the result will typically be dense. In that case, the memory requirements for A\B will probably be quite high, which might be even more of a problem than solving the linear system.

What are you doing with the matrix A\B once it's computed? Is there any way of doing this without explicitly computing it?

[Sam Sun]

Hi Christine,

You are absolutely correct. The RAM usage became so high as I needed to run it on the cluster. But I still've incentive to optimize it so to enable greater precision;

Here's what I need: the largest eigen value in real parts ('largestreal' or 'LR') for the following matrix:

L = C*A\B; where A, B, C are all large and sparse.

The backslash step takes really long-- so instead I tried to make a function that is a linear operator 'T' on x, a vector: i.e. Lx = Tx, so to avoid to compute L explicitly and use eig to find the function T's eigen pairs. It's again taking a great deal of time due to the backslash of A\(Bx) being called a great number of times.

Any other ideas of how to do this? Thanks!

[Christine Tobler]

You can use the decomposition object to store the QR decomposition of A, and then solve with a right-hand side inside of your function handle.

 dA = decomposition(A);
 T = @(x) C*(dA\(B*x));
 eigs(T, size(C, 1), 1, 'largestreal')

Of course, eigs expects a square matrix, so I'm assuming that size(C, 1) is equal to size(B, 2).

The decomposition object dA will need a bit more memory than just computing A\b once, but it should be much less than that used for B\A.

[Sam Sun]

This is a brilliant idea! Just amazing! Thanks!