Function that reorders the columns of a matrix

Question

The objective: Given an m by n matrix, I have written a function in Matlab that reorders the columns of the matrix to output linearly independent columns; other/redundant columns. Basically, it just pushes the redundant columns to the end of the matrix.

Here is my attempt at the code (which produces correct results):

function result = find_redundant_assets(A)

if (rank(A) == size(A, 2)) % if the rank is the same as the number of columns => full rank
    result = A;
else % there exists at least one linearly dependent column 
    A_red = [];
    for n = 1:(size(A,2)-1) % fix a column and compare to all neighbors to the right 
        for p = n+1:size(A,2) % loop through comparison columns 

            if ( isequal( A(:,p)' * round(A(:, n)'/ A(:,p)') == A(:,n)', ones(1, size(A,1)))) % TRUE iff the nth column is a scalar multiple of the pth column 
                A_red = [A_red, A(:,n), A(:,p)]; % add the linearly depedent rows
            end

            if (isequal(A(:,n)' * round(A(:, p)'/ A(:,n)') == A(:,p)', ones(1, size(A,1)))) % TRUE iff the pth column is a scalar multiple of the nth column
                A_red = [A_red, A(:,p), A(:,n)]; % add the linearly depedent row 
            end

        end
    end
   A_red = unique(A_red.', 'rows', 'stable')';
   result = fliplr([A, A_red]); 
   result = unique(result.', 'rows', 'stable')';
   result = fliplr(result);
end
end

This function gets pretty complicated, but it makes sense to me.

How can I improve my code? Does anyone know a better way to accomplish this task?

Also, I know this is not efficient code because I am constantly re-sizing the matrix.

Answer 1

I don't have Matlab here, and Octave's unique doesn't have the same options as Matlab's, so I can't run the entire code.

So I'll just point out a few things:

isequal( A(:,p)' * round(A(:, n)'/ A(:,p)') == A(:,n)', ones(1, size(A,1)));

Can be written:

all( A(:,p).' * round(A(:,n).'/A(:,p).') == A(:,n).');

You don't need to compare it to a vector of ones, you can just check if all entries are true. Also, conventiently, if [a, b, c, d] will only evaluate to true if all elements a, b, c, d are true. So, the if statement can simply be written:

if (A(:,p).' * round(A(:,n).'/A(:,p).') == A(:,n).');

Notice that I've added a dot . in front of all the '. That's because ' is the complex conjugate transpose, not the regular one.

You have initialized A_red = [], instead of giving it a size, most likely since you don't know the size. Luckily, you do know other things. For instance the maximum possible size, and the any zero-columns can be deleted in the end.

Initialize A_red like this (select an appropriate num_rows:

A_red = zeros(size(A,1), num_rows);`

And do this in the end to remove all the zero-columns:

result = result(:, all(result));

---

The code can probably be written a lot faster, cleaner and smoother. I'll leave that to someone else for now. This should at least improve the code =)