diff --git a/lectures/qr_decomp.md b/lectures/qr_decomp.md index f32312d8a..b2cb60c7d 100644 --- a/lectures/qr_decomp.md +++ b/lectures/qr_decomp.md @@ -26,33 +26,32 @@ This lecture describes the QR decomposition and how it relates to We'll write some Python code to help consolidate our understandings. -## Matrix Factorization +## Matrix factorization -The QR decomposition (also called the QR factorization) of a matrix is a decomposition of a matrix into the product of an orthogonal matrix and a triangular matrix. +In this lecture, we consider a real matrix $A \in \mathbb{R}^{n \times m}$ with $m \geq n$. -A QR decomposition of a real matrix $A$ -takes the form +A QR decomposition (also called a QR factorization) of $A$ takes the form $$ -A=QR +A = QR, \qquad +Q \in \mathbb{R}^{n \times n}, \qquad +R \in \mathbb{R}^{n \times m}, \qquad +Q^\top Q = I_n. $$ -where - -* $Q$ is an orthogonal matrix (so that $Q^TQ = I$) +The matrix $R$ has zeros below its main diagonal. -* $R$ is an upper triangular matrix +Hence, $R$ is upper triangular when $m = n$ and upper trapezoidal when $m > n$. +We'll use a *Gram-Schmidt process* to compute a QR decomposition. -We'll use a **Gram-Schmidt process** to compute a QR decomposition - -Because doing so is so educational, we'll write our own Python code to do the job +Because doing so is so educational, we'll write our own Python code to do the job. ## Gram-Schmidt process -We'll start with a **square** matrix $A$. +We'll start with a *square* matrix $A$. -If a square matrix $A$ is nonsingular, then a $QR$ factorization is unique. +If a square matrix $A$ is nonsingular, then a QR factorization is unique up to signs. We'll deal with a rectangular matrix $A$ later. @@ -60,15 +59,15 @@ Actually, our algorithm will work with a rectangular $A$ that is not square. ### Gram-Schmidt process for square $A$ -Here we apply a Gram-Schmidt process to the **columns** of matrix $A$. +Here we apply a Gram-Schmidt process to the *columns* of matrix $A$. In particular, let $$ -A= \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right] +A= \begin{bmatrix} a_1 & a_2 & \cdots & a_n \end{bmatrix} $$ -Let $|| · ||$ denote the L2 norm. +Let $|| \cdot ||$ denote the L2 norm. The Gram-Schmidt algorithm repeatedly combines the following two steps in a particular order @@ -76,16 +75,16 @@ The Gram-Schmidt algorithm repeatedly combines the following two steps in a par * **orthogonalize** the next vector -To begin, we set $u_1 = a_1$ and then **normalize**: +To begin, we set $u_1 = a_1$ and then *normalize*: $$ u_1=a_1, \ \ \ e_1=\frac{u_1}{||u_1||} $$ -We **orthogonalize** first to compute $u_2$ and then **normalize** to create $e_2$: +We *orthogonalize* first to compute $u_2$ and then *normalize* to create $e_2$: $$ -u_2=a_2-(a_2· e_1)e_1, \ \ \ e_2=\frac{u_2}{||u_2||} +u_2=a_2-(a_2\cdot e_1)e_1, \ \ \ e_2=\frac{u_2}{||u_2||} $$ We invite the reader to verify that $e_1$ is orthogonal to $e_2$ by checking that @@ -96,28 +95,28 @@ The Gram-Schmidt procedure continues iterating. Thus, for $k= 2, \ldots, n-1$ we construct $$ -u_{k+1}=a_{k+1}-(a_{k+1}· e_1)e_1-\cdots-(a_{k+1}· e_k)e_k, \ \ \ e_{k+1}=\frac{u_{k+1}}{||u_{k+1}||} +u_{k+1}=a_{k+1}-(a_{k+1}\cdot e_1)e_1-\cdots-(a_{k+1}\cdot e_k)e_k, \ \ \ e_{k+1}=\frac{u_{k+1}}{||u_{k+1}||} $$ -Here $(a_j \cdot e_i)$ can be interpreted as the linear least squares **regression coefficient** of $a_j$ on $e_i$ +Here $(a_j \cdot e_i)$ can be interpreted as the linear least squares *regression coefficient* of $a_j$ on $e_i$ * it is the inner product of $a_j$ and $e_i$ divided by the inner product of $e_i$ where $e_i \cdot e_i = 1$, as *normalization* has assured us. -* this regression coefficient has an interpretation as being a **covariance** divided by a **variance** +* this regression coefficient has an interpretation as being a *covariance* divided by a *variance* It can be verified that $$ -A= \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right]= -\left[ \begin{array}{c|c|c|c} e_1 & e_2 & \cdots & e_n \end{array} \right] -\left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1\\ 0 & a_2·e_2 & \cdots & a_n·e_2 -\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n·e_n \end{matrix} \right] +A= \begin{bmatrix} a_1 & a_2 & \cdots & a_n \end{bmatrix}= +\begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} +\begin{bmatrix} a_1\cdot e_1 & a_2\cdot e_1 & \cdots & a_n\cdot e_1\\ 0 & a_2\cdot e_2 & \cdots & a_n\cdot e_2 +\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n\cdot e_n \end{bmatrix} $$ -Thus, we have constructed the decomposision +Thus, we have constructed the decomposition $$ A = Q R @@ -126,28 +125,27 @@ $$ where $$ -Q = \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right]= -\left[ \begin{array}{c|c|c|c} e_1 & e_2 & \cdots & e_n \end{array} \right] +Q = \begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} $$ and $$ -R = \left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1\\ 0 & a_2·e_2 & \cdots & a_n·e_2 -\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n·e_n \end{matrix} \right] +R = \begin{bmatrix} a_1\cdot e_1 & a_2\cdot e_1 & \cdots & a_n\cdot e_1\\ 0 & a_2\cdot e_2 & \cdots & a_n\cdot e_2 +\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n\cdot e_n \end{bmatrix} $$ ### $A$ not square Now suppose that $A$ is an $n \times m$ matrix where $m > n$. -Then a $QR$ decomposition is +Then a QR decomposition is $$ -A= \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_m \end{array} \right]=\left[ \begin{array}{c|c|c|c} e_1 & e_2 & \cdots & e_n \end{array} \right] -\left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1 & a_{n+1}\cdot e_1 & \cdots & a_{m}\cdot e_1 \\ -0 & a_2·e_2 & \cdots & a_n·e_2 & a_{n+1}\cdot e_2 & \cdots & a_{m}\cdot e_2 \\ \vdots & \vdots & \ddots & \quad \vdots & \vdots & \ddots & \vdots -\\ 0 & 0 & \cdots & a_n·e_n & a_{n+1}\cdot e_n & \cdots & a_{m}\cdot e_n \end{matrix} \right] +A= \begin{bmatrix} a_1 & a_2 & \cdots & a_m \end{bmatrix}=\begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} +\begin{bmatrix} a_1\cdot e_1 & a_2\cdot e_1 & \cdots & a_n\cdot e_1 & a_{n+1}\cdot e_1 & \cdots & a_{m}\cdot e_1 \\ +0 & a_2\cdot e_2 & \cdots & a_n\cdot e_2 & a_{n+1}\cdot e_2 & \cdots & a_{m}\cdot e_2 \\ \vdots & \vdots & \ddots & \quad \vdots & \vdots & \ddots & \vdots +\\ 0 & 0 & \cdots & a_n\cdot e_n & a_{n+1}\cdot e_n & \cdots & a_{m}\cdot e_n \end{bmatrix} $$ which implies that @@ -162,7 +160,7 @@ a_{n+1} & = (a_{n+1}\cdot e_1) e_1 + (a_{n+1}\cdot e_2) e_2 + \cdots + (a_{n+1}\ a_m & = (a_m\cdot e_1) e_1 + (a_m\cdot e_2) e_2 + \cdots + (a_m \cdot e_n) e_n \cr \end{align*} -## Some Code +## Some code Now let's write some homemade Python code to implement a QR decomposition by deploying the Gram-Schmidt process described above. @@ -298,12 +296,12 @@ Q_scipy, R_scipy = adjust_sign(*qr(A)) Q_scipy, R_scipy ``` -## Using QR Decomposition to Compute Eigenvalues +## Using QR decomposition to compute eigenvalues Now for a useful fact about the QR algorithm. -The following iterations on the QR decomposition can be used to compute **eigenvalues** -of a **square** matrix $A$. +The following iterations on the QR decomposition can be used to compute *eigenvalues* +of a *square* matrix $A$. Here is the algorithm: @@ -311,7 +309,7 @@ Here is the algorithm: 2. Form $A_1 = R_0 Q_0 $ . Note that $A_1$ is similar to $A_0$ (easy to verify) and so has the same eigenvalues. -3. Form $A_1 = Q_1 R_1$ (i.e., form the $QR$ decomposition of $A_1$). +3. Form $A_1 = Q_1 R_1$ (i.e., form the QR decomposition of $A_1$). 4. Form $ A_2 = R_1 Q_1 $ and then $A_2 = Q_2 R_2$ . @@ -369,9 +367,9 @@ Compare with the `scipy` package. sorted(np.linalg.eigvals(A)) ``` -## $QR$ and PCA +## QR and PCA -There are interesting connections between the $QR$ decomposition and principal components analysis (PCA). +There are interesting connections between the QR decomposition and principal components analysis (PCA). Here are some. @@ -466,4 +464,4 @@ QPΛPQ = Q @ P_tilde @ Λ @ P_tilde.T @ Q.T ```{code-cell} ipython3 np.abs(QPΛPQ - XX).max() -``` \ No newline at end of file +```