Gram-Schmidt clásico para la matriz $A$

Question

Deje que $$A= \begin {bmatrix}1 & 1 & 1 \\ \epsilon & 0 & 0 \\ 0\ & \epsilon & 0 \\ 0 & 0 & \epsilon \end {bmatrix}.$$ En esto página esta matriz $A$ se utiliza para mostrar la inestabilidad del clásico algoritmo de Gram-Schmidt, utilizando el criterio de que $1+ \epsilon =1$ . Además, se puede demostrar que los vectores de salida de la GS clásica para $A$ no son ortogonales entre sí.

Parece que muchos sitios web brevemente parecen hablar sólo de los defectos del algoritmo al ejecutarlo en un ordenador. ¿Hay algún razonamiento más "general" de por qué el clásico algoritmo GS no siempre produce vectores ortonormales, incluso "en papel"?

¿Es porque la GS clásica (en este caso) no explica bien la aproximación $\epsilon +1=1$ ? ¿Alguien podría explicar esto un poco más a fondo?

Gracias

Answer 1

Tanto la Gram Schmidt clásica como la modificada son inestables. Si lees el texto de Trefethen describe la diferencia entre Householder y los dos primeros como lo siguiente.

Se trata de la Gram-Schmidt clásica y modificada, descrita Ortogonalización triangular $A \underbrace{R_{1} , R_{2} \cdots R_{n}}_{\hat{R}^{-1}} = \hat{Q} \tag{1}$

Abajo vemos a Householder, Triangularización ortogonal

$\underbrace{Q_{1} , Q_{2} \cdots Q_{n}}_{\hat{Q}^{*}}A = R \tag{2}$

¿Por qué son diferentes?

El número de condición de una matriz triangular puede ser cualquier cosa así que si tienes una serie de ellas entonces puede crecer muy grande sin embargo las matrices ortogonales tienen número de condición $1$ .

Al cambiar el $\epsilon$ se cambia el número de condición. Si se da cuenta de que $\epsilon$ está relacionado con los valores singulares. El primero es casi $1$ .

import numpy as np
import math

eps = math.exp(1e-3)-1
A = np.matrix([[1 ,1,1],[eps, 0 ,0 ], [0 ,eps, 0], [0 , 0 ,eps ]])
u, s, vt = np.linalg.svd(A)

s

Out[12]: array([1.73205110e+00, 1.00050017e-03, 1.00050017e-03])

eps
Out[13]: 0.0010005001667083846

Debido a la ortogonalización parece ser $\sqrt{3}$

Tenga en cuenta que

$$\kappa(A) = \frac{\sigma_{max}(A)}{\sigma_{min}(A)} = \frac{\sqrt{3}}{\epsilon} \tag{3}$$

Entonces notarías que como $\epsilon \to 0$ $\kappa \to \infty$

Clásica Gram Schmidt

El proceso de Gram Schmidt es el siguiente para los clásicos

$$v_{j} = a_{j} - (q_{1}^{*}a_{j})q_{1} -(q_{2}^{*}a_{j})q_{2} - \cdots - (q_{j-1}^{*}a_{j})q_{j-1} \tag{3}$$

podemos escribirlo así

$$q_{1} = \frac{a_{1}}{r_{11}} \tag{4}$$

$$q_{2} = \frac{a_{2} - r_{12}q_{1}}{r_{22}} \tag{5}$$

$$q_{3} = \frac{ a_{3} - r_{13} q_{1}- r_{23}q_{2} }{r_{33}} \tag{6}$$ $$q_{n} = \frac{a_{n} - \sum_{i=1}^{n-1} r_{in} q_{i} }{r_{nn} } \tag{7}$$

Ahora aquí está la Gram Schmidt modificada. Para empezar introducimos las proyecciones ortogonales

Gram Schmidt modificado

$$q_{1} = \frac{P_{1}a_{1}}{\| P_{1}a_{1}\|}, q_{2} = \frac{P_{2}a_{2}}{\| P_{2}a_{2}\|}, \cdots , q_{n} = \frac{P_{n}a_{n}}{\| P_{n}a_{n}\|} \tag{8}$$

Más concretamente $P_{j}$ es el proyector ortogonal. $P_{j}$ es el $m \times m$ matriz de rango $m -(j-1)$ que proyecta $\mathbb{C}^{m}$ en el espacio para $\langle q_{1}, \cdots , q_{j-1} \rangle$

El proyector $P_{j}$ puede representarse explícitamente. Aquí representamos $\hat{Q}_{j-1}$ como el $m \times (j-1)$ matriz que contiene las columnas de la proyección orhtogonal. Es decir

$$P_{j} = I - \hat{Q}_{j-1}\hat{Q}_{j-1}^{*} \tag{9}$$

entonces obtenemos

$$v_{j} = P_{j}a_{j} \tag{10}$$

¿Cómo es esto más estable?

Una nota más

Su matriz es famosa. Es llamado el Matriz de Lauchli

Tanto en CGS como en MGS

donde $1 + \epsilon^{2} =1$

$$v_{1} \to (1 , \epsilon, 0, 0) \tag{11}$$

$$r_{11} = \sqrt{1 + \epsilon^{2} } \approx 1 \tag{12}$$

$$q_{1} = \frac{v_{1}}{r_{11}} = (1 , \epsilon, 0, 0)\tag{13}$$ $$v_{2} = (1,0,\epsilon,0) \tag{14}$$ $$r_{12} = q_{1}^{T}a_{2} = q_{1}^{T}v_{2} = 1 \tag{15}$$ $$v_{2} = v_{2} - r_{12}q_{1} = (0,-\epsilon, \epsilon,0) \tag{16}$$ $$r_{22} = \sqrt{2}\epsilon \tag{17}$$ $$q_{2} = (0,\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}},0) \tag{18}$$

$$v_{3} = (1,0,0,\epsilon) \tag{19}$$ $$r_{13} = q_{1}^{t}v_{3} = 1 \tag{20}$$ $$v_{3} = v_{3} - r_{13}q_{1} = (0,-\epsilon,0,\epsilon) \tag{21}$$

Para CGS

$$r_{23} = q_{2}^{T}a_{3} =0 \tag{22}$$ $$v_{3} = v_{3} - r_{23}q_{2} = (0,-\epsilon,0,\epsilon) \tag{23}$$

$$r_{33} = \sqrt{2} \epsilon \tag{24}$$ $$q_{3} = \frac{v_{3}}{r_{33}} = (0,\frac{-1}{\sqrt{2}} ,0\frac{1}{\sqrt{2}} ) \tag{25}$$

Para MGS

$$r_{23} = q_{2}^{T}v_{3} =\frac{\epsilon}{\sqrt{2}} \tag{26}$$ $$v_{3} = v_{3} - r_{23}q_{2} = (0,\frac{-\epsilon}{2},\frac{-\epsilon}{2}, \epsilon ) \tag{27}$$

$$r_{33} = \frac{\sqrt{6}}{\epsilon 2} \tag{28}$$ $$q_{3} = \frac{v_{3}}{r_{33}} = (0,\frac{-1}{\sqrt{6}} ,\frac{1}{\sqrt{6}},\frac{2}{\sqrt{6}} ) \tag{29}$$