In?[1]:

import numpy as np

a1 = np.array([3, -6, 2])
a2 = np.array([6, -6, 1])
a3 = np.array([-1, 1, -1])

A = np.column_stack((a1, a2, a3))

In?[2]:

# Gram-Schmidt procedure
u1 = a1
u2 = a2 - (np.dot(u1, a2) / np.dot(u1, u1)) * u1
u3 = (
    a3 - (np.dot(u1, a3) / np.dot(u1, u1)) * u1 - (np.dot(u2, a3) / np.dot(u2, u2)) * u2
)

In?[3]:

# Noramlization of orthogonal basis
q1 = u1 / np.sqrt(np.dot(u1, u1))
q2 = u2 / np.sqrt(np.dot(u2, u2))
q3 = u3 / np.sqrt(np.dot(u3, u3))

# Matrix Q
Q = np.column_stack((q1, q2, q3))

In?[4]:

# Matrix R
R = np.zeros((3, 3))
for i in range(0, 3):
    for j in range(i, 3):
        R[i, j] = np.dot(Q[:, i], A[:, j])

In?[5]:

print(np.linalg.norm(Q@R - A, ord=2))

1.7922524720510307e-15

In?[6]:

print(Q)

[[ 0.42857143  0.85714286 -0.28571429]
 [-0.85714286  0.28571429 -0.42857143]
 [ 0.28571429 -0.42857143 -0.85714286]]

In?[7]:

print(R)

[[ 7.          8.         -1.57142857]
 [ 0.          3.         -0.14285714]
 [ 0.          0.          0.71428571]]

In?[?]: