310 lines
7.9 KiB
Plaintext
310 lines
7.9 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b2277afb",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Tutorial 3 - Orthogonalization & QR Decomposition"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "1e4a68da",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"In this tutorial, we study **orthogonalization** methods and the **QR decomposition**, which are central to numerical linear algebra.\n",
|
|
"\n",
|
|
"We will cover:\n",
|
|
"\n",
|
|
"- Orthogonal and orthonormal vectors\n",
|
|
"- QR decomposition\n",
|
|
"- Classical Gram-Schmidt\n",
|
|
"- Modified Gram-Schmidt\n",
|
|
"- Householder transformations\n",
|
|
"- Applications: least squares\n",
|
|
"- Examples with the `numethods` package\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "515dc61f",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 1. Why Orthogonalization?\n",
|
|
"\n",
|
|
"- Orthogonal vectors are easier to work with numerically.\n",
|
|
"- Many algorithms are more **stable** when using orthogonal bases.\n",
|
|
"- Key applications:\n",
|
|
" - Solving **least squares problems**\n",
|
|
" - Computing **eigenvalues**\n",
|
|
" - Ensuring numerical stability in projections\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "af61a18f",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 2. Definitions\n",
|
|
"\n",
|
|
"### Orthogonal and Orthonormal vectors\n",
|
|
"\n",
|
|
"Two vectors $u, v \\in \\mathbb{R}^n$ are **orthogonal** if\n",
|
|
"\n",
|
|
"$$ u \\cdot v = 0. $$\n",
|
|
"\n",
|
|
"A set of vectors $\\{q_1, \\dots, q_m\\}$ is **orthonormal** if\n",
|
|
"\n",
|
|
"$$ q_i \\cdot q_j = \\begin{cases} 1 & i = j, \\\\ 0 & i \\neq j. \\end{cases} $$\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "4f2330f1",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"### QR Decomposition\n",
|
|
"\n",
|
|
"For any $A \\in \\mathbb{R}^{m \\times n}$ with linearly independent columns, we can write\n",
|
|
"\n",
|
|
"$$ A = QR, $$\n",
|
|
"\n",
|
|
"- $Q \\in \\mathbb{R}^{m \\times n}$ has orthonormal columns ($Q^T Q = I$)\n",
|
|
"- $R \\in \\mathbb{R}^{n \\times n}$ is upper triangular\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "66567b83",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 3. Gram-Schmidt Orthogonalization\n",
|
|
"\n",
|
|
"### Classical Gram-Schmidt (CGS)\n",
|
|
"\n",
|
|
"Given linearly independent vectors $a_1, \\dots, a_n$:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"q_1 = \\frac{a_1}{\\|a_1\\|}\n",
|
|
"$$\n",
|
|
"$$\n",
|
|
"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\\|}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Matrix form:\n",
|
|
"\n",
|
|
"$$ A = QR, \\quad R_{jk} = q_j^T a_k. $$\n",
|
|
"\n",
|
|
"⚠️ CGS can lose orthogonality in finite precision arithmetic.\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "3fe97ce3",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Q (CGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])\n",
|
|
"R (CGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"from numethods import Matrix, Vector\n",
|
|
"from numethods import QRGramSchmidt, QRModifiedGramSchmidt, QRHouseholder, LeastSquaresSolver\n",
|
|
"\n",
|
|
"# Example matrix\n",
|
|
"A = Matrix([[2, -1], [1, 2], [1, 1]])\n",
|
|
"\n",
|
|
"# Classical Gram-Schmidt\n",
|
|
"qrg = QRGramSchmidt(A)\n",
|
|
"print(\"Q (CGS):\", qrg.Q)\n",
|
|
"print(\"R (CGS):\", qrg.R)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ba84b59a",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"### Modified Gram-Schmidt (MGS)\n",
|
|
"\n",
|
|
"Same idea, but orthogonalization is done step by step:\n",
|
|
"\n",
|
|
"```\n",
|
|
"for k = 1..n:\n",
|
|
" q_k = a_k\n",
|
|
" for j = 1..k-1:\n",
|
|
" r_jk = q_j^T q_k\n",
|
|
" q_k = q_k - r_jk q_j\n",
|
|
" r_kk = ||q_k||\n",
|
|
" q_k = q_k / r_kk\n",
|
|
"```\n",
|
|
"MGS is more stable than CGS.\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "e01e25ff",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Q (MGS): Matrix([[0.8164965809277261, -0.5520524474738834], [0.4082482904638631, 0.7590721152765896], [0.4082482904638631, 0.34503277967117707]])\n",
|
|
"R (MGS): Matrix([[2.449489742783178, 0.4082482904638631], [0.0, 2.41522945769824]])\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# Modified Gram-Schmidt\n",
|
|
"qrm = QRModifiedGramSchmidt(A)\n",
|
|
"print(\"Q (MGS):\", qrm.Q)\n",
|
|
"print(\"R (MGS):\", qrm.R)\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d893d189",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 4. Householder Reflections\n",
|
|
"\n",
|
|
"A more stable method uses **Householder matrices**.\n",
|
|
"\n",
|
|
"For a vector $x \\in \\mathbb{R}^m$:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"v = x \\pm \\|x\\| e_1, \\quad H = I - 2 \\frac{vv^T}{v^T v}.\n",
|
|
"$$\n",
|
|
"\n",
|
|
"- $H$ is orthogonal ($H^T H = I$).\n",
|
|
"- Applying $H$ zeros out all but the first component of $x$.\n",
|
|
"\n",
|
|
"QR via Householder:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"R = H_n H_{n-1} \\cdots H_1 A, \\quad Q = H_1^T H_2^T \\cdots H_n^T.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "15dfc35c",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Q (Householder): Matrix([[-0.8164965809277258, 0.552052447473883, -0.16903085094570333], [-0.40824829046386296, -0.7590721152765892, -0.5070925528371099], [-0.40824829046386296, -0.34503277967117707, 0.8451542547285166]])\n",
|
|
"R (Householder): Matrix([[-2.449489742783177, -0.408248290463863], [2.220446049250313e-16, -2.415229457698238], [2.220446049250313e-16, 2.220446049250313e-16]])\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# Householder QR\n",
|
|
"qrh = QRHouseholder(A)\n",
|
|
"print(\"Q (Householder):\", qrh.Q)\n",
|
|
"print(\"R (Householder):\", qrh.R)\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "2b6c612f",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 5. Applications of QR\n",
|
|
"\n",
|
|
"### Least Squares\n",
|
|
"\n",
|
|
"We want to solve\n",
|
|
"\n",
|
|
"$$ \\min_x \\|Ax - b\\|_2. $$\n",
|
|
"\n",
|
|
"If $A = QR$, then\n",
|
|
"\n",
|
|
"$$ \\min_x \\|Ax - b\\|_2 = \\min_x \\|QRx - b\\|_2. $$\n",
|
|
"\n",
|
|
"Since $Q$ has orthonormal columns:\n",
|
|
"\n",
|
|
"$$ R x = Q^T b. $$\n",
|
|
"\n",
|
|
"So we can solve using back-substitution.\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "25b399b7",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Least squares solution: Vector([1.0285714285714287, 0.828571428571429])\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# Least squares example\n",
|
|
"b = Vector([1, 2, 3])\n",
|
|
"x_ls = LeastSquaresSolver(A, b).solve()\n",
|
|
"print(\"Least squares solution:\", x_ls)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a94c88c8",
|
|
"metadata": {},
|
|
"source": [
|
|
"\n",
|
|
"## 6. Key Takeaways\n",
|
|
"\n",
|
|
"- CGS is simple but numerically unstable.\n",
|
|
"- MGS is more stable and preferred if using Gram-Schmidt.\n",
|
|
"- Householder QR is the standard in practice (stable and efficient).\n",
|
|
"- QR decomposition underlies least squares, eigenvalue methods, and more.\n"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.13.7"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|