numethods/tutorials/tutorial3_orthogonalization.ipynb

Tutorial 3 - Orthogonalization & QR Decomposition

In this tutorial, we study orthogonalization methods and the QR decomposition, which are central to numerical linear algebra.

We will cover:

• Orthogonal and orthonormal vectors
• QR decomposition
• Classical Gram-Schmidt
• Modified Gram-Schmidt
• Householder transformations
• Applications: least squares
• Examples with the numethods package

1. Why Orthogonalization?

• Orthogonal vectors are easier to work with numerically.
• Many algorithms are more stable when using orthogonal bases.
• Key applications:
  • Solving least squares problems
  • Computing eigenvalues
  • Ensuring numerical stability in projections

2. Definitions

Orthogonal and Orthonormal vectors

Two vectors \(u, v \in \mathbb{R}^n\) are orthogonal if

\[ u \cdot v = 0. \]

A set of vectors \(\{q_1, \dots, q_m\}\) is orthonormal if

\[ q_i \cdot q_j = \begin{cases} 1 & i = j, \\ 0 & i \neq j. \end{cases} \]

QR Decomposition

For any \(A \in \mathbb{R}^{m \times n}\) with linearly independent columns, we can write

\[ A = QR, \]

• \(Q \in \mathbb{R}^{m \times n}\) has orthonormal columns (\(Q^T Q = I\))
• \(R \in \mathbb{R}^{n \times n}\) is upper triangular

3. Gram-Schmidt Orthogonalization

Classical Gram-Schmidt (CGS)

Given linearly independent vectors \(a_1, \dots, a_n\):

\[ q_1 = \frac{a_1}{\|a_1\|} \] \[ q_k = \frac{a_k - \sum_{j=1}^{k-1} (q_j \cdot a_k) q_j}{\left\|a_k - \sum_{j=1}^{k-1} (q_j \cdot a_k) q_j\right\|} \]

Matrix form:

\[ A = QR, \quad R_{jk} = q_j^T a_k. \]

⚠️ CGS can lose orthogonality in finite precision arithmetic.

from numethods import Matrix, Vector
from numethods import QRGramSchmidt, QRModifiedGramSchmidt, QRHouseholder, LeastSquaresSolver

# Example matrix
A = Matrix([[2, -1], [1, 2], [1, 1]])

# Classical Gram-Schmidt
qrg = QRGramSchmidt(A)
print("Q (CGS):", qrg.Q)
print("R (CGS):", qrg.R)

Q (CGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])
R (CGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])

Modified Gram-Schmidt (MGS)

Same idea, but orthogonalization is done step by step:

for k = 1..n:
    q_k = a_k
    for j = 1..k-1:
        r_jk = q_j^T q_k
        q_k = q_k - r_jk q_j
    r_kk = ||q_k||
    q_k = q_k / r_kk

MGS is more stable than CGS.

# Modified Gram-Schmidt
qrm = QRModifiedGramSchmidt(A)
print("Q (MGS):", qrm.Q)
print("R (MGS):", qrm.R)

Q (MGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])
R (MGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])

4. Householder Reflections

A more stable method uses Householder matrices.

For a vector \(x \in \mathbb{R}^m\):

\[ v = x \pm \|x\| e_1, \quad H = I - 2 \frac{vv^T}{v^T v}. \]

• \(H\) is orthogonal (\(H^T H = I\)).
• Applying \(H\) zeros out all but the first component of \(x\).

QR via Householder:

\[ R = H_n H_{n-1} \cdots H_1 A, \quad Q = H_1^T H_2^T \cdots H_n^T. \]

# Householder QR
qrh = QRHouseholder(A)
print("Q (Householder):", qrh.Q)
print("R (Householder):", qrh.R)

Q (Householder): Matrix([[-0.8164965809277258, 0.552052447473883, -0.16903085094570333], [-0.40824829046386296, -0.7590721152765892, -0.5070925528371099], [-0.40824829046386296, -0.34503277967117707, 0.8451542547285166]])
R (Householder): Matrix([[-2.449489742783177, -0.408248290463863], [2.220446049250313e-16, -2.415229457698238], [2.220446049250313e-16, 2.220446049250313e-16]])

5. Applications of QR

Least Squares

We want to solve

\[ \min_x \|Ax - b\|_2. \]

If \(A = QR\), then

\[ \min_x \|Ax - b\|_2 = \min_x \|QRx - b\|_2. \]

Since \(Q\) has orthonormal columns:

\[ R x = Q^T b. \]

So we can solve using back-substitution.

# Least squares example
b = Vector([1, 2, 3])
x_ls = LeastSquaresSolver(A, b).solve()
print("Least squares solution:", x_ls)

Least squares solution: Vector([1.0285714285714287, 0.828571428571429])

6. Key Takeaways

• CGS is simple but numerically unstable.
• MGS is more stable and preferred if using Gram-Schmidt.
• Householder QR is the standard in practice (stable and efficient).
• QR decomposition underlies least squares, eigenvalue methods, and more.