0

I was trying to understand why LU factorization is more numerically stable than the method of inverting A directly. This answer talks about it. It says the below:

$$ \begin{array}{l}{\text { If } x \text { is computed by LU factorization, the residual can be bounded by }} \\ {\qquad\|b-A x\| \leq c u\|L\| U\|\| x \|}\end{array} $$ Can someone show a step-by-step way of proving this result?

I referred the textbook "Accuracy and Stability of Numerical Algorithms" by Nicholas J Higham , but I was not able to understand.

Community
  • 1
wanderer
  • 155
  • 3
  • 13

1 Answers1

0

The proof can be based on the following facts (see Theorems 9.3 and 9.4 in the first edition of the referenced book):

Fact 1 Let $T\in\mathbb{R}^{n\times n}$ be a nonsingular triangular matrix (upper or lower). Then the computed solution $\hat{x}$ of $Tx=b$ satisfies $$ (T+\Delta T)\hat{x}=b, \quad |\Delta T|\leq nu|T|+O(u^2). $$

Fact 2 Computed LU factors $\hat{L}$ and $\hat{U}$ of $A\in\mathbb{R}^{n\times n}$ satisfy $$ \hat{L}\hat{U}=A+\Delta A, \quad |\Delta A|\leq nu|\hat{L}||\hat{U}|+O(u^2). $$ Here the inequalities are understood component-wise.

We start by computing the LU factorization to obtain $\hat{L}$ and $\hat{U}$. Then we obtain $\hat{y}$ by solving $\hat{L}y=b$ such that $$ (\hat{L}+\Delta L)\hat{y}=b,\quad|\Delta L|\leq nu|\hat{L}|+O(u^2), $$ and $\hat{x}$ by solving $\hat{U}x=y$ such that $$ (\hat{U}+\Delta U)\hat{x}=\hat{y},\quad |\Delta U|\leq nu|\hat{U}|+O(u^2). $$ At the end, we have $$ (\hat{L}+\Delta L)(\hat{U}+\Delta U)\hat{x}=b. $$ Now for the residual, we get $$ \begin{split} b-A\hat{x} &=(\hat{L}+\Delta L)(\hat{U}+\Delta U)\hat{x}-(\hat{L}\hat{U}-\Delta A)\hat{x} \\&=(\Delta L\hat{U}+\hat{L}\Delta U+\Delta L\Delta U+\Delta A)\hat{x}=: E\hat{x}, \end{split} $$ where $$ |E|\leq|\Delta L||\hat{U}|+|\hat{L}||\Delta U|+|\Delta A|+|\Delta L||\Delta U| \leq 3nu|\hat{L}||\hat{U}|+O(u^2). $$ We absorbed $|\Delta L||\Delta U|$ in the higher order term $O(u^2)$.

Algebraic Pavel
  • 23,818