I am trying to write a function that takes either a square or non-square matrix as input and returns a subset of rows such that the returned matrix is full rank. I also need to track which linearly dependent rows were removed. Despite my simple example case, the matrix is not always binary valued.
There seem to be a lot of threads on this topic (e.g. How to find linearly independent rows from a matrix), which led me to trying a QR decomposition approach:
def reduce_to_full_rank(A):
Q, R = np.linalg.qr(A)
to_keep = (np.abs(np.diag(R)) > 1e-10)
return A[to_keep]
What I don't understand is why this approach seems to depend upon the order of the rows in A.
For example, given the input:
A = np.array(
[
[1,1,0],
[1,1,0],
[1,1,1],
]
)
The following code prints a rank of 2, which makes sense because [1,1,0] and [1,1,1] are independent.
Q, R = np.linalg.qr(A)
rank = (np.abs(np.diag(R)) > 1e-10).sum()
print(rank)
However, if I just change the order of A to:
A = np.array(
[
[1,1,1],
[1,1,0],
[1,1,0],
]
)
the exact same code prints a rank of 1 because there's only one non-zero diagonal element of R.
Please can anyone help me understand:
- Why does the number of non-zero diagonal elements in R depend upon the order of rows in A?
- If this is some sort of fundamental limitation of QR decomposition, how else can I solve this problem?
Thank you!
It probably depends on the hardware, indeed.
So, in my case
It is a correct Q/R decomposition, since, Q is orthogonal
is basically identity. So Q is indeed orthogonal.
And R is triangular.
And
is basically the original matrix.
But your Q/R is also correct (with bigger errors, but that is probably because of the copy&paste of numpy's output, which is truncated)
Likewise
is close to identity (again, the 1e-10 instead of 1e-17 is because of copy&paste in the forum)
So your Q is orthogonal as well
Your R is also triangular. And
So both are correct QR decomposition of the same matrix. No doubt about that. And that is math, not numpy (even tho I verified that with numpy, but at least, it is all the same numpy that says that both are OK. Plus, I am certain that any other linear algebra, or even just paper and pen would say the same)
So if two correct QR decomposition of the same matrix haven't the same number of 0 elements in the diagonal, that means that the postulate I overlooked at first (and assumed was true without thinking too much about it, and since it came from math.stackexchange...) is false. It is not a numpy problem. Just, if 2 correct QR decomposition of same rank matrices (and even same matrices) have not the same number of 0 elements in diagonal, it is proof that number of 0 elements in diagonal is not a criteria.
And if you think about it, it is quite logical.
What you known is that Q is made of independent rows and columns (at least in square case). And R columns are how to combine Q columns to create A=QR columns.
So, in my case, first column of A is just (to a multiplicative factor) the first column of Q. And same for your QR decomposition.
Second column of A is the same, since 2nd column of Q is the same as first. For both our decompositions. And that 0 in the diagonal element is indeed proof that we won't be using "new" vectors to create the 2nd column of A; only vectors that we have already used for the 1st. So, rank goes down -1. In both our cases.
Where scenario differs (and where it becomes important to think of the interpretation) is for the 3rd column.
In my case, 3rd column of R has 3 non-zero values (hence a non-0 on the diagonal element, but it doesn't really matter). So, that means that I will be using for the 3rd column of A "new" vectors, that I haven't used for the 2 first. Because it is the first time I have a non-zero in 2nd and 3rd rows. So the first time I will combine 2nd and 3rd columns of Q. So, that is a new dimension. No rank decrease.
In your case, 3rd column of R has only 2 non-zero values. And diagonal element is 0. But it doesn't matter. What matters is that, for the first time, 2nd row is non-zero. You are using to create the 3rd column of A the 2nd column of Q, and it is the first time. You haven't done that for the 2 first columns. So, that is a new dimension as well. So, no decrease of rank, despite the 0 on the diagonal value.
So, what matters is not the diagonal value being non-0. It is a new row being non 0. A row that was 0 in all previous columns being non 0.
In a full rank case, that is on the diagonal value, indeed. But in non full rank cases, that "new non-0" may be on the diagonal, or may be higher.
So, the only case where diagonal is pertinent, is to prove if matrix is full rank or not. If there are 0 on the diagonal of R, then matrix is not full rank. Otherwise it is.
So, short answer: number of non-zero elements in the diagonal of R is not the rank.