Sisteme liniare

Rezolvarea ecuației Ax = b este una dintre cele mai vechi probleme ale matematicii care continuă să ruleze în interiorul rețelei dumneavoastră neuronale.

Tip: Construire Limbaj: Python Cerințe preliminare: Faza 1, lecțiile 01 (Intuiție pentru algebra liniară), 02 (Vectori, matrice și operații), 03 (Transformări matriciale) Durată: ~120 de minute

Obiective de învățare

• Rezolvați Ax = b prin eliminare Gaussiană cu pivotare parțială și substituție inversă
• Factorizați matrice prin descompunerile LU, QR și Cholesky și explicați când este potrivită fiecare
• Deduceți ecuațiile normale pentru problema celor mai mici pătrate și legați-le de regresia liniară și regresia ridge
• Diagnosticați sistemele rău condiționate cu ajutorul numărului de condiționare și aplicați regularizarea pentru a le stabiliza

Problema

De fiecare dată când antrenați o regresie liniară, rezolvați o problemă liniară. De fiecare dată când calculați o ajustare prin metoda celor mai mici pătrate, rezolvați un sistem liniar. De fiecare dată când un strat al unei rețele neuronale calculează y = Wx + b, el evaluează o transformare liniară sau afină ce apare în asemenea sisteme. Când adăugați regularizare, modificați sistemul. Când folosiți procese gaussiene, factorizați o matrice. Când inversați o matrice de covarianță pentru distanța Mahalanobis, în practică ar trebui să rezolvați un sistem liniar.

Ecuația Ax = b apare peste tot. A este o matrice de coeficienți cunoscuți. b este un vector de ieșiri cunoscute. x este vectorul necunoscutelor pe care doriți să îl găsiți. În regresia liniară, A este matricea datelor, b este vectorul țintă, iar x este vectorul ponderilor. Întregul model se reduce la găsirea unui x pentru care Ax este cât mai aproape de b.

În această lecție veți construi de la zero principalele metode de rezolvare a acestei ecuații. Veți înțelege de ce unele metode sunt rapide, iar altele stabile, de ce unele funcționează numai pentru sisteme pătratice, în timp ce altele gestionează sisteme supradeterminate, și cum numărul de condiționare al matricei ajută la cuantificarea sensibilității răspunsului.

Conceptul

Ce înseamnă geometric Ax = b

Un sistem de ecuații liniare are o interpretare geometrică. Fiecare ecuație definește un hiperplan. Soluția este punctul sau mulțimea de puncte în care se intersectează toate hiperplanele.

2x + y = 5          Two lines in 2D.
x - y  = 1          They intersect at x=2, y=1.

Sunt posibile trei situații:

Pentru o matrice pătratică A, o soluție unică pentru orice b înseamnă că A este inversabilă. Absența unei soluții înseamnă că sistemul este incompatibil. Soluțiile multiple apar când sistemul este compatibil, iar A are un nucleu netrivial. Multe probleme de învățare automată au mai multe ecuații, adică puncte de date, decât necunoscute, adică parametri. Un asemenea sistem supradeterminat poate avea sau nu o soluție exactă; atunci când b nu aparține spațiului coloanelor lui A, intervine metoda celor mai mici pătrate.

Notă tehnică a traducerii: Originalul asociază generic „soluții infinite” cu existența unui spațiu nul și afirmă că majoritatea sistemelor cu mai multe ecuații decât necunoscute nu au soluție exactă. Orice matrice are un nucleu, posibil trivial; soluțiile multiple cer un nucleu netrivial și compatibilitatea sistemului. De asemenea, un sistem supradeterminat poate avea o soluție exactă dacă b se află în spațiul coloanelor lui A.

Imaginea coloanelor față de imaginea rândurilor

Există două moduri de a interpreta Ax = b.

Imaginea rândurilor. Fiecare rând al lui A definește o ecuație. Fiecare ecuație este un hiperplan. Soluția este locul în care se intersectează toate hiperplanele.

Imaginea coloanelor. Fiecare coloană a lui A este un vector. Întrebarea devine: ce combinație liniară a coloanelor lui A produce b?

A = | 2  1 |    b = | 5 |
    | 1 -1 |        | 1 |

Imaginea liniilor: rezolvați simultan 2x + y = 5 și x - y = 1.

Imaginea coloanelor: găsiți x1, x2 astfel încât:
  x1 * [2, 1] + x2 * [1, -1] = [5, 1]
  2 * [2, 1] + 1 * [1, -1] = [4+1, 2-1] = [5, 1]   verificare.

Imaginea coloanelor este mai fundamentală. Dacă b se află în spațiul coloanelor lui A, sistemul are o soluție. Dacă nu, căutați punctul cel mai apropiat din spațiul coloanelor. Acest punct cel mai apropiat corespunde soluției prin metoda celor mai mici pătrate.

Eliminarea Gaussiană

Eliminarea Gaussiană transformă Ax = b într-un sistem triunghiular superior Ux = c, pe care îl rezolvați prin substituție inversă. Este cea mai directă metodă.

Algoritmul:

1. Pentru fiecare coloană k (coloana pivotului):
   a. Găsiți elementul cu valoarea absolută maximă din coloana k, pe linia k sau sub ea (pivotare parțială).
   b. Interschimbați linia respectivă cu linia k.
   c. Pentru fiecare linie i de sub k:
      - Calculați multiplicatorul m = A[i][k] / A[k][k]
      - Scădeți din linia i linia k înmulțită cu m.
2. Aplicați substituția inversă: rezolvați începând de la ultima ecuație.

Exemplu:

Original:
| 2  1  1 | 8 |       R2 = R2 - (2)R1     | 2  1   1 |  8 |
| 4  3  3 |20 |  -->  R3 = R3 - (1)R1 --> | 0  1   1 |  4 |
| 2  3  1 |12 |                            | 0  2   0 |  4 |

                       R3 = R3 - (2)R2     | 2  1   1 |  8 |
                                       --> | 0  1   1 |  4 |
                                           | 0  0  -2 | -4 |

Back substitute:
  -2 * x3 = -4    -->  x3 = 2
  x2 + 2  = 4     -->  x2 = 2
  2*x1 + 2 + 2 = 8 --> x1 = 2

Eliminarea Gaussiană necesită O(n^3) operații. Pentru un sistem 1000 × 1000, ordinul de mărime este de aproximativ un miliard de operații în virgulă mobilă. Metoda este rapidă pentru un singur sistem, dar puteți reutiliza mai eficient calculul dacă trebuie să rezolvați mai multe sisteme cu aceeași matrice A.

Pivotarea parțială: de ce contează

Fără pivotare, eliminarea Gaussiană poate eșua sau poate produce rezultate profund inexacte. Dacă un pivot este zero, împărțiți la zero. Dacă este foarte mic, amplificați erorile de rotunjire.

Pivot nepotrivit:                 Cu pivotare parțială:
| 0.001  1 | 1.001 |            Mai întâi, interschimbați liniile:
| 1      1 | 2     |            | 1      1 | 2     |
                                 | 0.001  1 | 1.001 |
m = 1/0.001 = 1000              m = 0.001/1 = 0.001
R2 = R2 - 1000*R1               R2 = R2 - 0.001*R1
| 0.001  1     | 1.001   |      | 1      1     | 2     |
| 0     -999   | -999.0  |      | 0      0.999 | 0.999 |

x2 = 1.000 (corect)             x2 = 1.000 (corect)
x1 = (1.001 - 1)/0.001          x1 = (2 - 1)/1 = 1.000 (corect)
   = 0.001/0.001 = 1.000        Stabil, deoarece multiplicatorul este mic.

În aritmetica în virgulă mobilă cu precizie limitată, varianta fără pivotare poate pierde cifre semnificative. Pivotarea parțială selectează valoarea absolută maximă disponibilă în coloana pivot, la sau sub rândul curent, pentru a limita amplificarea erorilor.

Descompunerea LU

Descompunerea LU factorizează A într-o matrice triunghiulară inferioară L și o matrice triunghiulară superioară U: A = LU. Matricea L stochează multiplicatorii eliminării Gaussiene. Matricea U este rezultatul eliminării.

A = L @ U

| 2  1  1 |   | 1  0  0 |   | 2  1   1 |
| 4  3  3 | = | 2  1  0 | @ | 0  1   1 |
| 2  3  1 |   | 1  2  1 |   | 0  0  -2 |

De ce să factorizați în loc să efectuați doar eliminarea? După ce aveți L și U, rezolvarea lui Ax = b pentru orice vector b nou necesită numai O(n^2) operații:

Ax = b
LUx = b
Let y = Ux:
  Ly = b    (forward substitution, O(n^2))
  Ux = y    (back substitution, O(n^2))

Costul O(n^3) este plătit o singură dată la factorizare. Fiecare rezolvare ulterioară costă O(n^2). Dacă trebuie să rezolvați 1.000 de sisteme cu aceeași matrice A, dar vectori b diferiți, economiile pot fi foarte mari și cresc odată cu dimensiunea matricei.

Notă tehnică a traducerii: Originalul cuantifică economia printr-un „factor de 1000/3”, dar acest factor nu este general. Raportul real depinde de dimensiunea n, de constantele algoritmilor, de numărul de termeni din dreapta și de implementare; afirmația sigură este diferența dintre factorizarea unică O(n^3) și fiecare rezolvare ulterioară O(n^2).

Cu pivotare parțială, obțineți PA = LU, unde P este matricea de permutare care înregistrează schimburile de rânduri.

Descompunerea QR

Descompunerea QR factorizează A într-o matrice cu coloane ortonormale Q și o matrice triunghiulară superioară R: A = QR. Pentru o factorizare QR completă a unei matrice pătratice, Q este ortogonală; în forma economică pentru o matrice dreptunghiulară, coloanele sale sunt ortonormale.

O matrice ortogonală satisface Q^T Q = I. Coloanele sale sunt vectori ortonormali. Înmulțirea cu o matrice ortogonală păstrează lungimile și unghiurile.

A = Q @ R

Q are coloane ortonormale: Q^T Q = I
R este triunghiulară superior

Pentru a rezolva Ax = b:
  QRx = b
  Rx = Q^T b    (înmulțiți cu Q^T; nu este necesară inversarea)
  Aplicați substituția inversă pentru a obține x.

QR este mai stabilă numeric decât formarea ecuațiilor normale pentru problemele celor mai mici pătrate. Procesul Gram–Schmidt construiește Q coloană cu coloană:

Fiind date coloanele a1, a2, ... ale lui A:

q1 = a1 / ||a1||

q2 = a2 - (a2 . q1) * q1        (scădeți proiecția pe q1)
q2 = q2 / ||q2||                (normalizați)

q3 = a3 - (a3 . q1) * q1 - (a3 . q2) * q2
q3 = q3 / ||q3||

R[i][j] = qi . aj    pentru i <= j

Fiecare pas elimină componenta de-a lungul tuturor vectorilor q anteriori, lăsând numai noua direcție ortogonală.

Descompunerea Cholesky

Când A este simetrică, adică A = A^T, și pozitiv definită, adică toate valorile proprii sunt pozitive, o puteți factoriza ca A = L L^T, unde L este triunghiulară inferioară. Aceasta este descompunerea Cholesky.

A = L @ L^T

| 4  2 |   | 2  0 |   | 2  1 |
| 2  5 | = | 1  2 | @ | 0  2 |

L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 pentru k < i))
L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] pentru k < j)) / L[j][j]    pentru i > j

Cholesky necesită aproximativ jumătate din numărul de operații al unei factorizări LU dense și poate economisi spațiu prin exploatarea simetriei. Funcționează numai pentru matrice simetrice pozitiv definite, dar asemenea matrice apar frecvent:

• Matricele de covarianță sunt simetrice pozitiv semidefinite și devin pozitiv definite dacă regularizarea pozitivă elimină toate direcțiile nule.
• Matricele nucleu din procesele gaussiene sunt pozitiv semidefinite; adăugarea unui termen de zgomot strict pozitiv pe diagonală produce de regulă matricea pozitiv definită factorizată în practică.
• Hessiana unei funcții convexe este pozitiv semidefinită; este pozitiv definită numai în condiții suplimentare, de exemplu convexitate strictă locală fără direcții plate.
• A^T A este întotdeauna simetrică pozitiv semidefinită și este pozitiv definită exact când coloanele lui A sunt liniar independente.

Notă tehnică a traducerii: Originalul numește fără condiții matricea nucleu și Hessiana la minim „pozitiv definite”. În general, acestea pot fi doar pozitiv semidefinite. Cholesky standard cere o matrice pozitiv definită, motiv pentru care în procesele gaussiene se factorizează uzual K + σ²I sau se adaugă jitter numeric.

În procesele gaussiene, factorizați matricea nucleu regularizată K prin Cholesky, apoi rezolvați K alpha = y pentru a obține media predictivă. Factorul Cholesky vă oferă și log-determinantul necesar verosimilității marginale: log det(K) = 2 * sum(log(diag(L))).

Pătrate minime: când Ax = b nu are soluție exactă

Dacă A are dimensiunea m × n, cu m > n, sistemul este supradeterminat. O soluție exactă poate exista, dar în datele zgomotoase b nu aparține de obicei exact spațiului coloanelor lui A. În acest caz, minimizați eroarea pătratică:

minimizați ||Ax - b||^2

Aceasta este suma pătratelor reziduurilor:
  sum((A[i,:] @ x - b[i])^2 pentru i din range(m))

Pentru o soluție interioară, minimizatorul satisface ecuațiile normale:

A^T A x = A^T b

Deducere: dezvoltați ||Ax - b||^2 = (Ax - b)^T (Ax - b) = x^T A^T A x - 2 x^T A^T b + b^T b. Calculați gradientul în raport cu x și egalați-l cu zero: 2 A^T A x - 2 A^T b = 0.

Sistemul original (supradeterminat, 4 ecuații, 2 necunoscute):
| 1  1 |         | 3 |
| 1  2 | x     = | 5 |       Niciun x exact nu satisface toate cele 4 ecuații.
| 1  3 |         | 6 |
| 1  4 |         | 8 |

Ecuații normale:
A^T A = | 4  10 |    A^T b = | 22 |
        | 10 30 |            | 63 |

Soluție: x = [1.5, 1.7]

Aceasta este regresie liniară. x[0] este interceptul, iar x[1] este panta.

Rezolvarea corectă a ecuațiilor normale afișate este x = [1.5, 1.6]: x[0] este interceptul, iar x[1] este panta.

Notă tehnică a traducerii: Blocul original indică eronat panta 1.7. Din sistemul 4x₀ + 10x₁ = 22, 10x₀ + 30x₁ = 63 rezultă x₀ = 1.5 și x₁ = 1.6. Blocul a rămas neschimbat pentru fidelitate structurală.

Ecuațiile normale = regresie liniară

Legătura este exactă la nivelul problemei de optimizare. În regresia liniară, matricea datelor X are câte un rând pentru fiecare eșantion și câte o coloană pentru fiecare caracteristică. Vectorul țintă y are câte o intrare pentru fiecare eșantion. Vectorul ponderilor w satisface ecuațiile normale:

X^T X w = X^T y
w = (X^T X)^(-1) X^T y

Aceasta este formula închisă a regresiei liniare când X are rang complet pe coloane, astfel încât X^T X este inversabilă. În practică, apelurile la sklearn.linear_model.LinearRegression.fit() rezolvă aceeași problemă de minimizare prin algoritmi de pătrate minime, fără să fie necesară formarea explicită a ecuațiilor normale.

Notă tehnică a traducerii: Originalul poate sugera că scikit-learn calculează direct (X^T X)^(-1)X^T y. Pentru date dense, LinearRegression folosește scipy.linalg.lstsq; formarea ecuațiilor normale este mai puțin stabilă și formula cu inversă nu se aplică atunci când X nu are rang complet.

Adăugați termenul de regularizare lambda * I la matrice și obțineți regresia ridge:

(X^T X + lambda * I) w = X^T y
w = (X^T X + lambda * I)^(-1) X^T y

Regularizarea îmbunătățește condiționarea matricei și combate supraînvățarea prin micșorarea ponderilor spre zero. Matricea X^T X + lambda * I este întotdeauna simetrică pozitiv definită pentru lambda > 0, astfel încât o puteți rezolva prin Cholesky.

Pseudoinversa Moore–Penrose

Pseudoinversa A+ generalizează inversarea matricelor la matrice nepătratice și singulare. Pentru orice matrice A:

x = A+ b

unde A+ = V Sigma+ U^T    (calculată prin SVD)

Sigma+ se formează prin inversarea fiecărei valori singulare nenule și transpunerea rezultatului. Dacă A = U Sigma V^T, atunci A+ = V Sigma+ U^T.

Notă tehnică a traducerii: În aritmetica exactă se inversează valorile singulare nenule. Numeric, bibliotecile folosesc un prag de toleranță și tratează valorile suficient de mici drept zero, pentru a nu amplifica excesiv erorile de rotunjire.

A = U Sigma V^T        (SVD)

Sigma = | 5  0 |       Sigma+ = | 1/5  0  0 |
        | 0  2 |                | 0  1/2  0 |
        | 0  0 |

A+ = V Sigma+ U^T

Pseudoinversa oferă soluția de normă minimă dintre soluțiile care minimizează eroarea celor mai mici pătrate. Dacă sistemul are:

• O soluție unică: A+ b o furnizează.
• Nicio soluție exactă: A+ b furnizează soluția prin metoda celor mai mici pătrate.
• Infinit de multe soluții: A+ b o furnizează pe cea cu cea mai mică valoare ||x||.

Funcțiile NumPy np.linalg.lstsq și np.linalg.pinv folosesc algoritmi bazați pe descompuneri numerice robuste, iar pinv este definită prin SVD.

Numărul de condiționare

Numărul de condiționare măsoară sensibilitatea soluției la modificări mici ale intrării. Pentru o matrice pătratică inversabilă A, numărul de condiționare indus de o normă este:

kappa(A) = ||A|| * ||A^(-1)|| = sigma_max / sigma_min

Egalitatea cu raportul valorilor singulare maxime și minime este valabilă pentru norma spectrală, adică 2-norma.

Bine condiționat (kappa ~ 1):        Rău condiționat (kappa ~ 10^15):
Modificare mică în b -->             Modificare mică în b -->
modificare mică în x                 modificare uriașă în x

| 2  0 |   kappa = 2/1 = 2          | 1   1          |   kappa ~ 10^15
| 0  1 |   rezolvare sigură         | 1   1+10^(-15) |   soluție inutilizabilă

Reguli orientative:

• Un kappa apropiat de 1 indică o condiționare bună; un prag universal precum 100 nu poate garanta singur acuratețea.
• Pentru kappa ~ 10^k, puteți pierde până la aproximativ k cifre zecimale de precizie, peste erorile proprii algoritmului.
• Pentru kappa apropiat de inversul preciziei mașinii, de exemplu în jur de 10^16 în float64, matricea este efectiv singulară la acea precizie, iar soluția poate fi lipsită de cifre corecte.

Notă tehnică a traducerii: Originalul omite norma din egalitatea cu valorile singulare și declară orice kappa < 100 „sigur”. Raportul sigma_max / sigma_min definește numărul de condiționare în 2-normă, iar acuratețea efectivă depinde și de precizia datelor, perturbații, algoritm și toleranța aplicației.

În învățarea automată, condiționarea slabă apare când caracteristicile sunt aproape coliniare. Pentru sistemul ridge, adăugarea lui lambda * I la X^T X schimbă valorile proprii din sigma_i(X)^2 în sigma_i(X)^2 + lambda, astfel încât numărul de condiționare în 2-normă devine (sigma_max(X)^2 + lambda) / (sigma_min(X)^2 + lambda).

Notă tehnică a traducerii: Formula originală (sigma_max + lambda) / (sigma_min + lambda) este corectă dacă sigma desemnează valorile proprii ale unei matrice simetrice pozitiv semidefinite căreia i se adaugă lambda I. Pentru regularizarea lui X^T X, dacă sigma desemnează valorile singulare ale lui X, acestea trebuie ridicate la pătrat.

Metode iterative: gradientul conjugat

Pentru sisteme rare foarte mari, cu milioane de necunoscute, metodele directe precum LU sau Cholesky sunt prea costisitoare. Metodele iterative aproximează soluția îmbunătățind o estimare pe parcursul mai multor iterații.

Metoda gradientului conjugat (CG) rezolvă Ax = b când A este simetrică pozitiv definită. Găsește soluția exactă în cel mult n iterații în aritmetică exactă, dar în practică poate converge mult mai repede dacă valorile proprii ale lui A sunt grupate favorabil.

Schița algoritmului:
  x0 = aproximația inițială (adesea zero)
  r0 = b - A x0           (reziduu)
  p0 = r0                 (direcție de căutare)

  Pentru k = 0, 1, 2, ...:
    alpha = (rk . rk) / (pk . A pk)
    x_{k+1} = xk + alpha * pk
    r_{k+1} = rk - alpha * A pk
    beta = (r_{k+1} . r_{k+1}) / (rk . rk)
    p_{k+1} = r_{k+1} + beta * pk
    dacă ||r_{k+1}|| < toleranță: opriți

CG este folosit în:

• Optimizare la scară mare, de exemplu metoda Newton-CG
• Rezolvarea discretizărilor ecuațiilor cu derivate parțiale
• Metode cu nuclee în care matricea nucleu este prea mare pentru factorizare
• Sisteme iterative împreună cu precondiționarea

Rata de convergență depinde de distribuția valorilor proprii și poate fi mărginită cu ajutorul numărului de condiționare. Sistemele mai bine condiționate tind să converge mai repede, un alt motiv pentru care regularizarea și precondiționarea ajută.

Imaginea completă: ce metodă și când

Metodă	Cerințe	Cost	Caz de utilizare
Eliminare Gaussiană	A pătratică și nesingulară	O(n^3)	Rezolvarea unică a unui sistem pătratic
Descompunere LU	A pătratică și nesingulară	Factorizare O(n^3) + rezolvare O(n^2)	Rezolvări multiple cu aceeași A
Descompunere QR	Orice A, cu m >= n în cazul de pătrate minime supradeterminat	O(mn^2)	Pătrate minime, stabilitate numerică
Cholesky	A simetrică pozitiv definită	O(n^3/3)	Matrice de covarianță regularizate, procese gaussiene, regresie ridge
Ecuații normale	Sistem supradeterminat, de obicei m > n	O(mn^2 + n^3)	Regresie liniară când n este mic și condiționarea este acceptabilă
SVD / pseudoinversă	Orice A	O(mn^2)	Sisteme cu rang deficient, soluții de normă minimă
Gradient conjugat	A rară, simetrică pozitiv definită	O(k * (nnz + n))	Sisteme rare mari, unde k este numărul de iterații

Notă tehnică a traducerii: Tabelul original indică pentru gradientul conjugat costul O(n * k * nnz). O iterație necesită în principal un produs matrice–vector O(nnz) și operații vectoriale O(n), deci costul total pentru k iterații este O(k(nnz+n)), adesea abreviat O(k·nnz) pentru matrice rare cu cel puțin O(n) elemente nenule.

Legătura cu învățarea automată

Fiecare metodă din această lecție apare în sisteme de învățare automată din producție:

Regresie liniară. Soluția în formă închisă satisface ecuațiile normale X^T X w = X^T y. Problema poate fi rezolvată prin Cholesky dacă numărul caracteristicilor este mic și matricea are rang complet, prin QR dacă stabilitatea numerică este importantă sau prin SVD dacă matricea ar putea avea rang deficient.

Regresie ridge. Adaugă lambda * I la X^T X. Sistemul regularizat (X^T X + lambda * I) w = X^T y poate fi rezolvat prin Cholesky, deoarece X^T X + lambda * I este simetric pozitiv definit pentru lambda > 0.

Procese gaussiene. Media predictivă necesită rezolvarea lui K alpha = y, unde K este matricea nucleu cu termenul de zgomot sau regularizarea corespunzătoare. Factorizarea Cholesky a acestei matrice este abordarea standard. Log-verosimilitatea marginală folosește log det(K) = 2 sum(log(diag(L))).

Inițializarea rețelelor neuronale. Inițializarea ortogonală folosește descompunerea QR pentru a crea matrice de ponderi cu coloane ortonormale. Aceasta ajută la păstrarea scării semnalului în rețele adânci.

Precondiționare. Optimizatoarele la scară mare pot folosi Cholesky incompletă sau LU incompletă drept precondiționatori pentru metode iterative precum gradientul conjugat.

Ingineria caracteristicilor. Numărul de condiționare al lui X sau al lui X^T X indică posibila coliniaritate a caracteristicilor; rețineți că în 2-normă kappa(X^T X) = kappa(X)^2 pentru rang complet. Dacă kappa este mare, puteți elimina caracteristici redundante sau puteți adăuga regularizare.

linear-system-conditioning

Construiți soluția

Pasul 1: eliminare Gaussiană cu pivotare parțială

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

    for k in range(n):
        max_row = k + np.argmax(np.abs(Ab[k:, k]))
        Abk, max_row = Abmax_row, k

        if abs(Ab[k, k]) < 1e-12:
            raise ValueError(f"Matrix is singular or nearly singular at pivot {k}")

        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= m * Ab[k, k:]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]) / Ab[i, i]

    return x

Pasul 2: descompunerea LU

def lu_decompose(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()
    P = np.eye(n)

    for k in range(n):
        max_row = k + np.argmax(np.abs(U[k:, k]))
        if max_row != k:
            Uk, max_row = Umax_row, k
            Pk, max_row = Pmax_row, k
            if k > 0:
                L[[k, max_row], :k] = L[[max_row, k], :k]

        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]

    return P, L, U

def lu_solve(P, L, U, b):
    n = len(b)
    Pb = P @ b.astype(float)

    y = np.zeros(n)
    for i in range(n):
        y[i] = Pb[i] - L[i, :i] @ y[:i]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]

    return x

Notă tehnică a traducerii: Implementarea LU din bloc nu verifică dacă U[k, k] este zero sau aproape zero după pivotare și poate împărți la zero pentru o matrice singulară. Blocul executabil a fost păstrat exact; într-o implementare robustă, adăugați o verificare a pivotului analogă celei din gaussian_elimination.

Pasul 3: descompunerea Cholesky

def cholesky(A):
    n = A.shape[0]
    L = np.zeros_like(A, dtype=float)

    for i in range(n):
        for j in range(i + 1):
            s = A[i, j] - L[i, :j] @ L[j, :j]
            if i == j:
                if s <= 0:
                    raise ValueError("Matrix is not positive definite")
                L[i, j] = np.sqrt(s)
            else:
                L[i, j] = s / L[j, j]

    return L

Pasul 4: pătrate minime prin ecuațiile normale

def least_squares_normal(A, b):
    AtA = A.T @ A
    Atb = A.T @ b
    return gaussian_elimination(AtA, Atb)

def ridge_regression(A, b, lam):
    n = A.shape[1]
    AtA = A.T @ A + lam * np.eye(n)
    Atb = A.T @ b
    L = cholesky(AtA)
    y = np.zeros(n)
    for i in range(n):
        y[i] = (Atb[i] - L[i, :i] @ y[:i]) / L[i, i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - L.T[i, i+1:] @ x[i+1:]) / L.T[i, i]
    return x

Pasul 5: numărul de condiționare

def condition_number(A):
    U, S, Vt = np.linalg.svd(A)
    return S[0] / S[-1]

Folosiți soluția

Reuniți componentele pentru regresia liniară și regresia ridge pe date reale:

np.random.seed(42)
X_raw = np.random.randn(100, 3)
w_true = np.array([2.0, -1.0, 0.5])
y = X_raw @ w_true + np.random.randn(100) * 0.1

X = np.column_stack([np.ones(100), X_raw])

w_ols = least_squares_normal(X, y)
print(f"OLS weights (ours):    {w_ols}")

w_np = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"OLS weights (numpy):   {w_np}")
print(f"Max difference: {np.max(np.abs(w_ols - w_np)):.2e}")

w_ridge = ridge_regression(X, y, lam=1.0)
print(f"Ridge weights (ours):  {w_ridge}")

from sklearn.linear_model import Ridge
ridge_sk = Ridge(alpha=1.0, fit_intercept=False)
ridge_sk.fit(X, y)
print(f"Ridge weights (sklearn): {ridge_sk.coef_}")

Livrați rezultatul

Această lecție produce:

• code/linear_systems.py, care conține implementări de la zero pentru eliminarea Gaussiană, descompunerea LU, descompunerea Cholesky, metoda celor mai mici pătrate și regresia ridge
• O demonstrație funcțională că ecuațiile normale și LinearRegression din scikit-learn produc aceleași ponderi în cazul bine condiționat din exemplu

Exerciții

1. Rezolvați sistemul [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] folosind eliminarea Gaussiană implementată de dumneavoastră, rezolvitorul LU propriu și np.linalg.solve. Verificați că toate trei oferă același răspuns în limita toleranței aritmeticii în virgulă mobilă.

2. Generați o matrice aleatoare X de dimensiune 50 × 5 și ținta y = X @ w_true + noise. Rezolvați pentru w folosind ecuațiile normale, QR prin np.linalg.qr, SVD prin np.linalg.svd și np.linalg.lstsq. Comparați cele patru soluții. Măsurați numărul de condiționare al lui X^T X și explicați cum influențează metoda în care aveți încredere.

3. Creați o matrice aproape singulară făcând două coloane aproape identice, de exemplu column 2 = column 1 + 1e-10 * noise. Calculați-i numărul de condiționare. Rezolvați Ax = b cu și fără regularizare, adăugând 0.01 * I. Comparați soluțiile și reziduurile. Explicați de ce ajută regularizarea.

4. Implementați algoritmul gradientului conjugat pentru o matrice aleatoare simetrică pozitiv definită de dimensiune 100 × 100. Numărați iterațiile necesare pentru a converge la toleranța 1e-8. Comparați cu maximul teoretic de n iterații în aritmetică exactă.

5. Cronometrați rezolvitorul Cholesky propriu, rezolvitorul LU propriu și np.linalg.solve pe matrice simetrice pozitiv definite de dimensiunile 10, 50, 200 și 500. Reprezentați grafic rezultatele. Comparați raportul măsurat cu raportul teoretic de aproximativ 2 dintre numerele de operații ale factorizărilor Cholesky și LU dense; implementările Python și bibliotecile optimizate pot avea rapoarte diferite.

Termeni-cheie

Termen	Cum este numit în conversații	Ce înseamnă de fapt
Sistem liniar	„Rezolvați pentru x”	Un ansamblu de ecuații liniare Ax = b. Găsirea lui x înseamnă găsirea intrării care produce ieșirea b sub transformarea A.
Eliminare Gaussiană	„Reduceți rândurile”	Anularea sistematică a elementelor de sub diagonală prin operații pe rânduri, obținând un sistem triunghiular superior rezolvabil prin substituție inversă. Cost O(n^3).
Pivotare parțială	„Schimbați rândurile pentru stabilitate”	Înaintea eliminării din coloana k, mutați în poziția pivot rândul cu cea mai mare valoare absolută din acea coloană, dintre rândurile eligibile. Reduce riscul amplificării erorilor, dar nu poate face bine condiționată o matrice aproape singulară.
Descompunere LU	„Factorizați în matrice triunghiulare”	Scrieți A = LU, unde L este triunghiulară inferioară și stochează multiplicatorii, iar U este triunghiulară superioară și reprezintă matricea eliminată. Amortizează costul O(n^3) pe mai multe rezolvări.
Descompunere QR	„Factorizare ortogonală”	Scrieți A = QR, unde Q are coloane ortonormale, iar R este triunghiulară superioară. Este mai stabilă decât ecuațiile normale pentru pătrate minime.
Descompunere Cholesky	„Rădăcina pătrată a unei matrice”	Pentru A simetrică pozitiv definită, scrieți A = LL^T. Necesită aproximativ jumătate din operațiile LU. Se folosește pentru matrice de covarianță și nucleu regularizate și pentru regresia ridge.
Pătrate minime	„Cea mai bună ajustare când soluția exactă lipsește”	Minimizează suma pătratelor reziduurilor `
Ecuații normale	„Scurtătura din calculul diferențial”	A^T A x = A^T b, obținute prin anularea gradientului lui `
Pseudoinversă	„Inversare pentru matrice nepătratice”	A+ = V Sigma+ U^T prin SVD. Oferă soluția de normă minimă dintre soluțiile de pătrate minime pentru orice matrice, pătratică sau dreptunghiulară, singulară sau nesingulară.
Număr de condiționare	„Câtă încredere puteți avea în răspuns?”	În 2-normă, kappa = sigma_max / sigma_min. Măsoară sensibilitatea la perturbațiile intrării. Puteți pierde până la aproximativ log10(kappa) cifre de precizie.
Regresie ridge	„Pătrate minime regularizate”	Rezolvați (X^T X + lambda I) w = X^T y. Adăugarea lui lambda I îmbunătățește condiționarea și micșorează ponderile spre zero, reducând riscul de supraînvățare.
Gradient conjugat	„Ax=b iterativ pentru matrice mari”	Un rezolvitor iterativ pentru sisteme simetrice pozitiv definite. Converge în cel mult n pași în aritmetică exactă. Este practic pentru sisteme rare mari, unde factorizarea este prea costisitoare.
Sistem supradeterminat	„Mai multe date decât parametri”	m > n într-un sistem m × n. Poate avea o soluție exactă dacă b aparține spațiului coloanelor; altfel, metoda celor mai mici pătrate găsește cea mai bună aproximație. Multe probleme de regresie au această formă.
Substituție inversă	„Rezolvați de jos în sus”	Pentru un sistem triunghiular superior, rezolvați mai întâi ultima ecuație, apoi substituiți în sens invers. Cost O(n^2).
Substituție directă	„Rezolvați de sus în jos”	Pentru un sistem triunghiular inferior, rezolvați mai întâi prima ecuație, apoi substituiți înainte. Cost O(n^2). Este folosită în pasul L al rezolvărilor LU.

Lecturi suplimentare

• MIT 18.06: Linear Algebra (Gilbert Strang) — cursul de referință despre sisteme liniare și factorizări matriciale
• Numerical Linear Algebra (Trefethen și Bau) — referința standard pentru stabilitate numerică, condiționare și cauzele eșecurilor algoritmilor
• Matrix Computations (Golub și Van Loan) — referința enciclopedică pentru algoritmii matriciali
• 3Blue1Brown: Inverse Matrices — intuiție vizuală pentru sensul geometric al rezolvării lui Ax = b

---

Sursă: Linear Systems — original

Navigare: înapoi: 01.16 — Metode de eșantionare · Faza 1 — Fundamente matematice · Catalog complet · în continuare: 01.18 — Optimizare convexă.