[BUG] linalg.qr function seems to be wrong

[vega1013]

Describe the bug The qr functions seems to calculate wrong values if i'm not totally wrong. The proof "np.dot(q,r)" == A is true, but when i calculate the eigenvectors based on R they differ (its not a problem of scaling the eigenvectors) with numpy/Matlab.

If any additional information is needed please let me know.

ulab version: 4.4.0-2D-c

from ulab import numpy as np
A = np.array([[1, 1, 0], [0, 0, 2], [0, 0, -1]])
Q,E = np.linalg.qr(A)
>>> Q
array([[1.0, 0.0, 0.0],
       [0.0, -1.0, -0.0],
       [0.0, 0.0, 1.0]], dtype=float64)
>>> R
array([[1.0, 1.0, 0.0],
       [0.0, 0.0, -2.0],
       [0.0, 0.0, -1.0]], dtype=float64)

To Reproduce See the code example above.

Expected behavior In numpy the resulting q and r matrix are as follow:

>>> Q
array([[ 1.,  0.,  0.],
       [-0.,  1.,  0.],
       [-0., -0.,  1.]])
>>> R
array([[ 1.,  1.,  0.],
       [ 0.,  0.,  2.],
       [ 0.,  0., -1.]])

This is also the result in Matlab.

[v923z]

Thanks for bringing this up! I will try to sort it out.

[vega1013]

Some additional Info regarding this issue: If I use a QR-algorithm based on givens rotation I get the same results as in numpy/matlab (source is https://www.programmersought.com/article/84605916864/ with some adjustments for ulab):

def givens_rotation(A):
    # (r, c) = np.shape(A)
    (r, c) = A.shape
    # Q = np.identity(r)
    Q = np.eye(r)
    # R = np.copy(A)
    R = A.copy()
    # (rows, cols) = np.tril_indices(r, -1, c)
    (rows, cols) = my_tri_indices(r, -1, c)
    for (row, col) in zip(rows, cols):
        if R[row, col] != 0:  # R[row, col]=0 then c=1, s=0, R and Q remain unchanged
            # r_ = np.hypot(R[col, col], R[row, col])  # d
            r_ = np.sqrt(R[col, col] * R[row, col] + R[col, col] * R[row, col])  # d
            c = R[col, col]/r_
            s = -R[row, col]/r_
            # G = np.identity(r)
            G = np.eye(r)
            G[[col, row], [col, row]] = c
            G[row, col] = s
            G[col, row] = -s
            R = np.dot(G, R)  # R=G(n-1,n)*...*G(2n)*...*G(23,1n)*...*G(12)*A
            Q = np.dot(Q, G.T)  # Q=G(n-1,n).T*...*G(2n).T*...*G(23,1n).T*...*G(12).T
    return (Q, R)
def my_tri_indices(n, k, m):
    trimat = np.zeros((n,m))
    for i in range(n,-k-1,-1):
        trimat = trimat + np.eye(n,M=m,k=-i)
    trimat_indices = index_2D(trimat,1)
    return trimat_indices
def index_2D(A,c):
    k = 0
    a = np.empty(round(sum(sum(A))), dtype=np.uint8)
    b = np.empty(round(sum(sum(A))), dtype=np.uint8)
    for i, e in enumerate(A.tolist()):
        for j, ee in enumerate(e):
            if A[i][j] == c in e:
                a[k] = round(i)
                b[k] = round(j)
                k = k+1
    return (a,b)

The functions can be tested with the following code.

from ulab import numpy as np
A = np.array([[1, 1, 0], [0, 0, 2], [0, 0, -1]])
Q,R = givens_rotation(A)
eigval = np.diag(R)
V = np.eye(A.shape[0])
for i in range(1, A.shape[0] , 1):
    V[0:i, i] = np.dot( np.linalg.inv((R[i, i]* np.eye(i) - R[0:i, 0:i])), R[0:i, i])    
eigvec_QR = np.dot(Q,V)

The eigenvectors are calculated based on an answer in this thread: https://math.stackexchange.com/questions/3947108/how-to-get-eigenvectors-using-qr-algorithm

I hope that this information is helpful for someone.

[jimmo]

@vega1013

The qr functions seems to calculate wrong values if i'm not totally wrong.

QR decomposition isn't unique, and this is a valid QR decomposition. I would guess that perhaps numpy and R are trying to arrange one where Q is positive.

The functions can be tested with the following code.

It doesn't look like this code is computing the Schur decomposition (via the QR algorithm), it's only computing the QR decomposition.

The eigenvectors are calculated based on an answer in this thread:

I implemented the algorithm from that thread, using ulab's QR decomposition to do the Schur decomposition and got the expected eigenvectors for A (and it also works for the example in their matlab code too).

https://gist.github.com/jimmo/8e98661886f9c7785099871c5cd92871

[jimmo]

Note of course that to compute the QR algorithm you need a test of convergence rather than just doing 100 iterations :)

[v923z]

@vega1013

The qr functions seems to calculate wrong values if i'm not totally wrong.

QR decomposition isn't unique, and this is a valid QR decomposition. I would guess that perhaps numpy and R are trying to arrange one where Q is positive.

This has been bugging me for some time, because, while the ulab tests based on multiplying Q and R with the appropriate matrix always came back positive, I didn't always get the same Q matrix as numpy. So, you are, in effect, saying that numpy flips the sign of the matrix entries, if they are negative, aren't you?

Albeit this sign bears no relevance on the validity of the output, we could definitely fix this, if people think that aesthetic compatibility with numpy is important.

[jimmo]

So, you are, in effect, saying that numpy flips the sign of the matrix entries, if they are negative, aren't you?

I'm not sure if it's doing it intentionally, it could just be a result of whatever algorithm they use for the QR decomposition on this particular matrix.

Albeit this sign bears no relevance on the validity of the output, we could definitely fix this, if people think that aesthetic compatibility with numpy is important.

FWIW, Mathematica doesn't even give you square matrices for Q and R in this case :)

QRDecomposition[{{1, 1, 0}, {0, 0, 2}, {0, 0, -1}}]

{
  {{1, 0, 0}, {0, 2/Sqrt[5], -(1/Sqrt[5])}},
  {{1, 1, 0}, {0, 0, Sqrt[5]}}
}

[v923z]

Since there hasn't been any traffic on this thread for over a year now, and the issue is really just cosmetic, closing it now.