solve A * A^T *x = b by QR factorization

Hi guys,

I would like to solve the following system in Matlab with QR factorization of A^T

A * A^T *x = b

with A \in R^(m x n) , x \in R^m and b \in R^m. In my case is n>m.

I know how to do a QR factorization of A and solve Ax=b. I hope you can help me solve my problem.

Answer 1

[John D'Errico]

Can you solve this in two steps? Thus, if you know

A*A' * x = b

Does the associative law apply to matrix multiplication? (Yes.)

Then is it not true that this is just:

A* (A' * x) = b

Now, I could swear I recall you telling us that you already knew how to solve PART of that problem. That is, can you solve the problem

A*y = b

So IF you did that, then what could you say about y? How is y related to x?

Now solve the resulting problem. Even better, you already have a factorization for A', since you have a factorization of A. That is, if

A = Q*R

then just take the transpose. What is the transpose of the product of two matrices? Use a search engine if you do not know that already.

[Lukas Müller]

Since I need a QR factorization of A' I would say

A'=Q*R 

and I get

A=R'*Q'

In total:

R'*Q'*Q*R*x=b

which is equivalent to

R'*R*x=b since Q'*Q=I because Q is orthogonal.

Now I could calculate

pinv(R) and pinv(R') in matlab

and solve my equation by

x=pinv(R)*pinv(R')*b;

But is that efficient? Or is is better to solve the problem in two parts like you suggested (with A*y = b , y=A'*x)?

I'm dealing with large matrices.

[John D'Errico]

NO! You do not need to use pinv. You could do so, but if you were going to do that, then just compute pinv(A), and be done with it.

PA = pinv(A);

x = PA' * (PA * b);

However, since A is apparently large, computing the pseudo-inverse of A may be time consuming, more so than the QR. We are told that A has more columns than rows. So, computing the QR as:

[~,R] = qr(A');

will produce a matrix R that is upper triangular. Use the ~ form there to not return Q, since you don't need Q. It will look like this:

A = rand(2,5);
[~,R] = qr(A');
R
R =
      -1.5455     -0.91519
            0      0.41617
            0            0
            0            0
            0            0

Because Q is orthognal, you can reduce it to something of the form:

R'*R*x = b

Can you do this in one line of code, without an intermediate vector y? Of course. Just use backslash twice.

x = R \ (R' \ b);

In each case, the solve is efficient, since MATLAB can figure out that the matrix R (or R') does not need a factorization, so each of those calls is just an O(n^2) operation.

Note the careful use of parens there, as the alternative way of writing it without parens is not as efficient.

x = R \ R' \ b;

While the two would be in theory equivalent, the second form solves R\R' FIRST, then applies that to the vector b.

[Lukas Müller]

Thanks a lot :)

[John D'Errico]

To show that MATLAB does know how to solve triangular systems efficiently, here is an example:

A = rand(5000);
[L,U] = lu(A);
b = rand(5000,1);
timeit(@() U\(L\b))
ans =
      0.26981
timeit(@() A\b)
ans =
       2.6015
       

Once you have R from the QR, the solve itself is quite efficient.