What is a numerically stable way to remove linearly dependent rows/columns from a large matrix? I wish to remove linearly dependent rows/columns of fairly large, real, non-square matrices (around 500x500). Row echelon reduction of course can do this for me with Gaussian elimination, but is there a more numerically stable method available?
 A: In this context a useful tool is the rank-revealing QR factorization of the matrix $A$. It will allow you to determine how far $A$ is from a matrix of rank $r$.
Let $A \in \mathbb{R}^{m \times n}$ be any real matrix with no restrictions on the shape of $A$ (tall/square/wide). Then there exists an orthogonal matrix $Q \in \mathbb{R}^{m \times m}$, a permutations matrix $\Pi$ and an upper triangular matrix $R \in \mathbb{R}^{m \times n}$ such that
$$A\Pi = QR.$$ Moreover, the diagonal entries of $R$ are nonnegative and nonincreasing. Let us partition $Q$ and $R$ conformally, and write
$$ A \Pi = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_{11} & R_{12} \\ 0 & R_{22} \end{bmatrix} = \begin{bmatrix} Q_1 R_{11} & Q_1 R_{12} + Q_2 R_{22} \end{bmatrix}$$
Now if say $\|R_{22}\|_F$ is insignificant compared with the smallest singular value of $R_{11}$, then $A\Pi$ is exceedingly close to a matrix whose rank equals the dimension $r$ of $R_{11}$. In fact, the distance to a matrix of rank $r$ is at most $$\|R_{22}\|_2 \leq \|R_{22}\|_F.$$ In this case, you may wish to retain the columns of $A\Pi$ that are given by $Q_1R_{11}$ and discard the rest.
The use of orthogonal transformations in any kind of serious numerical work is preferred over any other kind of transformations simply because they preserve the 2-norm and the Frobenius norm.
Here is a MATLAB implementation of a rank revealing QR factorization based on Givens rotations. I include the dependency and a minimal working example.
function [Q, R, pi]=ccGivensQRCP(A)
% ccGivensQR  Computes QR factorization with column pivoting
% 
% CALL SEQUENCE: [Q, R, pi]=ccGivensQRCP(A)
%
% INPUT:
%    A      a real matrix m by n matrix
%
% OUTPUT:
%    Q, R   real matrices such that A(:,pi) = QR where
%             Q  is an m by m orthogonal matrix,
%             R  is an m by n upper triangular matrix
%             pi is a permutation of length m
%
% Moreover, the diagonal elemeents of R are nonincreasing
% 
% MINIMAL WORKING EXAMPLE: ccGivensQRCP_mwe1

% PROGRAMMING by Carl Christian Kjelgaard Mikkelsen (spock@cs.umu.se)
%   2022-10-04  Code adapted from ccGivensQR

% Extract the dimension of the matrix
[m, n]=size(A);

% Initialize the output
Q=eye(m,m); R=A; pi=1:n;

% Allocate space for the column norms of A
colnorm=zeros(1,n);
% Compute the norm of the n columns of A
for j=1:n
    colnorm(j)=norm(R(:,j),2);
end

% Loop over the columns of R
for j=1:n
    % Find the column that has the largest norm
    idx=j;
    for i=j+1:n
        if colnorm(i)>colnorm(idx)
            idx=i;
        end
    end
    % At this point we know that colnorm(i)=max(colnorm)
    
    if idx>j 
        % Swap columns j and idx     
        R(:,[j idx])=R(:,[idx j]);
        % Swap the column norms
        colnorm([j idx])=colnorm([idx j]); 
        % Register the swap as well
        pi([j idx])=pi([idx j]);
    end
    % At this point, we know that the dominant column is leading    
        
    % Process the subdiagonal entries of R(:,j) from the bottom up
    % from the bottom up
    for i=m-1:-1:j
        % Givens rotation needed to nullify R(i+1,j) using R(i,j)
        [c, s, r]=ccGivens(R(i,j),R(i+1,j));
        % Define the matrix that represents Givens rotation 
        G=[c -s; s c]; 
        % Apply the rotation to rows i:i+1 ignoring leading zeros
        R(i:i+1,j+1:end)=G'*R(i:i+1,j+1:end);
        % Explicitly set the entries of R(i:i+1,j)
        R(i:i+1,j)=[r; 0];
        % Accumulate the transformation into Q
        Q(:,i:i+1)=Q(:,i:i+1)*G;
    end
    
    % Compute the norms of the columns of the matrix R(j+1:end,j+1:end)
    % Certainly, this can be done faster, but the slow procedure does 
    % not invite subtractive cancellation.
    for k=j+1:n
        colnorm(k)=norm(R(j+1:end,k),2);
    end
end

The dependency:
function [c, s, r]=ccGivens(x,y)

% ccGivens  Computes a Givens rotation
%
% CALL SEQUENCE: [c, s, r]=ccGivens(x,y)
%
% INPUT:
%   x    real number
%   y    real number
%
% OUTPUT:
%   c,s  the cosine and sine such that
%          [c -s; s c]'*[x; y] = [r; 0]
%   r    the norm of the vector [x; y]
%
% MINIMAL WORKING EXAMPLE: ccGivens_mwe1

% PROGRAMMING by Carl Christian Kjelgaard Mikkelsen (spock@cs.umu.se)
%   2020-08-20  Initial programming and testing

% Compute the norm of the vector [x,y]
r=norm([x y]);

% Is this a special case?
if abs(x)~=0
    % No, apply the standard formula
    c=x/r; s=y/r;
else
    % Yes, rotate 90 degrees
    c=0; s=1;
end

A minimal working example:
% Compute a QR factorization with column pivoting.

% Set the dimension of the tall matrix
m=7; n=6; 

% Select the rank r<=n;
r=3; k=n-r;

% Set the random seed
rng(2022);

% //////////////////////////////////////////////////////////
% Generate a random matrix of rank r
% //////////////////////////////////////////////////////////

% Almost certainly A will have r linearly independent columns 
A=rand(m,r); 
% Generate a random matrix and extend A to n columns
B=rand(k,k); A=[A A*B]; 
% The extra columns are linear combinations of the previous columns
% so we are not increasing the rank

% Shuffle the columns of A randomly
p=randperm(n); A=A(:,p);

% At this point A is almost certainly an m by n matrix of rank r.

% Compute the QR factorization
[Q, R, pi]=ccGivensQRCP(A);

% /////////////////////////////////////
% Verify the factorization
% /////////////////////////////////////

% Compute the residual S
S=A(:,pi)-Q*R; 

% Compute the relative residual
rres=norm(S)/norm(A);

% Check orthogonality
orth=norm(eye(m,m)-Q'*Q,2);

% Display results
fprintf('Relative residual = %8.4e\n',rres);
fprintf('Orthogonality     = %8.4e\n',orth);

% Change the format so that tiny nonzeros are emphasized
format short e

% Display the structure of R
display(R);
% Notice that the diagonal entries decay
% Notice that R has a bunch of tiny entries
% Notice that it is clear that R is close to a matrix of rank r.

A: A numerically stable way to remove linearly dependent rows/columns from a large matrix is to use a QR factorization.
