Usually one cares more about identifying linear dependent variables in the model matrix rather than just the rank. This way, if there are linearly dependent variables, these can be dropped such that one obtains a model matrix that is full rank. If the approach to solving the problem is through inversion of the normal matrix, I believe the best method would be.

- Compute the QR decomposition of X.
- The diagonal elements of R identify linearly independent variables as non-zero elements and linearly dependent as those with value zero.
- Then one can compute the Cholesky decomposition of the normal matrix… For calculating the Gram matrix one can use R^{\top}R.

The QR decomposition is a preferred approach to singular-values decomposition (which `rank`

uses) and has can be tailored to have similar tolerances.

Edit:

Actually I was thinking from a perspective of using iterative solvers. If you are using the Cholesky decomposition for solving the problem, then it is better not to do the QR decomposition as it is \mathcal{O}\left(m n^{2}\right) rather than Cholesky which is \mathcal{O}\left(n^{3}\right). The R in the QR decomposition is the same as the upper triangular of the Cholesky factor from the Gram matrix. You can identify which columns are linearly independent checking which diagonal elements of the L are non-zero.

As an aside, I still need to verify if the R approach from the QR decomposition can be linked to the singular-values decomposition for calculating the effective degrees of freedom, (i.e., \text{Tr}\left[ X \left(X^{\top}X + \Gamma\right)^{-1} X^{\top}\right])