def lstqr(A, b):
    # Initialize
    n = len(b)
    x = np.zeros(n)
    r = np.zeros(n)
    p = np.zeros(n)
    q = np.zeros(n)
    r[0] = np.linalg.norm(b)
    p[0] = b[0] / r[0]
    for i in range(1, n):
        s = 0
        for j in range(0, i):
            s += A[i][j] * A[i][j] * r[j]
        r[i] = np.linalg.norm(b[i] - s)
        q[i] = (b[i] - s) / r[i]
    for i in range(n - 1, -1, -1):
        s = 0
        for j in range(i + 1, n):
            s += A[i][j] * A[i][j] * r[j]
        x[i] = (q[i] - s) / r[i]
    return x