Overall task: Solve linear system of equations
In the discussion afterwards, we assume that \(\boldsymbol{A}\) is invertible, \(\boldsymbol{A}\in \mathbb{R} ^{n\times n}\) but \(n\) is very large.
+In scientific computing, there are three important issues:
+Here are two examples for linear systems in scientific computing.
+Consider:
+where \(f(x)\) is given. We want to know \(u\) using the computer.
+The important step is to do the discretization.
+We split the interval \([0,1]\) into: \(0=x_0<x_1<\cdots <x_{i-1}<x_i<x_{i+1}<\cdots <x_N=1\), where \(x_{i+1}-x_i\triangleq h=\frac{1}{N}, \left( i=0,1,\cdots ,N-1 \right)\). We just need to know the solutions on these points for \(u(x_i)\).
+Let \(u\left( x_i \right) \triangleq u_i\).
+Using the finite difference method, we know that:
+if \(\varepsilon\) is small. We take \(\varepsilon=h\), then \(u\prime\left( x \right) \approx \frac{u\left( x+h \right) -u\left( x \right)}{h}\). We can get \(\frac{u\left( x_i+h \right) -u\left( x_i \right)}{h}=\frac{u_{i+1}-u_i}{h}\).
+Similarly, we can get \(u''\left( x \right) \approx \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}\). Also \(u''\left( x \right) \propto O\left( h^2 \right)\).
+Therefore, \(-\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}=f\left( x_i \right)\)
+To sum up, \(\boldsymbol{u}=\left[ u_1,u_2,\cdots ,u_{N-1} \right]^T\) is a linear system that satisfies:
+We can get:
+The degrading image model is:
+Sometimes \(k\left( x \right) =\frac{1}{c}e^{-\frac{x^2}{2}}\) which is Gaussian blur. Usually:
+To simplify, we ignore the noise: let \(\eta \left( x \right) =0\).
+Given \(g(x)\), we want to find \(u(x)\). Assume that \(k(x)\) is given and the image/signal is 1D. We can model the linear system as:
+\(\boldsymbol{A}\) is called Toeplitz matrix.
+ + + + + + +The range of vector norm is \(\boldsymbol{x}\in \mathbb{R} ^n, \left\| \boldsymbol{x} \right\| \in \mathbb{R} ^+\cup \left\{ \mathbf{0} \right\}\).
+Example:
+\(L^2\)-Norm:
+\(L^1\)-Norm:
+\(L^\infty\)-Norm:
+\(L^p\)-Norm:
+Weighted Norm: ( \(\boldsymbol{W}\) is a given matrix)
+Geometric Intuition:
+Unit disk with different norms: \(\left\{ \boldsymbol{x}\in \mathbb{R} ^n|\left\| \boldsymbol{x} \right\| \le 1 \right\}\)
+\(L^2\) is like a circle. \(L^1\) is like a diamond. \(L^\infty\) is like a square. \(L^p (1<p<2)\) is like the shape between circle and diamond. \(L^p (p>2)\) is like the shape between circle and square.
+Two different ways to define the matrix norm for \(\boldsymbol{A}\in \mathbb{R} ^{m\times n}\):
+Way 1: View matrix as a vector (reshape):
+Way 2: Induced Matrix Norm (preferred way to define the matrix norm)
+Also:
+Example:
+If \(\left\| \boldsymbol{x} \right\| _1\) is taken:
+If \(\left\| \boldsymbol{x} \right\| _{\infty}\) is taken:
+What about induced 2-norms? It turns out that we cannot find a formula for it!
+ + + + + + +Let's start with an example:
+Given \(x\in \mathbb{R} ^+\) , compute \(\sqrt{x}\) . We are considering the problem not the algorithm here.
+We denote the problem as:
+Rather than considering \(\frac{\left\| \delta f \right\|}{\left\| fx \right\|}\) , we can consider the relative ratio:
+For \(f\left( x \right) =\sqrt{x}\) :
+It means that the output does not amplify the permutation/mistake of the input.
+Assume \(\delta x\) is small, we define the condition number as:
+Here \(f\) is any problem.
+++If \(\boldsymbol{A}\in \mathbb{R} ^{n\times n}\) , then matrix condition number can be defined as:
+\[ +\kappa \left( \boldsymbol{A} \right) =\left\| \boldsymbol{A} \right\| \cdot \left\| \boldsymbol{A} \right\| ^{-1} +\]+
Why is this definition important?
+Consider \(\boldsymbol{Ax}=\boldsymbol{b}\), perturbations may occur to \(\boldsymbol{A}, \boldsymbol{x}, \boldsymbol{b}\) .
+\(\boldsymbol{x}\rightarrow \boldsymbol{x}+\delta \boldsymbol{x}\) , \(\boldsymbol{A}\) is given.
+Introduce \(\boldsymbol{y}=\boldsymbol{Ax}\) , \(\boldsymbol{A}\) is invertible, then:
+Perturbation happens to \(\boldsymbol{b}\) :
+Using case 1, we can easily have the condition number:
+Given \(\boldsymbol{A}, \boldsymbol{b}\) , and let \(\boldsymbol{A}\rightarrow \boldsymbol{A}+\delta \boldsymbol{A}\) . Then it incurs \(\boldsymbol{x}\rightarrow \boldsymbol{x}+\delta \boldsymbol{x}\) . We can get:
+Usually \(\delta \boldsymbol{A}\cdot \delta \boldsymbol{x}\) is smaller than \(\delta \boldsymbol{x}\) or \(\delta \boldsymbol{A}\), then we can also get:
+Based on the three conditions above, we can derive the conclusion:
+Here are a few more examples:
+If:
+then \(\kappa \left( \boldsymbol{A} \right) =\frac{\lambda _1}{\lambda _n}\)
+Let's look at another example. If:
+then:
+The \(\kappa \left( \boldsymbol{A} \right)\) satisfies \(\kappa \left( \boldsymbol{A} \right) \geqslant 1\) .
+If \(\kappa \left( \boldsymbol{A} \right)\) is large, the question is called ill-conditioned. If \(\kappa \left( \boldsymbol{A} \right)\) is small, it is called well-conditioned.
+ + + + + + +Stability is for the algorithms.
+Two concepts:
+ +Problem: \(f: X\mapsto Y\)
+Algorithm: \(\tilde{f}: X\mapsto Y\)
+The algorithm \(\tilde{f}\) is used to solve the problem \(f\)
+Definition (Backward Stability): An algorithm \(\tilde{f}\) is called backward stable if for each \(x\in X\) :
for some \(\tilde{x}\) with \(\frac{\left\| x-\tilde{x} \right\|}{\left\| x \right\|}\leqslant O\left( \varepsilon _{\mathrm{machine}} \right)\)
+++In words, a backward stable algorithm gives exactly the right solution to a nearly right problem.
+
Theorem: If \(\tilde{f}\) is a backward stable algorithm for a problem \(f: X\mapsto Y\) with condition number \(\kappa \left( f \right)\) , then the relative error satisfies:
Proof: \(\tilde{f}\) is backward stable: \(\forall x\in X\) , \(\exists \tilde{x}\in X\) such that \(\tilde{f}\left( x \right) =f\left( \tilde{x} \right) , \frac{\left\| x-\tilde{x} \right\|}{\left\| x \right\|}\leqslant O\left( \varepsilon _{\mathrm{machine}} \right)\) . Then we can get:
Remark: If \(\boldsymbol{Ax}=\boldsymbol{b}\) is ill-conditioned, one must expect to lose \(\log _{10}\kappa \left( \boldsymbol{A} \right)\) digits in computing the solution with a backward stable method except under very special cases.
Question: How can we know an algorithm is backward stable or not? (tough and tedious)
Example: Forward substitution is backward stable.
Consider \(\boldsymbol{Rx}=\boldsymbol{b}\) , where \(\boldsymbol{R}\in \mathbb{R} ^{m\times m}\) is a lower triangular matrix. The forward substitution steps are:
+where \(r_{ij}\) is the element of \(\boldsymbol{R}\) .
+Denote \(+,-,\times ,\div\) as exact operations and \(\oplus ,\ominus ,\otimes ,\oslash\) as operations with error.
+Case \(m=1\) : Exact version:
+Computed version:
+Introduce \(\left( 1+\varepsilon _1 \right) \left( 1+\varepsilon _1\prime \right) =1\) . We can rewrite as:
+The outcome can be viewed as an exact solution for a slightly perturbed problem.
+Case \(m=2\) :
+Step 1:
+Step 2:
+Introduce \(\left( 1+\varepsilon _i \right) \left( 1+\varepsilon _i\prime \right) =1\) , then:
+Then:
+We can claim that solving \(\boldsymbol{Rx}=\boldsymbol{b}\) by a computer with a numerical solution \(\hat{\boldsymbol{x}}\) can be viewed as \(\left( \boldsymbol{R}+\delta \boldsymbol{R} \right) \hat{\boldsymbol{x}}=\boldsymbol{b}\) exactly. Also:
+For general \(m\) , \(\hat{\boldsymbol{x}}\) satisfies:
+Finally, we can claim that forward substitution algorithm computes the solution \(\hat{\boldsymbol{x}}\) of \(\boldsymbol{Rx}=\boldsymbol{b}\) satisfying
+for some \(\delta \boldsymbol{R}\) satisfying
+Also, forward and backward substitutions are all backward stable.
+Problem: \(f: X\mapsto Y\)
+Algorithm: \(\tilde{f}: X\mapsto Y\)
+Definition: \(\tilde{f}\) is called stable if for \(\forall x\in X\) ,
for some \(\tilde{x}\) satisfying \(\frac{\left\| \tilde{x}-x \right\|}{\left\| x \right\|}=O\left( \varepsilon _{\mathrm{machine}} \right)\)
+Conclusion: If \(\tilde{f}\) is backward stable, then it is also stable.
++In words, a stable algorithm gives nearly the right solution to nearly the right algorithm.
+
Theorem: Let \(\boldsymbol{A}=\boldsymbol{LU}\) of a nonsingular matrix \(\boldsymbol{A}\in \mathbb{R} ^{m\times m}\) be computed by Gaussian Elimination without pivoting, then:
for some \(\delta \boldsymbol{A}\in \mathbb{R} ^{m\times m}\) satisfying \(\frac{\left\| \delta \boldsymbol{A} \right\|}{\left\| \boldsymbol{L} \right\| \left\| \boldsymbol{U} \right\|}=O\left( \varepsilon _{\mathrm{machine}} \right)\)
+If \(\left\| \boldsymbol{L} \right\| \left\| \boldsymbol{U} \right\| =O\left( \left\| \boldsymbol{A} \right\| \right)\) , then Gaussian Elimination without pivoting is stable. Otherwise, the method is not stable.
+Theorem: Gaussian Elimination with partial pivoting is backward stable:
for some \(\delta \boldsymbol{A}\in \mathbb{R} ^{m\times m}\) satisfying \(\frac{\left\| \delta \boldsymbol{A} \right\|}{\left\| \boldsymbol{A} \right\|}=O\left( \rho \cdot \varepsilon _{\mathrm{machine}} \right)\) where
+is called growth factor.
+Theorem: If \(\boldsymbol{A}\) is symmetric positive definite (SPD), then Cholesky Factorization is backward stable:
for some \(\delta \boldsymbol{A}\) satisfying \(\frac{\left\| \delta \boldsymbol{A} \right\|}{\left\| \boldsymbol{A} \right\|}=O\left( \varepsilon _{\mathrm{machine}} \right)\)
+ + + + + + +Describing the small perturbation properties
+Conditioning: for the problem
+Stability: for the algorithm
+ + + + + + +Eigenvalue problem phase 1: \(\boldsymbol{A}\rightarrow \boldsymbol{H}\,\,\left( \boldsymbol{A}\sim \boldsymbol{H} \right)\), \(\boldsymbol{H}\) is a upper Hessenberg matrix. We get \(\boldsymbol{A}=\boldsymbol{XHX}^{-1}\).
+For example, let \(\boldsymbol{A}\in \mathbb{R} ^{n\times n}\), the steps for phase 1 can be shown as:
+Algorithm: Householder Reduction to upper Hessenberg form:
+For \(k=1:m-2\):
+End
+The cost of the algorithm is \(O\left( \frac{10}{3}m^3 \right)\).
+Why does symmetric positive definite (SPD) matrix is easier to compute? Most stable situation in physics.
+If \(\boldsymbol{A}\) is symmetric, the cost is reduced to \(O\left( \frac{4}{3}m^3 \right)\), and \(\boldsymbol{H}\) is tridiagonal and also symmetric.
+Theorem (backward stability): Let the Hessenberg reduction \(\boldsymbol{A}=\boldsymbol{QHQ}^*\) be computed by the householder algorithm, then we have \(\tilde{\boldsymbol{Q}}\tilde{\boldsymbol{H}}\tilde{\boldsymbol{Q}}^*=\boldsymbol{A}+\delta \boldsymbol{A}\) where for some \(\delta \boldsymbol{A}\) satisfying \(\frac{\left\| \delta \boldsymbol{A} \right\|}{\left\| \boldsymbol{A} \right\|}=O\left( \varepsilon _{\mathrm{machine}} \right)\).
We discuss some classical eigenvalue solvers below. Assume \(\boldsymbol{A}\) is SPD, and eigenvalues \(\lambda _1,\lambda _2,\cdots \lambda\) satisfy \(\lambda _1\geqslant \lambda _2\geqslant \cdots \geqslant \lambda _m\). Their corresponding eigenvectors are \(\boldsymbol{q}_1,\boldsymbol{q}_2,\cdots \boldsymbol{q}_m\).
+Power Iteration algorithm:
+Pick \(\boldsymbol{v}^{\left( 0 \right)}\) as some vector with \(\left\| \boldsymbol{v}^{\left( 0 \right)} \right\| =1\).
+For \(k=1,2,\cdots\):
+End
+Claim of convergence: \(\lambda ^{\left( k \right)}\rightarrow \lambda _1,\boldsymbol{v}^{\left( k \right)}\rightarrow \boldsymbol{q}_1\) as \(k\rightarrow +\infty\).
+Proof: We can write \(\boldsymbol{v}^{\left( k \right)}=c_k\boldsymbol{A}^k\boldsymbol{v}^{\left( 0 \right)}\). we can also get \(\boldsymbol{v}^{\left( 0 \right)}=a_1\boldsymbol{q}_1+a_2\boldsymbol{q}_2+\cdots +a_m\boldsymbol{q}_m\). Therefore,
+We know that \(\left| \frac{\lambda _i}{\lambda _1} \right|<1, i\geqslant 2\). Therefore, \(\left( \frac{{\lambda _i}^k}{{\lambda _1}^k} \right) \rightarrow 0\) as \(k\rightarrow +\infty\). We can get the result.
+ + + + + + +Overall task: Given
Find \(\lambda\) and \(\boldsymbol{x}\).
+What are the methods to compute eigenvalues and eigenvectors?
+Traditional Method:
+Example: Let
+Then \(\lambda _{1,2}=1\), \(p_{\boldsymbol{A}}\left( z \right) =\left( z-1 \right) ^2=z^2-2z+1\). The coefficients encodes in the computer are \(\left[ \begin{matrix} + 1& -2& 1\\ +\end{matrix} \right]\). Assume that the input is \(\left[ \begin{matrix} + 1& -2& 1-10^{-16}\\ +\end{matrix} \right]\). Then the polynomial becomes \(z^2-2z+\left( 1-10^{-16} \right) =\left( z-\left( 1-10^{-8} \right) \right) \cdot \left( z-\left( 1+10^{-8} \right) \right) =0\). We get \(\lambda _1=1-10^{-8}, \lambda _2=1+10^{-8}\). Therefore, when the input error is \(10^{-16}\), the output error is \(O(10^{-8})\). The root finding procedure amplifies the error!
+Root finding is not stable, and algorithms using root finding procedures are not a stable algorithms.
+Actually, root finding in computers is implemented as eigenvalue solving algorithms.
+Claim: Eigenvalue problem is equivalent to root finding for polynomials.
Given \(p\left( z \right) =z^m+a_{m-1}z^{m-1}+\cdots +a_1z+a_0\). Introduce a matrix:
+\(\boldsymbol{A}\in \mathbb{R} ^{m\times m}\) is called the companion matrix for \(p(z)\). We can prove that:
+Then getting the roots for the original polynomial is euqivalent to getting the eigenvalues for the corresponding companion matrix.
+There is fundamental difficulty in the eigenvalue computation: no explicit formula for a general matrix \(\boldsymbol{A}_{m\times m}\) when \(m\geqslant 5\).
+Theorem: For any \(m\geqslant 5\), there is a polynomial \(p(z)\) of degree \(m\) with rational coefficients that has a real root \(p(r)=0\) with the property that \(r\) can NOT be written using any expression involving rational number additions, subtractions, multiplications, divisions or \(k\) -th roots.
++All the methods for eigenvalues must be iterative!
+
Even if working with exact arithmetic, there could be no computing program that would produce the exact roots of an arbitrary polynomial of degree \(m\geqslant 5\) in a finite number of steps. All eigenvalue methods must be iterative.
+We want to convert the original matrix to a matrix that is easy to recognize the eigenvalues (upper/lower triangular or diagonal), and the transformation should not change the eigenvalues of matrices (this is the concept of similar).
+Commonly used numerical methods for eigenvalue problem have two phases:
+Phase 1: A direct method to produce an upper Hessenberg matrix \(\boldsymbol{H}\) (an upper Hessenberg matrix has zero entries below the first subdiagonal, and a lower Hessenberg matrix has zero entries above the first superdiagonal):
+The procedure can end within finite number of steps. This is done by similar transforms. In linear algebra, \(\boldsymbol{A}\) and \(\boldsymbol{B}\) are similar ( \(\boldsymbol{A}\sim \boldsymbol{B}\) ), if there exists \(\boldsymbol{X}\in \mathbb{R} ^{m\times m}\) invertible such that \(\boldsymbol{XAX}^{-1}=\boldsymbol{B}\). We can let:
+Because \(\det \left( \boldsymbol{X} \right) \ne 0, \det \left( \boldsymbol{X}^{-1} \right) \ne 0\), we know that the characteristic polynomials of \(\boldsymbol{A}\) and \(\boldsymbol{B}\) are the same, and \(\boldsymbol{A}\) and \(\boldsymbol{B}\) have the same set of eigenvalues.
+Phase 2: An iterative procedure that is also based on similarity transforms to produce a formally infinite sequence of upper Hessenberg matrices that converges to a triangular form.
+Why we divide the algorithms into two phases? The cost concern is the reason for the two-phase strategy. Phase 1 has \(O(m^3)\) flops, while in Phase 2, it normally converges to \(O\left( \varepsilon _{machine} \right)\) within \(O(m)\) steps, each step requiring at most \(O(m^2)\) steps to finish.
+Building blocks: QR Factorization
+ + + + + + +QR Algorithm is the recommended algorithm for eigenvalue and eigenvector computation.
+The matrix \(\boldsymbol{A}\) below has been computed after phase 1.
+Pure QR Algorithm:
+Initialize: \(\boldsymbol{A}^{\left( 0 \right)}=\boldsymbol{A}\);
+For \(k=1,2,\cdots\):
+End
+Task 1: Verify that QR algorithm performs similarity transformation, that is, \(\boldsymbol{A}^{\left( k \right)}\sim \boldsymbol{A}^{\left( k-1 \right)}\).
+If \(\boldsymbol{A}^{\left( k \right)}\) converges to \(\mathrm{diag}\left( \lambda _1,\cdots ,\lambda _n \right)\) as \(k\rightarrow +\infty\), then \(\bar{\boldsymbol{Q}}^{\left( k \right)}\rightarrow \boldsymbol{Q}\) which contains the eigenvectors of \(\boldsymbol{A}\).
+Task 2: Why do we need phase 1 before the algorithm performs?
+Let us assume \(\boldsymbol{A}\) is symmetric for simplicity, then \(\boldsymbol{H}=\boldsymbol{QAQ}^T\) is tridiagonal. Starting from \(\boldsymbol{H}\), the cost for each iteration is of \(O(m)\) (using Given Rotation, for example).
+If we don't perform phase 1, starting from \(\boldsymbol{A}\), the cost for each iteration is \(O(m^3)\) (using Householder QR factorization, for example).
+Question: What is the cost per iteration if \(\boldsymbol{H}\) is and upper Hessenberg matrix? (HW Practice)
+Task 3: QR algorithm is convergent.
++++
Theorem: Let the pure QR algorithm be applied to a real symmetric matrix \(\boldsymbol{A}\) whose eigenvalues satisfy \(\left| \lambda _1 \right|>\left| \lambda _2 \right|>\cdots >\left| \lambda _m \right|\), and the corresponding eigenmatrix \(\boldsymbol{Q}\) has all nonsingular leading principal submatrices. Then as \(k\rightarrow +\infty\), \(\boldsymbol{A}^{\left( k \right)}\) converges linearly with a constant rate given by \(\max_j \frac{\left| \lambda _{j+1} \right|}{\left| \lambda _j \right|}\) to diagonal matrix \(\mathrm{diag}\left( \lambda _1,\lambda _2,\cdots ,\lambda _n \right)\), and \(\bar{\boldsymbol{Q}}^{\left( k \right)}\) (with signs of its columns adjusted as necessary) converges to \(\boldsymbol{Q}\) at the same rate.
In the theorem above, \(\boldsymbol{A}=\boldsymbol{Q}^T\mathbf{\Lambda }\boldsymbol{Q}\) where \(\boldsymbol{Q}\in \mathbb{C} ^{m\times m}\) is unitary.
+The sketch of the proof (in three steps):
+Step 1: QR algorithm is essentially equivalent to the so-called simultaneous iteration.
+Simultaneous Iteration:
+Pick \(\hat{\boldsymbol{Q}}^{\left( 0 \right)}\in \mathbb{R} ^{m\times n}\) with orthonormal columns;
+For \(k=1,2,\cdots\):
+End
+Claim: If we pick \(\hat{\boldsymbol{Q}}^{\left( 0 \right)}=\mathbf{I}\), then the simultaneous iteration becomes the QR algorithm. (How to verify?)
+Step 2: Simultaneous iteration is a block power iteration.
+Step 3: The block power iteration is convergent with the rate given by
+if \(\left| \lambda _1 \right|>\left| \lambda _2 \right|>\cdots >\left| \lambda _j \right|>\left| \lambda _{j+1} \right|>\cdots >\left| \lambda _m \right|\geqslant 0\). Then we get \(0<c<1\).
+A special case for block power iteration to show the convergence (this is the essential idea):
+Consider \(\left[ \begin{matrix} + {\boldsymbol{v}_1}^{\left( 0 \right)}& {\boldsymbol{v}_{\left( 2 \right)}}^{\left( 0 \right)}\\ +\end{matrix} \right]\) where \({\boldsymbol{v}_1}^{\left( 0 \right)}\) and \({\boldsymbol{v}_{\left( 2 \right)}}^{\left( 0 \right)}\) are unit vectors and orthogonal to each other. Then:
+Also:
+\(\boldsymbol{Q}\) spans the space generated by the eigenvectors \(\boldsymbol{q}_1, \boldsymbol{q}_2\). We can further prove that \({\boldsymbol{v}_1}^{\left( k \right)}\rightarrow \boldsymbol{q}_1, {\boldsymbol{v}_2}^{\left( k \right)}\rightarrow \boldsymbol{q}_2\) as \(k\rightarrow +\infty\)
+We assume \(\boldsymbol{A}\) is SPD below.
+Phase 1 is the same as before: \(\left( \boldsymbol{Q}^{\left( 0 \right)} \right) ^T\boldsymbol{A}^{\left( 0 \right)}\boldsymbol{Q}^{\left( 0 \right)}=\boldsymbol{A}\) where \(\boldsymbol{A}^{\left( 0 \right)}\) is a tridiagonal matrix;
+For \(k=1,2,\cdots\):
+End
+Two key modifications:
+How can we pick the shift?
+Method 1: we can pick \(\mu ^{\left( k \right)}={\boldsymbol{A}_{m,m}}^{\left( k \right)}\) as the last element of \(\boldsymbol{A}^k\).
+As \(k\rightarrow +\infty\), \(\boldsymbol{q}^{\left( k \right)}\rightarrow \boldsymbol{q}_m, \mu ^{\left( k \right)}\rightarrow \lambda _m\)
+This strategy may fail in some situations. For example, if \(\boldsymbol{A}=\left[ \begin{matrix} + 0& 1\\ + 1& 0\\ +\end{matrix} \right]\), we get \(\lambda _{1,2}=\pm 1\). \(\boldsymbol{AI}=\boldsymbol{A},\boldsymbol{RQ}=\boldsymbol{A}\) remains unchanged. The shift is zero.
+Method 2 (Wilkerson's Shift):
+Let:
+as a lower-right 2 by 2 matrix. Pick the eigenvalue of \(\boldsymbol{B}\), closer to \(a_m\), as the shift \(\mu ^{\left( k \right)}\), In the case of a tie, just pick and arbitrary one:
+where \(\delta =\frac{a_{m-1}-a_m}{2}\). If \(\delta =0\), set \(\delta =\pm 1\) arbitrarily.
+Claims:
+QR Factorization definition:
+\(\boldsymbol{Q}\) is a unitary matrix and \(\boldsymbol{R}\) is an upper triangular matrix.
+Gram-Schmidt is an algorithm to produce \(\boldsymbol{Q}\) and \(\boldsymbol{R}\).
+Given \(\left[ \begin{matrix} + \boldsymbol{a}_1& \boldsymbol{a}_2& \cdots& \boldsymbol{a}_m\\ +\end{matrix} \right]\). Let \(\mathrm{span}\left\{ \boldsymbol{a}_1, \boldsymbol{a}_2, \cdots , \boldsymbol{a}_m \right\} \triangleq \left< \boldsymbol{a}_1, \boldsymbol{a}_2, \cdots , \boldsymbol{a}_m \right>\). We need to find \(\boldsymbol{q}_1, \boldsymbol{q}_2, \cdots , \boldsymbol{q}_m\) such that \(\left< \boldsymbol{q}_1, \boldsymbol{q}_2, \cdots , \boldsymbol{q}_m \right> =\left< \boldsymbol{a}_1, \boldsymbol{a}_2, \cdots , \boldsymbol{a}_m \right>\) and the inner products \(\left( \boldsymbol{q}_i, \boldsymbol{q}_j \right) =0\) if \(i\ne j\), \(\left( \boldsymbol{q}_i, \boldsymbol{q}_i \right) =1\).
+Example steps (only 2-norm is used below):
+Given \(\boldsymbol{A}=\left[ \begin{matrix} + \boldsymbol{a}_1& \cdots& \boldsymbol{a}_j& \cdots& \boldsymbol{a}_m\\ +\end{matrix} \right]\),
+For \(j=1:m\):
+End
+The outcome: \(\boldsymbol{Q}=\left[ \begin{matrix} + \boldsymbol{q}_1& \boldsymbol{q}_2& \cdots& \boldsymbol{q}_m\\ +\end{matrix} \right]\) unitary, and \(\boldsymbol{R}=\left[ r_{ij} \right]\) upper triangular.
+Claim: The classical Gram-Schmidt algorithm above is not stable.
Given \(\boldsymbol{A}=\left[ \begin{matrix} + \boldsymbol{a}_1& \cdots& \boldsymbol{a}_j& \cdots& \boldsymbol{a}_m\\ +\end{matrix} \right]\),
+For \(i=1:m\):
+End
+For \(i=1:m\):
+End
+Claim: The modified Gram-Schmidt algorithm is stable.
Questions:
+The cost of the modified Gram-Schmidt is \(O(2mn^2)\) for \(\boldsymbol{A}\in \mathbb{R} ^{m\times n}\).
+Theorem: Every matrix \(\boldsymbol{A}\in \mathbb{R} ^{m\times n}\,\,\left( m\geqslant n \right)\) has a QR factorization.
If column vectors of \(\boldsymbol{A}\) are linearly independent(full rank), then
+This is reduced QR factorization(default). The other one:
+is called full QR factorization(rarely used).
+Theorem: If \(\boldsymbol{A}\in \mathbb{R} ^{m\times n}\,\,\left( m\geqslant n \right)\) is full rank, the reduced QR factorization is unique if we require \(r_{ii}>0, i=1,\cdots ,n\).
We can use QR factorization to solve linear systems (least square solution):
+Matrix explanation of Gram-Schmidt process (Using upper triangular matrices to generate orthogonal matrix):
+Two types of orthogonal linear transformation:
+Actually, there are other ways to computer QR factorization using unitary matrices to generate upper triangular matrix (Using unitary transformation):
+where \(\boldsymbol{Q}_i\) are unitary. There are two ways for this process: Householder transformation and Givens rotation.
+The overall procedure (recursive):
+Example step 1:
+where \(\boldsymbol{Q}_1\) is unitary.
+
We can get:
+Question: why not choose \(\mathbf{H}_2\) and get \(\boldsymbol{v}=\left\| \boldsymbol{x} \right\| \boldsymbol{e}_1-\boldsymbol{x}\)?
+We have two choices: \(\boldsymbol{Q}_1\boldsymbol{x}=-\left\| \boldsymbol{x} \right\| \boldsymbol{e}_1, \boldsymbol{Q}_1\boldsymbol{x}=\left\| \boldsymbol{x} \right\| \boldsymbol{e}_1\).
+In practice, we select:
+\(\boldsymbol{x}_1\) is the first entry of \(\boldsymbol{x}\). We choose it for the stability reason. This results in a larger \(\left\| \boldsymbol{v} \right\|\).
+Example step 2:
+We can get:
+Here is the algorithm for Householder QR factorization for a matrix \(\boldsymbol{A}\in \mathbb{R} ^{m\times n}\) :
+For \(k=1:m\):
+End
+Remarks:
+The core element of Givens Rotation:
+Transformation pattern:
+Introduction to Eigenvalue Problems
+ + + + + + + + + +