Stability of Householder triangularization

from numpy import triu,abs
from numpy.linalg import qr, norm
from numpy.random import normal

Construct a random upper-triangular \(R\), a random unitary \(Q\), form \(A=Q R\) and compute its QR decomposition.

n = 50    # n x n matrix
R = triu(normal(0,1,(n,n)))
Q, X = qr(normal(0,1,(n,n)))
A = Q@R
Q2, R2 = qr(A)

How accurate are the computed QR factors ?

print(norm(Q2 - Q))
print(abs(Q2-Q).max())

1.332256804201302
0.23351568658096772

print(norm(R2 - R) / norm(R))
print(abs(R2-R).max() / abs(R).max())

0.04063321315995649
0.25101293318026785

Very inaccurate !!! But the product

print(norm(A - Q2@R2) / norm(A))
print(abs(A - Q2@R2).max() / abs(A).max())

6.735162924623873e-16
7.861716015013622e-16

is accurate to machine precision !!! The errors in the factors seem to be very special and cancel out in the product.

If we form factors which are equally or more accurate but have random errors

Q3 = Q + 1e-5 * normal(0, 1, (n, n))
R3 = R + 1e-5 * normal(0, 1, (n, n))
print(norm(A - Q3@R3) / norm(A))
print(abs(A - Q3@R3).max() / abs(A).max())

7.364730264449354e-05
6.77367981403941e-05

they dont provide such a good approximation.