LAPACK/ScaLAPACK Development

Hello, everyone! I apologize if this is a basic question. I am looking for a way to find a maximal set of linearly independent columns. I know that dgesvd can be used to find the rank of a matrix. Can it also be used to find the actual indices of the linearly independent columns?

Thank you very much for your time!

Best regards,
Susan

You might want to use QR with pivoting for that. (DGEQP3.)

Thank you! Okay, I just want to be absolutely sure that I have this right. The Lapack call dgeqp3 takes the following:

* JPVT (input/output) INTEGER array, dimension (N)
* On entry, if JPVT(J).ne.0, the J-th column of A is permuted
* to the front of A*P (a leading column); if JPVT(J)=0,
* the J-th column of A is a free column.
* On exit, if JPVT(J)=K, then the J-th column of A*P was the
* the K-th column of A.

Since I don't know anything about the rank and columns of my matrix beforehand, I initialize the array JPVT to zero, and then pass it to db3qp3. On exit, I walk the array JPVT, and find *any* column with a non-zero entry. These columns form a maximal linearly independent set.

Or, do I need to know the rank of A, and find only the columns that have been permuted to the 1 -> rank(A) position?

Thank you very much! I am sure this is a very basic question, so I apologize for taking time!

Very best regards,
Susan

Well I indeed forgot myself how this routine worked ...
OK so I re-read the routine, I think I am up to speed now.

JPVT is an array of INTEGER. It is input/ouput.
In input, what matters is whether it is different or zero or not.
If it is equal to 0 then the column j is treated in the general case.
The general case is the following:
(1) you find the column with greatest norm, say column k
(2) you permute the k-th column with the first
(3) you factorize with the first column that is: you orthogonalize all columns from 2 to n with column 1
And you loop on these steps (1)-(2)-(3).
So at each step you extract the column with largest norm from the column space.

Now you might want to specify a few columns that no matter how large their norms are, you want to orthonormalize against them first. Then you set their JPVT in input to a value different from 0 (say 1). The first thing the routine DGEQP3 will do is to orthonormalize against these columns using standar QR factorization (without pivoting). This is the special case.

OK so that said, yes, just set up JPVT to 0 for JPVT(1) to JPVT(N).

To get the rank, you have several ways to go. You however need to realize that your question is really tricky in finite precision arithmetic (see remark after).

One way to go is the following: the columns of the R factor are ordered from the largest in norm to the smallest. So you can scan the columns of R, actually you simply need to scan the diagonal of R, and when you see a big drop in the diagonal. There is your rank, and you get the columns in JPVT(1) to JPVT(RANK).

Well, so of course all this is very tricky with finite precision arithmetic. The concept of rank is tricky in itself. To define it properly (to define the numerical rank) you need the SVD. But then you can not relate the rank to a given set of column of the matrix A has you have properly in your first post. So please try the "one way to go" and if it works for you, that's great. Otherwise ....

--julien

If possible, can you please quantify the big drop in the diagonal of R.

Not really ... this really depends on your application.
If the singular values of the matrix are 1, 1, 1, 1e-3, 1e-3, 1e-10, 1e-10, what is the rank of the matrix?
Some will tell you that the numerical rank is 7, others that the numerical rank is 5, and others that the numerical rank is 3.
We cannot answer this question for the users.
Julien.