Say we're given a set of $d$ vectors $S=\{\mathbf{v}_1,\dots,\mathbf{v}_d\}$ in $\mathbb{R}^n$, with $d\leq n$ (obviously). We want to test in an efficient way if S is linearly independent. Now, write the coefficient matrix $\mathbf{A}=[\mathbf{v}_1 \dots\mathbf{v}_d]$ (the $\mathbf{v}_i$ are considered to be column vectors).

A non-efficient way would be to compute all minors of rank $d$ of $\mathbf{A}$, and check that they are non-zero (up to some tolerance, as always when we do numerical linear algebra). Another way would be using Gram–Schmidt orthogonalization, but I recall that Gram–Schmidt orthogonalization is numerically unstable. Which is the correct alternative? Singular Value Decomposition (if all singular values are strictly positive, then the vectors are independent)? QR factorization?

  • 1
    $\begingroup$ To the extent that this is a numerical algebra question, and the terms "robust" and "efficient" are not defined, this is not a well-posed question. In your situation I would try to find the Cholesky decomposition of the $d\times d$ Gram matrix of your vectors. $\endgroup$ Dec 17, 2019 at 18:57
  • $\begingroup$ If you want stability, use SVD. That's what Matlab and other general-purpose tools do. $\endgroup$ Dec 19, 2019 at 17:11
  • $\begingroup$ @AlgebraicPavel they use QR factorization, usually $\endgroup$
    – DeltaIV
    Dec 20, 2019 at 15:35
  • 1
    $\begingroup$ @DeltaIV Sure, some also use LU. However, SVD is currently (AFAIK) the most robust and recommended. It is used by Matlab, armadillo, numpy,... just to name a few. $\endgroup$ Dec 20, 2019 at 20:10

2 Answers 2


Let $A$ be not the $n\times d$-matrix, but the $n\times n$-matrix $$A=[v_1 \ldots v_d\ \vec{0}\ldots\ \vec 0]^t.$$ Reduce $A$ to Smith Normal Form. This takes $O(n\cdot d)$ row operations. The number of non-zero rows of the resulting matrix equals the number of these vectors that were linearly independent.

  • 1
    $\begingroup$ How do you compute the Jordan Form in $O(n\cdot d)$ steps? Computing the eigenvalues alone seems harder. $\endgroup$ Dec 17, 2019 at 19:00
  • 2
    $\begingroup$ Explaining how would take a while, so I changed my answer to say "Smith Normal Form" instead of "Jordan Canonical Form". Its reduction also takes $O(n\cdot d)$ multiplications, and gives the same result -- but in this case, WikiPedia actually gives an example of how to do it. Please, take a look at en.wikipedia.org/wiki/Smith_normal_form#Example for how they did it. Thanks! $\endgroup$
    – Alex
    Dec 17, 2019 at 19:05
  • $\begingroup$ $O(nd)$ complexity? So you say that finding a matrix rank has essentially the same complexity as just writing it down? Interesting, please tell me more about it! $\endgroup$ Dec 18, 2019 at 19:09
  • $\begingroup$ @AlgebraicPavel Look at the algorithm in that WikiPedia article. I believe it's pretty clear how it works... If you can't figure it out, you can ask a separate question and link it here, and I can explain it. $\endgroup$
    – Alex
    Dec 19, 2019 at 8:06
  • $\begingroup$ I accepted the answer, but I have to say that the number seems weird to me too. I used to use Gaussian elimination or Gauss-Jordan to compute the Jordan form, and that's $O(n^3)$ for a square matrix. According to your formula, it would be $O(n^2)$. $\endgroup$
    – DeltaIV
    Dec 19, 2019 at 9:37

The standard numerical approaches would be computing a (rank revealing) QR decomposition or SVD. Both Numpy and Matlab use the SVD, and both algorithms are $nd^2$[1]. In general, whenever stability is a concern, the SVD tends to be preferred. However, even plain QR with householder reflections or modified Gram-Schmidt may give reasonable results in many cases.

As you allude to in the original question, there is always a trade off between accuracy and computation time, and the "best" algorithm depends on your application. The notion of "rank" itself is not well defined in finite precision, and in defining numerical rank it would be reasonable to take any of the equivalent characterizations of exact rank, and then relax them.

In general, one could define $A$ as numerically rank $k$ if all singular values $k+1, \ldots, n$ are less than some tolerance (this is what numpy and matlab do). Of course, even if you define it this way, how close your numerical method comes to computing the exact SVD of your finite precision matrix is another concern; i.e. two different implementations of an SVD algorithm would give different results.

On the other hand, if $1\ll k\ll n$ and you are satisfied with an approximation to the rank, then perhaps a randomized algorithm would be best.

I am not familiar with Smith normal form, but I would be very hesitant to use it in a general purpose numerical method as it seems like it would be susceptible to the same issues as Gaussian elimination.

[1] a rank-$r$ SVD can be computed somewhat more cheaply, so if you have an upper bound on the rank of your matrix then you may be able to do things somewhat faster.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .