6 - More on Solving Linear Systems

ตอนนี้รู้แล้วว่าถ้าเจอ Linear System ที่เป็น Lower Triangular ให้ใช้ Forward Substitution
ถ้าเป็น Upper Triangular ควรใช้ Backward Substitution
แล้วถ้า General Case ก็สามารถใช้ GEPP + 4 steps procedure

มันมี special case ของ matrix ที่เราสามารถทำได้ดีกว่าใช้ GEPP อีก

Positive Definite and Positive Semidefinite Matrices

Definitions

• A symmetric matrix $A \in \mathbb{R}^{n \times n}$ is positive definite if:
  $\boxed{x^T Ax > 0}$
  for all nonzero $x \in \mathbb{R}^n$

• A symmetric matrix $A \in \mathbb{R}^{n \times n}$ is positive semidefinite if:
  $\boxed{x^T Ax \geq 0}$
  for all $x \in \mathbb{R}^n$

• Note: Both are usually defined for only symmetric matrices (as unsymmetric positive definite/semidefinite matrices are not that useful in practice)

Positive Definite / Positive Semidefinite นี่เป็น special case ของ matrix ที่ “สวย” และ “แก้สมการได้ง่ายกว่าใช้ Gaussian Elimination (GEPP)” เพราะเรามี Cholesky decomposition โดยเฉพาะ

Cholesky Factorization

Overview

• Solving linear system: $Ax = b$
• If the matrix $A$ is symmetric positive definite, its LU factorization can be arranged so that $U = L^T$

Theorem

Any symmetric positive definite matrix $A \in \mathbb{R}^{n \times n}$ can be factored:
$\boxed{A = LL^T}$
where $L$ is lower triangular (not necessarily unit lower triangular).

• This factorization is known as Cholesky factorization
• Cholesky factorization can be computed from GE, but there is a more efficient way

Computing Cholesky Factorization on a 2×2 Matrix

Given: $A = LL^T$
$\begin{bmatrix} a_{11} & a_{21} \\ a_{21} & a_{22} \end{bmatrix} = \begin{bmatrix} l_{11} & 0 \\ l_{21} & l_{22} \end{bmatrix} \begin{bmatrix} l_{11} & l_{21} \\ 0 & l_{22} \end{bmatrix} = \begin{bmatrix} l_{11}^2 & l_{11}l_{21} \\ l_{11}l_{21} & l_{21}^2 + l_{22}^2 \end{bmatrix}$

Solution steps:

• $l_{11}^2 = a_{11} \Rightarrow \boxed{l_{11} = \pm\sqrt{a_{11}}}$

• $l_{11}l_{21} = a_{21} \Rightarrow \boxed{l_{21} = a_{21}/l_{11}}$

• $l_{21}^2 + l_{22}^2 = a_{22} \Rightarrow \boxed{l_{22} = \pm\sqrt{a_{22} - l_{21}^2}}$

• By convention, we always choose the positive values for the diagonal entries of $L$

• The same idea can be extended to a more general $n$-by-$n$ matrix

Exercise 1 Solution

Compute the Cholesky factorization of:
$\begin{bmatrix} 9 & -15 \\ -15 & 26 \end{bmatrix}$
Solution:

• $l_{11} = \sqrt{a_{11}} = \sqrt{9} = 3$
• $l_{21} = a_{21}/l_{11} = -15/3 = -5$
• $l_{22} = \sqrt{a_{22} - l_{21}^2} = \sqrt{26 - (-5)^2} = \sqrt{26 - 25} = 1$

Therefore:
$L = \begin{bmatrix} 3 & 0 \\ -5 & 1 \end{bmatrix}$

And: $\begin{bmatrix} 9 & -15 \\ -15 & 26 \end{bmatrix} = LL^T$

Cholesky Factorization Algorithm

for k = 1 to n
    a_kk = sqrt(a_kk)
    for i = k + 1 to n
        a_ik = a_ik / a_kk
    end
    for j = k + 1 to n
        for i = k + 1 to n
            a_ij = a_ij - a_ik * a_jk
        end
    end
end

Good Points of Cholesky Factorization

• The algorithm is well-defined. The square roots required are all of positive numbers (unless $A$ is not positive definite)
• No pivoting required for numerical stability (unlike GE)
• Is only $n^3/3 + O(n^2)$ operations. Cheaper than GE/GEPP by a factor of 2 (GE/GEPP are $2n^3/3 + O(n^2)$ operations)

Theorem: In exact arithmetic, the Cholesky factorization algorithm runs to completion if and only if the matrix $A$ is symmetric positive definite.

Solving Linear System with Cholesky

To solve $Ax = b$ when $A$ is symmetric positive definite:

1. Factor $A = LL^T$
2. Solve $Lw = b$ for $w$ (forward substitution)
3. Solve $L^T x = w$ for $x$ (back substitution)

Verify:

• Let $x$ be the output of the above procedure

• See that: $Ax = (LL^T)x = L(L^T x) = Lw = b$

• This algorithm is backward stable for solving symmetric positive definite linear system

Cholesky Factorization in MATLAB

>> R = chol(A);        % Return R = L^T (so A = R^T*R)
>> L = chol(A,'lower'); % Return L (A = L*L^T)

Improving Accuracy of Linear Systems

Example Problem

Consider the linear system:
$\begin{bmatrix} 1 & 0 \\ 0 & \epsilon \end{bmatrix} x = \begin{bmatrix} 1 \\ \epsilon \end{bmatrix}$

• The matrix has condition number (in p-norm) is $1/\epsilon$, which is very ill-conditioned if $\epsilon$ is very small
• If we multiply the second row by $1/\epsilon$, then the system becomes:
  $\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} x = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$
  which is perfectly well-conditioned

Iterative Refinement

Iterative Refinement is another means of potentially improving the accuracy of a computed solution.

Process

Suppose we have computed an approximate solution $x_0$ to the system $Ax = b$.

1. We compute the residual: $r_0 = b - Ax_0$
2. If $r_0$ is large (say, because the system is ill-conditioned or we used an iterative method), we can solve: $As_0 = r_0$
3. Take: $x_1 = x_0 + s_0$ as a new and "better" approximate solution

Verification

$x_1$ is a solution because:
$Ax_1 = A(x_0 + s_0) = Ax_0 + As_0 = (b - r_0) + r_0 = b$

Process Continuation and Limitations

• This process can be repeated successively to reduce the residual as needed
• However, iterative refinement requires double the storage as we need to store both $A$ (to compute the residual) and its factors (to solve the subsequent systems)
• Moreover, if the residual is already small, the residual must usually be computed with higher precision than that used in computing $x_0$

Solving a Tridiagonal Linear System

Problem Setup

Consider solving $Ax = b$ where $A$ is a tridiagonal matrix: