Àlgebra lineal numèrica: condicionament, LU, Cholesky, QR i SVD

A classe d'àlgebra lineal, resoldre un sistema és aplicar Gauss i ja està. A l'ordinador no: hi ha arrodoniment, hi ha matrius que amplifiquen els errors i hi ha una diferència enorme entre «el residu és petit» i «la solució és bona». Aquest tema és el que necessites per no confiar cegament en el que et retorna un programa.

Condicionament: una propietat del problema

Per a un sistema amb matriu invertible, el número de condició en una norma subordinada és el producte de la norma de la matriu per la de la seva inversa, i sempre val com a mínim $1$. La seva utilitat és aquesta cota:

$$\frac{\lVert \delta x\rVert}{\lVert x\rVert} \leq \kappa(A)\,\frac{\lVert \delta b\rVert}{\lVert b\rVert}.$$

És a dir, un error relatiu a les dades es pot amplificar fins a $\kappa(A)$ vegades a la solució.

Error cap endavant, cap enrere i residu

L'error cap endavant compara la solució calculada amb l'exacta. L'error cap enrere pregunta una altra cosa: quina pertorbació de les dades faria que la solució calculada fos exacta. Un algorisme és estable cap enrere si aquesta pertorbació és petita, i això és el màxim que es pot demanar.

El residu és el que queda en substituir. Com que la diferència amb la solució exacta és la inversa aplicada al residu,

$$\frac{\lVert x - \widehat{x}\rVert}{\lVert x\rVert} \leq \kappa(A)\,\frac{\lVert r\rVert}{\lVert b\rVert}.$$

Un residu petit no garanteix una solució bona si el condicionament és gran.

LU amb pivotatge

L'eliminació gaussiana amb pivotatge parcial produeix una factorització $PA = LU$. Per resoldre, es fa una substitució endavant i una enrere. Factoritzar costa de l'ordre de $n^{3}$ operacions, però cada nou terme independent només costa $n^{2}$: si has de resoldre molts sistemes amb la mateixa matriu, factoritza un cop i reutilitza.

El pivotatge evita dividir per pivots nuls o massa petits en comparació amb els altres candidats de la columna. A la pràctica és molt estable, encara que la seva garantia teòrica depèn del creixement dels elements i no és absoluta.

Cholesky i QR

Si la matriu és simètrica i definida positiva, existeix una única factorització com a producte d'una triangular inferior amb diagonal positiva pel seu transposat. Cholesky aprofita l'estructura i fa aproximadament la meitat del treball de LU, a més de no necessitar pivotatge. Una matriu simètrica però indefinida no admet aquesta factorització amb diagonal real positiva, i intentar-la falla precisament en arribar a una arrel de número negatiu; això la converteix en un test de definició positiva.

Per a mínims quadrats amb més equacions que incògnites i rang columna complet, la factorització QR reduïda escriu la matriu com el producte d'una amb columnes ortonormals per una triangular invertible, i la solució surt de resoldre un sistema triangular amb el transposat aplicat al terme independent. La part del terme independent perpendicular a la imatge és el residu irreductible: no hi ha cap solució que el faci desaparèixer.

SVD, rang i regularització

La descomposició en valors singulars escriu qualsevol matriu com un producte de dues ortogonals amb una diagonal de valors singulars no negatius ordenats. És l'eina més informativa que hi ha: dona el rang, les bases dels quatre subespais fonamentals, el condicionament (quocient entre el valor singular més gran i el més petit) i la millor aproximació de rang baix.

La pseudoinversa inverteix només els valors singulars no nuls, i aplicada al terme independent dona la solució de mínims quadrats de norma mínima, també quan manca rang.

El problema és que els valors singulars molt petits amplifiquen el soroll: dividir per ells multiplica l'error. Per això el rang numèric necessita una tolerància que tingui en compte l'escala, la precisió de la màquina i el soroll de les dades; comparar exactament amb zero no té cap sentit en coma flotant.

La regularització de Tikhonov afegeix una penalització a la norma de la solució. En coordenades singulars, l'efecte és substituir cada factor $1/\sigma_i$ per $\sigma_i/(\sigma_i^{2}+\lambda)$, que està acotat: els valors singulars grans amb prou feines es toquen i els petits es domestiquen. Es guanya estabilitat a canvi d'introduir biaix, i triar el paràmetre és una decisió, no un càlcul automàtic.

Mètodes iteratius i autovalors

Per a matrius grans i disperses no es factoritza: s'itera. Amb una descomposició de la matriu en dues parts, la iteració convergeix des de qualsevol punt de partida si i només si el radi espectral de la matriu d'iteració és menor que $1$. Jacobi pren la diagonal; Gauss i Seidel prenen la part triangular inferior inclosa la diagonal. La dominància diagonal estricta és una condició suficient clàssica, no necessària.

El gradient conjugat és el mètode de referència per a matrius simètriques definides positives i grans. En aritmètica exacta acabaria en com a molt $n$ passos, però a la pràctica el que importa és la velocitat de convergència, que depèn del condicionament; per això es fa servir sempre amb precondicionament.

El mètode de la potència normalitza iteracions del producte per la matriu i convergeix a la direcció pròpia dominant, sempre que hi hagi un únic autovalor de mòdul màxim i el vector inicial hi tingui component no nul·la. Si hi ha empat de mòduls, o si has tingut la mala sort de començar en un vector ortogonal a la direcció dominant, la garantia falla.

Exercici resolt: condicionament i Cholesky d'una diagonal

Per a la matriu diagonal amb entrades $2$ i $5$, troba el número de condició en norma euclidiana i la factorització de Cholesky.

Valors singulars. Per a una diagonal amb entrades positives, els valors singulars són les mateixes entrades: $5$ i $2$. Per tant

$$\kappa_2(A) = \frac{5}{2} = 2{,}5.$$

És un condicionament excel·lent: no es perd pràcticament cap xifra.

Cholesky. Busquem una triangular inferior amb diagonal positiva el producte de la qual pel seu transposat doni la matriu. Amb una diagonal és immediat:

$$L = \operatorname{diag}(\sqrt2,\ \sqrt5),$$

perquè el quadrat de cada entrada recupera l'original.

Comprovació de la hipòtesi: la matriu és simètrica i té els dos autovalors positius, o sigui definida positiva. Si una entrada hagués estat negativa, Cholesky hauria fallat en intentar arrelar-la, i això mateix hauria demostrat que no és definida positiva.

Exercici resolt: per què truncar la SVD ajuda amb soroll

Per què una SVD truncada pot millorar una reconstrucció a partir de dades sorolloses?

El mecanisme. La pseudoinversa divideix per cada valor singular. Les direccions amb valor singular molt petit són les que la matriu «aixafa» gairebé del tot; en invertir, es multipliquen per un número enorme, i amb elles el soroll que hi hagi en aquella direcció.

Què fa el truncament. Descartar els valors singulars per sota d'un llindar elimina precisament aquestes direccions. El resultat és molt més estable.

El preu. També descartes la part del senyal real que vivia en aquelles direccions, de manera que introdueixes biaix. És exactament el mateix compromís entre biaix i variància que apareix a l'estadística, i el llindar és una decisió de regularització que s'ha de justificar, no un valor objectivament correcte.

Una manera habitual de triar-lo és buscar un «salt» clar a l'espectre de valors singulars, o fixar-lo a partir de la magnitud estimada del soroll.