Algorithm to extract linearly dependent columns in a matrix

Is there any general or standard approach to extract columns that are linearly dependent from the given matrix ?

Thanks and any help is apperciated !

Answer 1

[John D'Errico]

It is difficult to answer this. Consider a simple matrix:

format short g
A = rand(4,5)
A =
      0.73796      0.36725     0.042904    0.0056783      0.45642
      0.13478      0.54974      0.60439      0.39645      0.36555
       0.4903      0.25921      0.63407      0.77345      0.95769
      0.43373      0.27658       0.6462      0.25778      0.23421

The matrix (since it is random) will be of full rank, thus 4 in this case.

EVERY column is linearly dependent. That is, We can write every column as a linear combination of the other 4 columns. I can argue problems exist with other matrices too.

A = magic(6)
A =
    35     1     6    26    19    24
     3    32     7    21    23    25
    31     9     2    22    27    20
     8    28    33    17    10    15
    30     5    34    12    14    16
     4    36    29    13    18    11
     
rank(A)
ans =
     5

So A here has rank 5. But again, we can write any column of A as a linear combination of the other 5. So which column is linearly dependent? They all are!

Perhaps you might get something out of the null space vector(s).

Anull = null(A);
Anull = Anull/Anull(1)
Anull =
            1
            1
         -0.5
           -1
           -1
          0.5

This gives us the linear combination of importance as:

A(:,1) + A(:,2) - 0.5*A(:,3) - A(:,4) - A(:,5) + 0.5*A(:,6) = 0

We can now solve for ANY of those columns, in terms of the others.

How it helps you, I don't really know, because I have no idea what you really want to do.

If I had to guess, what you really need is to learn enough about linear algebra, and perhaps what a pivoted QR decomposition might provide. Because once you have that pivoted QR, you also have enough to do almost anything you want to do.

[Q,R,E] = qr(A,0)
Q =
    -0.017647      0.64381     -0.20699      0.22818      0.49018          0.5
     -0.56472     -0.11589     -0.45002      0.59421     -0.33476  -3.7638e-16
     -0.15883      0.52674       -0.446     -0.47355     -0.15542         -0.5
     -0.49413   -0.0017098      0.37629      0.14508      0.58583         -0.5
    -0.088237      0.52963      0.62872      0.19923     -0.52605  -4.2484e-16
     -0.63531     -0.11878      0.13728     -0.55665    -0.059779          0.5
R =
      -56.666      -16.377      -42.107       -33.53       -33.53      -35.224
            0       53.915       18.611       30.231       30.676       29.048
            0            0       32.491      -7.9245      -8.9182      -11.289
            0            0            0       10.101       5.6138     -0.56312
            0            0            0            0       5.1649      -5.1649
            0            0            0            0            0  -2.0136e-15
E =
     2     1     3     6     4     5

We can use the QR to tell us much about the problem, as well as efficiently and stably provide all you need.

First, notice that Q(6,6) is essentially zero, compared to the other diagonal elements.

As well, the QR tells us that it decided column 5 (that is the last element of E) is the one it wants to drop out.

[HN]

Thank you so much,

What if the matrix is a random 6 by 3 ? will still visually checking the value of diagonal elements work ? Or will only follow the rank estimator at this time ?

Thank you !

[John D'Errico]

A random 6x3 matrix will almost never be rank deficient, at least, not if it is truly random. A slightly better measure of rank deficiency is rank. See how it will work in an example. I'll create a 6x4 array to make things interesting.

A = [rand(6,1),rand(6,1)*rand(1,3)];

The matrix A is known to be a rank 2 matrix based on how I created it. As such, again, I could pick column 1 and ANY of the other three columns, thus discarding any 2 of the columns 2:4. I'll show you how this will work.

rank(A)
ans =
     2

So we learn from rank the matrix is rank 2. Therefore we know we need to find two independent columns.

[Q,R,E] = qr(A,0);

diag(R)
ans =
      -1.5743
      0.77783
  -1.3837e-16
  -3.8753e-17

E
E =
     1     3     2     4

Did it work? QR does agree with rank. The pivoted QR is pretty much as good as rank, but rank gives you a hard number, whereas here we are counting the number of diagonal elements that are close to zero. So rank is a much better tool to just predict the numerical rank of a matrix.

Again, look at the diagonal elements of R. Here we see two elements that are essentially zero compared to the others.

Therefore we discard the last TWO columns indicated by the last two elements of E, as being the "most" dependent columns. What remains are columns 1 and 3.

A(:,E(1:2))
ans =
      0.55352      0.95123
       0.9532      0.68952
      0.76997      0.17197
      0.53268     0.037538
      0.60707      0.25574
       0.1352      0.26301

So the pivoted QR gives you what you need, and depending on what you needed to do, the pivoted QR might also provide other useful information. For example, you can use it for a projection into a rank 2 subspace. I would use it (and have used it) for example if I were working with triangles in 3 or 4 dimensions, where I then need to convert the problem into a 2-dimensional one.

[HN]

Thank you so much John D'Errico . what is case if our A is 3 x 6 ? R has only three diagonal elements to check if its value is comparatively close to zero to be dicarded ? Example:

    A =
         0    1.0000         0   -0.0073         0   -0.2499
   -0.8660   -0.5000         0   -0.0320    0.0554   -0.2416
    0.8660   -0.5000         0    0.0283    0.0491   -0.2434
    
[Q, R, E] = qr(A,0)

Q =

   -0.8165    0.0000    0.5774
    0.4082   -0.7071    0.5774
    0.4082    0.7071    0.5774


R =

   -1.2247    0.0000    0.0060    0.0045    0.0427         0
         0    1.2247   -0.0013    0.0427   -0.0045         0
         0         0   -0.4243   -0.0063    0.0603         0


E =

     2     1     6     4     5     3  

The result is correct since I know the rank and independent columns from the physical meaning. However, discarding dependents based on the diagonal element of R seems to be confusing. How this can be explained ? How the procedure can be explained in general, like pseudocode or something ?

[Bruno Luong]

@HN, QR is equivalent (in practice MATLAB uses given rotation) to the Gram-Schmidt orthorgonalization process on column vectors of the matrix.

It build a chain of growing subspaces

V1 = span({c1})

V2 = span({c1 c2})

...

V6 = span({c1 c2 ... c6})

where c1, ... c6 are columns of the matrix A.

When operates with permutation (3rd output argument p requested)

[Q,R,p] = qr(A,'vector')

at every iteration of the process, QR selects the "best" columns in the sense that the new vector ck has the largest orthogonal components (per unit) to the previous space, and this orthogonal component is R(k,k), element of the diagonal of R. When the matrix rank r is reached, R(k+1,k+1) = 0, because of numerical error we never test == 0 but with some tolerance. The (r) previous vectors form an (r) independent columns vectors.

That's the meaning of diag(R).

The index of the columns selected as above avec p, meaning

c1 = A(:,p(1))

c2 = A(:,p(2))

c6 = A(:,p(6))

p(1:r) are therefore the indices of independent colums of A .

[HN]

Thank you and you are amazing! Quite generous expert I've ever seen!

Thank you again !

[John D'Errico]

Actually, I would assume QR uses Householder transformations, not that it really matters. In order to know, I would need to check LAPACK, which will be what MATLB uses for the QR.

[HN]

Is LAPACK written in mulitiple language ?

The following document tells Matlab has used Householder algorithm . http://www.seas.ucla.edu/~vandenbe/133A/lectures/qr.pdf

[HN]

Bruno Luong,

Would you give me one more help on how E (output of qr) is obtained? and its mathematical relation with independent columns ?

Thank you

[Bruno Luong]

I already explained in words in my previous comments:

"at every iteration of the process, QR selects the "best" (remaining) columns in the sense that the new vector ck has the largest orthogonal components (per unit) to the previous space"

The third output E is a permutation matrix (or vector) that reflects this selection.

When pivoting is requested, the diagonal of R is automatically sorted in descending order in absolute value due to the selection.

If you really are willing to learn about the details, check my vectorized QR implementation here and read the if block under IF PIVOTING

[HN]

Thank you for directing me to the detail. I am now studying how it is implemented.

Thank you !

[HN]

Bruno Luong, May I know the reason why R is used used both in the input and output of the function " MultipleQR()" ? Changing to different variable name causes some indexing issue. That's why I am asking . Thanks

[John D'Errico]

LAPACK was originally written in Fortran as I recall, but it was quickly put into c++.

You can find it on netlib.

http://www.netlib.org/lapack/

Personally, my favorite source for these basic algorithms was the LINPACK User's Guide. This was a well written and highly readable resource when I was just learning numerical linear algebra.

https://epubs.siam.org/doi/book/10.1137/1.9781611971811

[HN]

Thank you John ! I will read it !

[Bruno Luong]

"May I know the reason why R is used used both in the input and output of the function " MultipleQR()" ? "

This is an implement detail I chose for mainly two reasons.

• It saves memory by not duplicate the matrix
• The QR factorization can be rewritten as Q'*A*P = R

The later identity meaning I search sucessive the basis orthonormal basis Q and row-permutation A*P of A such that the projection of the second on the basis is new matrix is upper diagonal (R). The the implementation just assume A is R before iteration (thus non upper triangular) and update R while the iterations until it reaches the triangular form, using Householder transformation at each iteration to make elements bellow the diagonal become 0s.

[ahmadreza hormozi]

i cant understand

[Nguyen Le]

How fast is this method compare to build a chain of growing subspaces, check for the rank and only keep the columns that increases the rank of the subspaces

[Matt J]

@Nguyen Le I suspect that method would be a lot slower because it could require O(N) svd operations to compute the rank. Additionally, it would be difficult to decide what numerical tolerance parameter to choose for the submatrix rank computations. Consider for example,

A=[1e-20 1;
   1e-20 1];

B=[1e-20 2e-20;
   1e-20 2e-20];

Numerically, the first column of A should be considered to have rank 0, whereas the same column in matrix B should be considered to have rank 1. But there's no way to know that without somehow choosing rank()'s tolerance parameter adaptively based on the entire matrix. With the selections below, we get the opposite results from what we would want:

rank(A(:,1))

ans = 1

rank(B(:,1),1e-10)

ans = 0

Answer 2

[Matt J]

See this FEX contribution

https://www.mathworks.com/matlabcentral/fileexchange/77437-extract-linearly-independent-subset-of-matrix-columns

[HN]

Thank you Matt J

I just found it a bit earlier before you post it here and it clears everything for me !

Thank you

Answer 3

[Bruno Luong]

Test matrix (10 x 6) with rank 4

M = rand(10,4)*rand(4,6)

Automatic selection of independent columns of M

[Q,R,p] = qr(M,'vector');
dr = abs(diag(R));
if dr(1)
    tol = 1e-10;
    r = find(dr>=tol*dr(1),1,'last');
    ci = p(1:r) % here is the index of independent columns
else
    r = 0;
    ci = [];
end

% Submatrix with r columns (and full column rank).
Mind=M(:,ci)

% Those three rank estimation should be equals 
% if it's not then the cause if MATLAB selection of tolerance for rank differs with the above
% and usage of more robust SVD algorithm for rank estimation
rank(Mind)
rank(M)
r

[HN]

Thank you Bruno Luong

[HN]

Bruno Luong,

Is this problem actually the theoretically impossible problem or something that doesn't have thoretical explantion? Your program above works well for any arbitrary matrix I tested and I heard many people saying it is impssible to extract linear independent columns of matrices. How would this constradiction come to agreement? I am just wondering since I am not mathematician !

Regards !

[Bruno Luong]

I can't comment on other people saying without the context. You better returns the question to them.

I don't understand your question of "doesn't have thoretical explantion".

I, John and Matt show you HERE with more or less the same method with explanation that is doable, providing the limitation of numeric rank/vector dependency is a gray notion in finite-precision calculation.

This algorithm works with solid and rigorous algebra logic background, it's not an emperical algorithm that is impossible to explain.

[HN]

I compltely agree John's explantion tells its doable. I am actually working in the field of robotics and usually had a fat complementary subspace matrix. So, I am looking for a method to recognize dependent and independent columns from it and many experts in my field told it is impossible. Moreover, Matt in his similar submission told "but to this day, I have not been able to find a rigorous mathematical explanation for why it works". I actually left similar question to Matt but didn't get a reply yet !

That is why I asked you !

Thank you

[Bruno Luong]

OK here is MY claim: The algorithm I post above works and gives the right answer (otherwise I don't post an anwer to your original question).

If other people statements disturb you please direct your question to THEM. I have no further comment.

[Matt J]

I have used this solution for a long time and never seen it fail, but it seems to rely on the assumption that the rows of the matrix abs(R) from the QR decomposition are maximized at the diagonal. Otherwise, why is it sufficient for the code to look only at the diagonal elements dr? Again, I have never seen this fail to be true, but it would be nice to have a source from mathematical literature that confirms this property of the QR decomposition.

[HN]

Yeah, it works well but you may check Burno's implementation for detail. I learnt QR of mathworks is based on householder algorithm.

Thanks !

[Matt J]

but it seems to rely on the assumption that the rows of the matrix abs(R) from the QR decomposition are maximized at the diagonal. Otherwise, why is it sufficient for the code to look only at the diagonal elements dr?

As an example of why this might not always be true, consider the example from this wiki page

https://en.wikipedia.org/wiki/QR_decomposition#Example

When we do a QR decomposition without permutation, some of the off-diagonals dominate the diagonals,

>> M=[12 -51 4;6 167 -68; -4 24 -41];
>> [Q,R]=qr(M,0); R

R =

  -14.0000  -21.0000   14.0000
         0 -175.0000   70.0000
         0         0  -35.0000

So, the question is why doesn't this happen when we include permutation? Is it theoretically guaranteed to always be this way, and if so why?

>> [Q,R,p]=qr(M,0); R

R =

  176.2555  -71.1694    1.6680
         0   35.4389   -2.1809
         0         0  -13.7281

[Bruno Luong]

Matt: "but it seems to rely on the assumption that the rows of the matrix abs(R) from the QR decomposition are maximized at the diagonal."

Yes this statement is true

abs(R(i,i)) >= abs(R(i,j)) for all j>=i

But this is not only an assumption, it is a consequence of the pivoting selection of QR algo with permutaion/.

Proof

At iteration #i, the algorithm selects the remaining (n-i+1) columns of A that maximizes the projection on orthogonal of the subspace V_{i-1} (notation used in my explanation below John's answer), which is abs(R(i,i)).

For simplification let denote W_{k} := orthogonal of the subspace V_{k}.

So later on at the iteration step j > i, abs(R(i,j)) is the norm of the projection of A(:,p(j)) on span(Q(:,i)).

As span(Q(:,i)) is included in W_{i-1} by construction, thus

norm(projection x on span(Q(:,i))) <= norm(projection x on W_{i-1}) for any vector x in R^n.

Apply this for x = A(:,p(j)), we get

    abs(R(i,j)) <= norm(projection A(:,p(j)), on W_{i-1}) <= norm(projection A(:,p(i)), on W_{i-1}) = abs(R(i,i))

This is the proof of why diag(R) dominates off-diagonal terms row-wise, therefore you can remove the statement as "assumption".

[Bruno Luong]

By similar proof we can show that QR factorization with pivotting satisfied

abs(R(j,j)) <= abs(R(i,i)) for j>=i

In other word abs(diag(R)) is decending sorted.

It can easily see that

the absolute value of determinant of A(:,p(1)) projected on span{Q(:,1)} is abs(R(1,1))

the absolute value of determinant of A(:,p(1:2)) projected on span{Q(:,1:2)} is abs(R(1,1)*R(2,2))

...

the absolute value of determinant of A(:,p(1:n)) projected on span{Q(:,1:n)} is abs(R(1,1)*R(2,2)*...*R(n,n)

We also know when determinant ~= 0, column vectors are independant.

All these observations lead to this fact:

A(:,p(1:r)) are independent with r is the last index such that abs(R(r,r)) > 0.

r being the rank of A, since the product after n will result 0.

QED

[Matt J]

Fantastic!!!!

[Henry Wolkowicz]

Are you sure about this algorithm? Are you sure it does what you claim?

I have a 500x2000 matrix and I used both [Q,R,E]=qr(X,0); [~,R,E]=qr(X,0)

and then I plotted plot(abs(diag(R(:,1:500)))) and it definitely!!! was NOT monotonically decreasing!!!

In fact there were some elements on the diagonal that were 'tiny' and so the rank estimation was wrong as the tiny elements were in the middle and not at the bottom.

Has the QR changed at all in recent years???

[Henry Wolkowicz]

I now see the problem, i.e. qr behaves differently with sparse input matrix X. It is fine with X=full(X)

[Bruno Luong]

When you apply qr with permutation on sparse matrix S

[Q,R,p] = qr(S,'vector')

MATLAB returns the permutation to have R having "good" sparse pattern, and does not make the diagonal

abs(diag(R))

decreasing (this requirement usually makes R fully filled with non-zeros elements).