Extracting the largest set of linearly independent vectors from a set of vectors in matlab

Question

I am new to working with matlab. I couldn't find an obvious way to to extract the largest subset of linearly independent vectors from a given set of vectors.

So given a set V = [v1 v2 -- vn] where dim(vi) >> n (i=1,2,3,....) I need to find a random set of r linearly independent vectors "vi" (where r is maximal), i.e removing ALL (n-r) linearly dependent "vi" vectors from V. Sorry if this is ambiguous but can't find a better way to state it.

Answer 1

I'm going to assume that your vectors are all n-dimensional, and that we can concatenate them all into a single matrix. If we go into matrix and linear algebra, what you are looking for is the Column Space of a matrix. Simply put, the column space is defined as the set of columns in your matrix that can uniquely produce another vector in n-dimensional space. Or, it is the set of all possible linear combinations of the column vectors. As such, if you want to find the largest set of linearly independent vectors, all you have to do is determine what the column space of your matrix is.

Therefore, given your matrix V which is of size n x m, where we have m columns / vectors, with each column being of size n x 1 (or n rows), you would call the rref or the Row-Reduced Echelon Form (RREF) command. This reduces your matrix down to its row-reduced echelon form. This is the beginning of finding the column space for a matrix. You would call it like so:

[R, RB] = rref(V);

R would contain the RREF form of V and RB would contain the indices or column numbers of R that form the column space. Therefore, if you want to produce your linearly independent vectors, you simply have to do:

VMax = V(:,RB);

VMax will contain only those columns of V that formed the column space, and hence those that are linearly independent vectors. If you want to determine how many independent vectors we have, you simply have to count how many values of RB we have:

r = numel(RB);

---

Here's a quick example to illustrate my point. Supposing I have this matrix:

>> V = [1 1 2 0; 2 2 4 9; 3 3 6 7; 4 4 8 3]

V =

     1     1     2     0
     2     2     4     9
     3     3     6     7
     4     4     8     3

The second column / vector is simply the first vector. The third column / vector is simply the first vector plus the second vector, or it could be twice the first or second vector. Either way, this is not a linearly independent vector as this is based off of the first vector. The same goes with the second vector. The last vector is independent of the other three as there is no way we can create combinations or a scaling that can produce this last vector from the other three. If we called rref, this is what will happen:

>> [R, RB] = rref(V)

R =

     1     1     2     0
     0     0     0     1
     0     0     0     0
     0     0     0     0


RB =

     1     4

R contains the row reduced echelon form, while RB tells us which columns are linearly independent or form the column space of A. As you can see, only columns 1 and 4 are linearly independent, which makes perfect sense. If you also take a look at the last two rows of R, we can see that they consist of all zeroes. This alludes to the rank of a matrix, where it is simply the total number of rows (or columns depending on what you're doing) that are non-zero. This also tells you how many vectors form the column space.

To complete your task, simply do:

>> VMax = V(:,RB)

VMax =

     1     0
     2     9
     3     7
     4     3

As you can see, each column of VMax is a linearly independent vector from V, which also forms the column space of V.

Now, your next task is to randomly choose linearly independent vectors from this column space each time you run the algorithm. Bear in mind that because there are multiple ways to produce the column space, the solution by rref will only give you one such column space. If I am interpreting your question correctly, you wish to generate random column spaces and choose a subset of those vectors each time. Thanks to Luis Mendo (and also a gentle prod from knedlsepp), what you can do is randomly rearrange or permute the columns, then run rref on this permuted matrix.

As such, you can do something like this:

ind = randperm(size(V,2));
Vperm = V(:,ind);
[R,RB] = rref(Vperm);
VMax = Vperm(:,RB);

randperm will generate a random permutation array from 1 up to the argument that you specify to randperm. In our case, this is the total number of columns in your matrix V. We use this array to randomly shuffle the columns of V, store this into Vperm and run our code that we have done before. By doing this random shuffling, you are biasing the input to rref so that you are forcing it to choose different basis vectors, but if you have several vectors that are linearly dependent, there will be a case where we will choose one of these linearly dependent vectors to build our basis.

Answer 2

This naive approach to the problem of searching for a basis of the vector space spanned by your vectors is similar to Luis Mendo's. It should however have better theoretical worst case complexity of a maximum of n rank-computations, as opposed to n choose r. The downside is, that the best case will perform r rank-computations, as opposed to 1.

The idea is simple:

• Compute the rank r of V. This is the maximum number of vectors in your set.
• Define the set of your chosen vectors as empty.
• For each vector V(:,k) in V (randomly chosen):
  • If the vector V(:,k) is linearly independent from the vectors in the set chosen:
    • Add it to the set of chosen vectors.
  • If the set chosen already has r elements:
    • Return the chosen vectors as the solution.

The resulting code will look like this:

function chosen = randBasis(V)
% V is a matrix of column vectors
% chosen are the indices of the randomly selected column vectors,
% so that span(V(:,chosen)) == span(V)

n = size(V,2);
r = rank(V);
chosen = []; 
currentRank = 0;
for k = randperm(n)
    if rank(V(:, [chosen, k])) > currentRank
        currentRank = currentRank+1;
        chosen = [chosen, k];
    end
    if currentRank==r
        break
    end
end

Additional remarks:

This approach can be beneficial when you are dealing with sparse matrices, as the check for linear dependency can be replaced by a solution of a linear system instead of the rank computation. As the rank computation involves a call to svd, it is quite expensive in comparison. Despite rayryeng's rref-solution being the more elegant one, for some reason this naive approach seems to be faster even for non-sparse matrices.

Sped up version using mldivide:

function chosen = randBasis(V)
n = size(V,2);
if issparse(V)
    r = sprank(V)
else
    r = rank(V);
end
chosen = []; 
for k = randperm(n)
    if isLinearlyIndependent(V(:,chosen),V(:,k))
        chosen = [chosen, k];
    end
    if numel(chosen)==r
        break
    end
end
end

function b = isLinearlyIndependent(V, v)
    warning('off','MATLAB:singularMatrix');
    bestapprox = V\v;
    bestapprox(~isfinite(bestapprox)) = 0;
    b = ~(norm(V*bestapprox-v,'inf')<=eps*norm(V,'inf'));
    warning('on','MATLAB:singularMatrix');
end

Beware:

Both these iterative approaches suffers from a non-obvious numerical problem: Consider the matrix:

V = diag(10.^(-40:10:0))
V =
        1e-40            0            0            0            0
            0        1e-30            0            0            0
            0            0        1e-20            0            0
            0            0            0        1e-10            0
            0            0            0            0        1e+00

Mathematically you would say this matrix has rank five. Yet numerically only the last two columns are nonzero compared to the others. This is why rank(V) gives 2. The problem with this iterative approach is that adding a linearly independent vector to our set can in fact result in a linearly dependent set! Imagine this algorithm chooses the first vector. Even though all the other vectors are linearly independent mathematically and numerically, the subset {V(:,1), V(:,5)} does not have numerical rank 2, and yet could be chosen by our function. This is why this algorithm will only work well for vectors of similar norm.

The conclusion is:

• Numerical mathematics is tricky.
• If you want stability, you have to pay the price of computational cost.
• rref is the way to go for this problem.

Answer 3

I would try something along these lines. It will probably not be fast, though:

1. Compute the rank of the matrix formed by your vectors. Call that r.
2. Generate all combinations of r numbers taken from 1:n.
3. Sort those combinations randomly.
4. Take the first combination and check if the vectors indicated by that combination are independent (maximum rank). If they aren't, test next combination. By construction it is guaranteed that some combination will succeed.

Code:

vectors = randi(3,4,3); %// example. Each vector is a row
n = size(vectors,1);
r = rank(vectors); %// step 1
combs = nchoosek(1:n,r); %// step 2
c = size(combs,1);
combs = combs(randperm(c),:); %// step 3
for s = 1:r %// step 4
    pick = combs(s,:);
    if rank(vectors(pick,:))==r
        break %// we're done. Now `pick` contains the selected vectors
    end
end
result = vectors(pick,:); %// each row is a vector

Answer 4

A partial solution on an example ... Sorry, I did not notice your comment about "different linearly independent sets" ... Do you mean "different maximal linearly independent sets"? I apologize I have found only one of the maximal linearly independent sets ... But the solution using [R, RB] = rref(M) gives the same answer

+ echo ('on')
+ more ('off')
+ format ('bank')
+ M = [0, 0, 1; 0, 1, 1; 0, 1, 0]'
M =

   0.00   0.00   0.00
   0.00   1.00   1.00
   1.00   1.00   0.00

+ [U, S, V] = svd (M)
U =

   0.00   0.00  -1.00
  -0.71   0.71   0.00
  -0.71  -0.71   0.00

S =

Diagonal Matrix

   1.73      0      0
      0   1.00      0
      0      0  -0.00

V =

  -0.41  -0.71  -0.58
  -0.82   0.00   0.58
  -0.41   0.71  -0.58

+ U (:, 1:2) * S (1:2, 1:2) * V (1:2, 1:2)'
ans =

   0.00   0.00
   0.00   1.00
   1.00   1.00

A quote from Applications_of_the_SVD

The left-singular vectors corresponding to the non-zero singular values of M span the range of M.