Some question about QR Factorization, thanks

[Hongliu]

Dear Sir:

I need solve a large sparse unsymmetrical overdetermined systems, Ax=B, where A is m by n, and m>n. A is a full rank matrix. But it is possible that some equations are inconsistent. For example, x1+x2=1, and x1+x2=2 are inconsistent. I need to detect the inconsistency in Ax=B first, if it exists, then needs not solve the equation; If it does not exist, solve the system.

I will detect the inconsistent in this way: equation Ax=B can be wrote as Rx= QTB, where rank of R is k. If the last (k+1), (k+2)… (m-k) of QTB are not equal to 0, then inconsistency exists. Otherwise, inconsistency does not exist.

I have two questions

1. Is my method of detect the insistences is correct? (obtain R by calling DGEQRF and obtain QtB by calling DORMQR ). Or other more effective routine can be used?

2. If QR factorization is sparse QR factorization in DGEQRF? If yes, any reference about it?
Thanks

[sven]

Dear Hongliu,

What is the application that leads to your equations? Since you are expecting your equations to be consistent, is it possible to model your problem so that you do not obtain duplicate equations in the first place?

But, yes, if the residual, r, of the least squares problem are zero then your equations are consistent. Numerically you should not expect to obtain exactly zero, instead norm(r)/(norm(A).norm(B)) would be of order relative machine precision.

LAPACK does not contain routines for sparse matrices - DGEQRF is intended for dense matrices. You might find a sparse QR routine on netlib at: http://www.netlib.org/

Best wishes,

Sven Hammarling.

[Hongliu]

Thanks Sven for reading and kingly reply.

Bellow is my application:
A is a matrix of m by n, and m<=n. The row vectors of matrix A are expressed by Ai=[a11 a12 … a1n].
The first (m-1) row vectors {A1, A2 … Am-1} are independent. I need to check if the Am is linearly dependent with some of other row vectors of matrix A, i.e. the mth row vector is linear combinations of the other row vectors of matrix A. the row vectors that are linearly dependent with mth row vector need to be indicated also.

I plan to use Equation (1) to check the dependence between {A1, A2 … Am-1} and Am.
[A1’, A2’ …Am -1 ‘] [ x1 x2 …xm-1]’ = Am ‘ (1)
where [A1’, A2’ …Am -1 ‘] is full rank matrix with rank m-1. Let B= Am ‘, X=[x1 x2 …xm-1]’. Equation (1) can be written as A’X=B. If Am is linear independent with other vectors, then should X has no solution.
After sparse QR factorization, RX= Q’B, where rank of R is m-1. If the last (mth, (m+1)th… nth) row vectors of Q’B are not nearly equal to 0, then inconsistency exists. Otherwise, inconsistency does not exist, solve X by RX= Q’B.

Equation (1) can have a least-squares solution also, but I hope the system can detect the inconsistence at earlier stage of solving progress, because if inconsistence exit then need not get the value of X. Hopefully it can save the running time. The size of matrix A’ is about 500 by 400.

I wrote too much this time, thank for your reading again.

[sven]

Dear Hongliu,

For problems of the size you mention, treating them as dense and using an LAPACK routine should be fine. Nowadays, 400 by 500 is not really very large.

Best wishes,

Sven.

[Hongliu]

Dear Prof. Hammarling

Thank you for your reply. I can call the routine DGEQRF of from C now.

But I have a question about what you said. You mentioned:

“But, yes, if the residual, r, of the least squares problem are zero then your equations are consistent. Numerically you should not expect to obtain exactly zero, instead norm(r)/(norm(A).norm(B)) would be of order relative machine precision.”

Why not use relative inaccuracy of r: norm(r)/norm(B) or relative inaccuracy of x : k(A).norm(r)/norm(B) to estimate the consistent ? Where k(A) is the condition number of matrix A. I can not understand why do you use “norm(r)/(norm(A).norm(B))”. Any references about it? Hope it will not take too much of your time.

Thanks and Regrads,
Hongliu

[sven]

Dear Hongliu,

I had a typo in my previous posting, norm(r)/(norm(A).norm(B))
should be norm(r)/(norm(A).norm(x)), which I hope makes more
sense!

This is a measure of the backward error in the solution and shoud
be small when you use a stable method to solve the problem,
independent of the condition of the problem. See Chapter 4 of the
LAPACK Users' Guide for background, as well as more detailed,
information and further references.

Best wishes,

Sven.

[Hongliu]

Dear Prof. Hammarling

It really makes more sense for me now!
Thank you very much, I appreciate.

Best wishes,
Hongliu