Cholesky Factorization and SPD Matrices

Least Squares Problems. Let \(A \in \mathbb {R}^{m \times n}\) with \(m > n\) have rank \(n\) and \(b \in \mathbb {R}^m\). The Least Squares problem is to find \(x \in \mathbb {R}^n\) which makes \(\|Ax - b\|_2\) as small as possible.

1. Show that \(x\) solves this Least Squares problem if and only if \(A^T Ax = A^T b\).
2. What is meant by the QR factorization of matrix \(A\)?
3. Show how to use the QR factorization of \(A\) to solve the Least Squares problem.

(QR and Cholesky Factorizations). Let \(A\) be a nonsingular square matrix and let \(A = QR\) and \(A^*A = U^*U\) be the QR factorization of \(A\) and the Cholesky factorization of \(A^*A\), respectively. Assume that the usual normalizations \(r_{jj}, u_{jj} > 0\) are in effect. Is it true or false that \(R = U\)? Justify your answer with mathematical details.

(Properties of the SPD Matrices) Let \(A \in \mathbb {R}^{n \times n}\) be a symmetric positive definite matrix (spd). Show that there exists a unique upper triangular matrix \(H \in \mathbb {R}^{n \times n}\) with positive diagonal entries such that, \(A = H^T H\).

(QR and Cholesky Factorizations). Let \(A\) be a nonsingular square matrix and let \(A = QR\) and \(A^*A = U^*U\) be the QR factorization of \(A\) and the Cholesky factorization of \(A^*A\), respectively. Assume that the usual normalizations \(r_{jj}, u_{jj} > 0\) are in effect. Is it true or false that \(R = U\)? Justify your answer with mathematical details.

Least Squares Problems. Let \(A \in \mathbb {R}^{m \times n}\) with \(m > n\) have rank \(n\) and \(b \in \mathbb {R}^m\). The Least Squares problem is to find \(x \in \mathbb {R}^n\) which makes \(\|Ax - b\|_2\) as small as possible.

1. Show that \(x\) solves this Least Squares problem if and only if \(A^T Ax = A^T b\).
2. What is meant by the QR factorization of matrix \(A\)?
3. Show how to use the QR factorization of \(A\) to solve the Least Squares problem.

Matrix Factorizations

Let \(A\) be a nonsingular square real matrix and let \(A=QR\) and \(A^TA=H^TH\) be QR and Cholesky factorizations, respectively, with the usual normalizations \(r_{jj},h_{jj}>0\). Is it true or false that \(R=H\)? Justify your answer.

CG Method

a) Formulate the Conjugate Gradient (CG) algorithm for the solution of the linear system \(Ax=b\), where \(A\in \mathbb {R}^{n\times n}\) is symmetric positive definite matrix.

b) Show that for any initial guess \(x_0\in \mathbb {R}^n\), Conjugate Gradient algorithm converges to the solution \(x\) of the linear system \(Ax=b\) in at most \(n\) steps.

Positive Definite Matrices The Cholesky decomposition of a (symmetric and) positive definite matrix is defined to be a factorization of the form \(A=LL^T\) where \(L\) is lower triangular. Show that the Cholesky decomposition exists. Discuss whether or not it is unique.

Suppose \(A\) is an \(m\times n\) real matrix with \(m>n\), and \(b\in \mathbb {R}^m\). For any vector \(v\in \mathbb {R}^m\), let \[ \|v\|_2=\sqrt {\sum _{i=1}^m v_i^2} \] as usual. Describe how to solve the Least Squares Problem \[ \|Ax-b\|_2=\min \] using the \(QR\) factorization of \(A\). Briefly state how you would compute the \(QR\) factorization.

\(QR\) and Cholesky Factorizations

Let \(A\) be a nonsingular square matrix and let \(A=QR\) and \(A^*A=U^*U\) be the \(QR\) factorization of \(A\) and the Cholesky factorization of \(A^*A\), respectively. Assume that the usual normalizations \(r_{jj},u_{jj}>0\) are in effect. Is it true or false that \(R=U\)? Justify your answer.

QR factorization Show how to compute the QR factorization of a real \(m\times n\) matrix \(A\) (with \(m\geq n\)) using Householder reflections, and describe how to use the QR factorization to solve the Least Squares problem \[ \|Ax-b\|_2=\min . \tag {1} \]

Least Squares

Suppose \(A\in \mathbb {R}^{m\times n}\) has linearly independent columns, and \(C\in \mathbb {R}^{\ell \times n}\) has linearly independent rows. Let \(d\in \mathbb {R}^\ell \) and \(b\in \mathbb {R}^m\). Show how two full-rank QR factorizations can be used to solve the linearly constrained least squares problem: \[ \text {Find } x^*\in \{x\in \mathbb {R}^n:Cx=d\} \text { which minimizes } \|b-Ax^*\|_2. \]

Sensitivity Consider a \(6\times 6\) symmetric positive definite matrix \(A\) with singular values \(\sigma _1=1000\), \(\sigma _2=500\), \(\sigma _3=300\), \(\sigma _4=20\), \(\sigma _5=1\), \(\sigma _6=0.01\).

a) Suppose you use a Cholesky factorization package on a computer with a machine epsilon \(10^{-14}\) to solve the system \(Ax=b\) for some nonzero vector \(b\). How many digits of accuracy do you expect in the computed solution? Justify your answer in terms of condition and stability. You may assume that the entries of \(A\) and \(b\) are exactly represented in the computer’s floating-point system.

b) Suppose that instead you use an iterative method to find an approximate solution to \(Ax=b\) and you stop iterating and accept iterate \(x^{(k)}\) when the residual \(r^{(k)}=Ax^{(k)}-b\) has 2-norm less than \(10^{-9}\). Give an estimate of the maximum size of the relative error in the final iterate? Justify your answer.

Iterative Methods for Linear Systems Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \), \[ \alpha _n=\frac {r_{n-1}^T r_{n-1}}{p_{n-1}^T A p_{n-1}}, \] \[ x_n=x_{n-1}+\alpha _n p_{n-1}, \] \[ r_n=r_{n-1}-\alpha _n A p_{n-1}, \] \[ \beta _n=\frac {r_n^T r_n}{r_{n-1}^T r_{n-1}}, \] \[ p_n=r_n+\beta _n p_{n-1}. \] End.

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &= \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} = \operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &= \operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} = \operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^T Ax)^{1/2}\).

a) Prove the following: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N\to 0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Householder Reflections A Householder Reflection \(H\) is an \(n\times n\) orthogonal matrix of the form \[ H=I-2uu^T \tag {1} \] where \[ \lVert u\rVert _2=1 \tag {2} \] and \(I\) is the \(n\times n\) identity matrix. Let \(a\in \mathbb {R}^n\) be a given (non-zero) vector. Show how to pick \(u\) in (1) such that \[ Ha=\lVert a\rVert _2 e_1 \tag {3} \] where \(e_1\) is the first unit vector. Describe how to use this technique to compute the \(QR\) factorization of an \(n\times n\) matrix \(A\).

The \(QR\) factorization Describe the \(QR\) factorization of an \(m\times n\) matrix \(A\), how to compute it, and what it can be used for.

Krylov Subspace Methods Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \) \[ \begin {aligned} \alpha _n &= r_{n-1}^T r_{n-1}/p_{n-1}^T A p_{n-1},\\ x_n &= x_{n-1}+\alpha _n p_{n-1},\\ r_n &= r_{n-1}-\alpha _n A p_{n-1},\\ \beta _n &= r_n^T r_n/r_{n-1}^T r_{n-1},\\ p_n &= r_{n-1}+\beta _n p_{n-1}. \end {aligned} \] End

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &\equiv \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} =\operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &=\operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} =\operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^TAx)^{1/2}\).

a) Prove the following:

Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N=0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Krylov Subspace Methods Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \) \[ \begin {aligned} \alpha _n &= r_{n-1}^T r_{n-1}/p_{n-1}^T A p_{n-1},\\ x_n &= x_{n-1}+\alpha _n p_{n-1},\\ r_n &= r_{n-1}-\alpha _n A p_{n-1},\\ \beta _n &= r_n^T r_n/r_{n-1}^T r_{n-1},\\ p_n &= r_{n-1}+\beta _n p_{n-1}. \end {aligned} \] End

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &\equiv \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} =\operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &=\operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} =\operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^TAx)^{1/2}\).

a) Prove the following:

Theorem: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N=0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

The QR Factorization. Let \(A\) be an \(m\times n\) matrix, where \(m\ge n\). What is the QR factorization of \(A\)? How do you compute it? What do you use it for?

Gaussian Elimination is LU Factorization. What is the LU factorization of a square matrix? How do you solve a linear system \(Ax=b\) by Gaussian Elimination? Show that for a general \(4\times 4\) system the matrix part of Gaussian Elimination is equivalent to computing the LU factorization of \(A\).

a) What is the best least squares approximation to \(f(x)=\cos 3x\) by functions of the form \[ a_0+a_1\cos x+a_2\cos ^2 x+b_1\sin x+b_2\sin ^2 x \] with the inner product given by \[ (v,w)_C=\frac {1}{2\pi }\int _0^{2\pi }v(x)\overline {w(x)}\,dx \]

b) What if the inner product is \[ (v,w)_D=\frac {1}{J}\sum _{j=0}^{J-1}v(x_j)\overline {w(x_j)} \] with \(x_j=jh\), \(h=2\pi /J\).

c) Briefly discuss the normal equations (do they form a well- or ill-conditioned system?), the Gram-Schmidt process, QR factorization in the context of parts a) and b). Observe the connection with problem #2.